Anaïs Crestetto
email: anais.crestetto@univ-nantes.fr

Nicolas Crouseilles
email: nicolas.crouseilles@inria.fr

Mohammed Lemou
email: mohammed.lemou@univ-rennes1.fr

A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling

published or not. The documents may come

A particle micro-macro decomposition based numerical scheme for collisional kinetic equations in the diffusion scaling

Anaïs Crestetto, Nicolas Crouseilles, Mohammed Lemou

Introduction

Particle systems appearing in plasma physics or radiative transfer can be described at different scales. When the system is far from its thermodynamical equilibrium, a kinetic description is necessary. Particles are then represented by a distribution function f which depends on time t ≥ 0, position x ∈ R d and velocity v ∈ V ⊂ R d , d ≥ 1. The distribution f (t, x, v) satisfies a collisional kinetic equation. Particle methods are often used for simulating kinetic problems, especially in realistic 3-dimensional situations, d = 3. However, they are affected by numerical noise due to their probabilistic character. A simple way to reduce this noise is to increase the number of particles, but then the numerical cost increases as well. Other standard kinetic descriptions, as phase space grid methods, may require too much memory in the two or three dimensional framework. Otherwise, macroscopic descriptions depending only on t and x can be sufficient if the system stays near its thermodynamical equilibrium, and are less expensive since their unknown does not depend on the velocity variable anymore. Beside the noisy character of standard particle methods, there is an additional difficulty in kinetic descriptions which is linked to the presence of various scales in the system. Multi-scale phenomena may indeed appear in plasma devices or radiative transfer applications, depending on some physical parameters as for example the mean free path of particles or the Knudsen number denoted here by ε. This multi-scale character is often represented by stiff terms in the kinetic equation, and the general challenge is to construct efficient numerical methods for these multiscale kinetic equations: this means that, without numerically resolving the stiffness, the numerical method must solve accurately the kinetic regime, must have the right asymptotics in the high-stiffness limit (the so-called asymptotic preserving property) and its computational cost should decrease as the system approaches the equilibrium (a time diminishing property). Note that direct numerical methods whose parameters resolve the smallest scale of size ε are impossible to use, since they automatically involve an extremely high computational cost.

Several strategies have been proposed to overcome this strong constraint. Domain decomposition methods can be applied when we have different regions with different values of the scaling parameter, see [START_REF] Degond | Fluid simulations with localized Boltzmann upscaling by Direct Simulation Monte-Carlo[END_REF][START_REF] Golse | A domain decomposition analysis for a two-scale linear transport problem[END_REF][START_REF] Tallec | Coupling Boltzmann and Navier-Stokes by half-fluxes[END_REF][START_REF] Tiwari | A particle-particle hybrid method for kinetic and continuous equations[END_REF]. When the different scales are less clearly delimited, we have to develop kinetic schemes that naturally reduce to good approximations of the macroscopic problem when the system goes near its equilibrium, and overcome the stiffness. Such schemes are often called Asymptotic-Preserving (AP), see [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF][START_REF] Krycki | Asymptotic preserving numerical schemes for a nonclassical radiation transport model for atmospheric clouds[END_REF][START_REF] Klar | An Asymptotic-Induced Scheme For Nonstationary Transport Equations In The Diffusive Limit[END_REF][START_REF] Dimarco | Implicit-Explicit Runge-Kutta schemes for the Boltzmann-Poisson system for semiconductors[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF][START_REF] Lemou | Relaxed micro-macro schemes for kinetic equations[END_REF][START_REF] Lemou | Micro-Macro Schemes for Kinetic Equations Including Boundary Layers[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Jin | Uniformly accurate diffusive relaxation schemes for multiscale transport equations[END_REF][START_REF] Naldi | Numerical schemes for kinetic equations in diffusive regimes[END_REF][START_REF] Buet | Asymptotic preserving scheme and numerical methods for radiative hydrodynamic models[END_REF]]. Mainly, the numerical cost remains comparable to the one of the non-stiff kinetic problem, even when ε ≪ 1.

Our goal is to design an efficient AP scheme, using particles, for the following kinetic radiative transport equation (RTE) in the diffusion scaling

∂ t f + 1 ε v∂ x f = 1 ε 2 (ρM -f ), f (t = 0, x, v) = f 0 (x, v), (1.1)
where

x ∈ Ω ⊂ R, ρ(t, x) = 1 2 V f (t, x, v)dv, V = [-1, 1], M (v) = 1 ∀v ∈ V and f 0 (x, v
) is a given initial condition. Periodic boundary conditions are considered. It is well-known (see [START_REF] Larsen | Asymptotic solution of neutron transport problems for small mean free paths[END_REF][START_REF] Degond | Diffusion limit for non homogeneous and non-microreversible processes[END_REF]) that when ε goes to zero, the distribution function f (t, x, v) converges towards ρ(t, x)M (v), where ρ satisfies the following diffusion equation

∂ t ρ - 1 3 ∂ xx ρ = 0, ρ(t = 0, x) = 1 2 V f 0 (x, v)dv. (1.2)
An extension to the Vlasov-Poisson-BGK case is presented in Section 5. The kinetic equation is coupled to a Poisson equation for the electric field denoted by E(t, x). More precisely, we consider

∂ t f + 1 ε v∂ x f + 1 ε E∂ v f = 1 ε 2 (ρM -f ), (1.3)
∂ x E = ρ -1, (1.4)
where

x ∈ Ω ⊂ R, ρ(t, x) = V f (t, x, v)dv, V = R, M (v) = 1 √ 2π exp -v 2 2
is the absolute Maxwellian and we consider periodic boundary conditions. Note that an additional condition Ω Edx = 0 is imposed to obtain a well-posed problem. When ε goes to zero, the asymptotic model is a drift-diffusion equation satisfied by ρ(t, x) (see [START_REF] Ben Abdallah | Diffusion Approximation for the one dimensional Boltzmann-Poisson system[END_REF])

∂ t ρ -∂ x (∂ x ρ -Ē ρ) = 0, ρ(t = 0, x) = R f 0 (x, v)dv, (1.5)
where Ē is linked to ρ by the Poisson equation ∂ x Ē = ρ -1. The strategy will be the use of the micro-macro decomposition (see [START_REF] Liu | Boltzmann Equation: Micro-Macro Decompositions and Positivity of Shock Profiles[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF]). It consists in writing the distribution function as the sum of the equilibrium part ρ(t, x)M (v) and a rest g(t, x, v). One can then derive a system of two equations: a kinetic one for the rest g(t, x, v) and a macroscopic one for the equilibrium ρ(t, x)M (v). AP micro-macro schemes for (1.1) have been proposed in [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF]. These schemes consist in a semi-implicit phase space grid method for the kinetic part, coupled to a classical spatial grid method for the macro part. Our strategy in this work follows the strategy of [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF] in the case of a fluid scaling: we use particles to sample the kinetic part whereas an Eulerian solver is used to discretize the macro unknown. The main motivation of this strategy lies in the fact that the micro part g converges to zero when ε goes to zero, so that very few particles can sample it. As a consequence, in this regime, the cost of the global micro-macro solver is almost the same as the cost of an asymptotic solver for (1.2).

In this work, we focus on a diffusion type scaling (as in (1.1) or (1.3)) so that an additional scale is involved compared to the fluid scaling considered in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF]. In [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF], a diffusion scaling was studied, but using a fully grid based solver. Hence, the stiffest term (of order 1/ε 2 ) is considered implicit in time in the micro equation, which enables to stabilize the transport term (of order 1/ε) and then to derive an AP scheme for (1.1) and (1.3). The use of particles for the micro part prevents from a similar strategy since a splitting between the transport term (of order 1/ε) and the source term needs to be done. Then, a uniformly stable scheme is hard to obtain in this context. To overcome this difficulty, a suitable formulation of the original model (1.1) is performed, so that the stiff transport term (1/ε) v∂ x g becomes (ε/∆t)(1 -e -∆t/ε 2 )v∂ x g where ∆t > 0 denotes a fixed time step of a numerical time discretization, which will be used to solve this equation (see [START_REF] Crestetto | Asymptotic-Preserving scheme based on a Finite Volume/Particle-In-Cell coupling for Boltzmann-BGK-like equations in the diffusion scaling, Finite Volumes for Complex Applications VII -Elliptic, Parabolic and Hyperbolic Problems[END_REF]). This reformulation is correct up to ∆t 2 (for fixed ε > 0) and ensures that the transport speed remains finite even when ε → 0. This formulation is the starting point of the design of micromacro-particle based numerical schemes which enjoy the AP property and for which the numerical cost diminishes as ε goes to zero. This approach is extended to the second-order (in time) and to the Vlasov-Poisson-BGK case (1.3)- (1.4).

The sequel of the paper is organized as follows. In Section 2, we recall the formal derivation of the asymptotic model of (1.3) and (1.1). The first-order (in time) reformulation of (1.1) is presented in Subsection 3.1 and its Lagrangian discretization in Subsection 3.2. Its extension to a second-order in time model is detailed in Section 4: the continuous model is presented in Subsection 4.1 and its discretization is developed in Subsection 4.2. Section 5 proposes an extension of our strategy to the Vlasov-Poisson-BGK system. Finally, Section 6 is devoted to numerical simulations.

Diffusion asymptotics

In this section, we recall the main steps of the derivation of the model obtained from (1.1) when ε goes to zero. To do so, we consider the micro-macro decomposition (see [START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF][START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Liu | Boltzmann Equation: Micro-Macro Decompositions and Positivity of Shock Profiles[END_REF]) of f :

f (t, x, v) = ρ(t, x)M (v) + g(t, x, v), with ρ(t, x) = f , M (v) = 1 and the rest g satisfies g = 0. Here f = 1 2 V f (v)dv, with V = [-1, 1]. The following micro-macro model is equivalent to the original model (1.1)      ∂ t ρ + 1 ε ∂ x vg = 0, ∂ t g + 1 ε (I -• ) [v∂ x (ρM + g)] = - 1 ε 2 g. (2.1)
Since vM = 0, the micro equation can be rewritten as

∂ t g + 1 ε vM ∂ x ρ + 1 ε (I -• )(v∂ x g) = - 1 ε 2 g. (2.2)
When ε goes to zero, one gets from (2.2), g = -ε(vM ∂ x ρ) + O(ε 2 ), which gives using v 2 M = 1/3 the following diffusion equation satisfied by the limit ρ

∂ t ρ - 1 3 ∂ xx ρ = 0.
The same calculations enable to derive the micro-macro model equivalent to

(1.3)      ∂ t ρ + 1 ε ∂ x vg = 0, ∂ t g + 1 ε (I -• ) [v∂ x (ρM + g) + E∂ v (ρM + g)] = - 1 ε 2 g, (2.3)
where M (v) is now the Maxwellian equilibrium and f = V f (v)dv, with V = R. When ε goes to zero, one gets g = -ε(vM ∂ x ρ -vM Eρ) + O(ε 2 ), which gives, using v 2 M = 1 the following drift-diffusion equation satisfied by the limit ρ

∂ t ρ -∂ xx ρ + ∂ x ( Ē ρ) = 0, ∂ x Ē = ρ -1.
(2.4)

3 First-order in time reformulation and its discretization

In this part, a first-order reformulation of the micro part is proposed, which enables to avoid the stiff transport term in space. The strategy is presented in the case of the equation (1.1) and its corresponding micro-macro model (2.1).

3.1 First-order in time reformulation.

We start with (1.1) (with periodic boundary condition in space) and consider the micro-macro decomposition of f = ρ + g (here

M (v) = 1 for all v ∈ [-1, 1]
) and the micro-macro model (2.1). First, using the relation ∂ t (e t/ε 2 g) = e t/ε 2 (g/ε 2 + ∂ t g), we rewrite the micro part (2.2) as

∂ t (e t/ε 2 g) = - e t/ε 2 ε F (ρ, g) , (3.1)
where F (ρ, g) is given by

F (ρ, g) = v∂ x ρ + v∂ x g -∂ x vg . (3.2)
We denote ∆t > 0 the time step, t n = n∆t with n ∈ N. Then, a second stage consists in integrating (3.1) on [t n , t n+1 ] to get

g(t n+1 ) = e -∆t/ε 2 g(t n ) -ε(1 -e -∆t/ε 2 )F (ρ(t n ), g(t n )) + O(∆t 2 ).
To derive a continuous (in time) equation, we make appear a discrete time derivative on the lefthand side

g(t n+1 ) -g(t n ) ∆t = - 1 -e -∆t/ε 2 ∆t g(t n ) -ε 1 -e -∆t/ε 2 ∆t F (ρ(t n ), g(t n )) + O(∆t), (3.3)
which can be rewritten, up to terms of order O(∆t), as

∂ t g(t n ) = - 1 -e -∆t/ε 2 ∆t g(t n ) -ε 1 -e -∆t/ε 2 ∆t F (ρ(t n ), g(t n )) , ∀n.
We finally obtain the first-order reformulation of (1.1)

∂ t ρ + 1 ε ∂ x vg = 0, (3.4)
∂ t g = - 1 -e -∆t/ε 2 ∆t g -ε 1 -e -∆t/ε 2 ∆t F (ρ, g) , (3.5)
with F (ρ, g) given by (3.2). We remark that the micro equation does not contain any stiff term and satisfies the following property: for all fixed ε > 0, equation (3.5) is consistent with the initial micro equation (2.2) as ∆t goes to zero. Then, it has a suitable form for a numerical discretization using either a deterministic scheme (described in Appendix A) or a particle scheme. Note that the main goal of the paper is to present an AP particle scheme, the cost of which decreases when ε → 0, but in the numerical tests, we will use the deterministic AP scheme based on (3.5) for comparison.

Lagrangian discretization.

This subsection is devoted to the derivation of an AP-particle based numerical scheme for (3.4)- (3.5).

We propose now a Lagrangian discretization of (3.5). More precisely, we adopt a weighted particle method, see [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF], and consider a set of N p ∈ N macro particles. The position of particle k,

1 ≤ k ≤ N p , is denoted by x k (t) ∈ Ω = [0, L x ], with L x > 0, its velocity by v k (t) ∈ V = [-1, 1] and its weight by ω k (t) ∈ R. Let L v = |V | = 2.
The function g is then assumed to be of the form

g (t, x, v) ≈ Np k=1 ω k (t)δ(x -x k (t))δ(v -v k (t)), (3.6)
where δ denotes the Dirac mass function. Weights ω k (t) are related to the distribution function g through

ω k (t) = g(t, x k (t), v k (t)) L x L v N p . (3.7)
Initially, particles are uniformly distributed in the phase-space domain [0, L x ] × V and their weights are computed following (3.7). Note that another approach can also be chosen; for example (unweighted) Monte Carlo methods can also be used (see [START_REF] Crouseilles | Asymptotically stable and time diminishing schemes for rarefied gas dynamics[END_REF][START_REF] Crestetto | Micro-macro decomposition based numerical scheme using Monte-Carlo technique for collisional kinetic equations in the diffusion scaling[END_REF]) and in this case, collisions are taken into account through the change of particles velocities, whereas in our weighted particle method collisions are taken into account through the variation of the weights ω k .

The density ρ is computed on a uniform spatial grid defined by x i = i∆x, i = 0, . . . , N x , N x ∈ N ⋆ and ∆x = L x /N x . We denote by ρ n i the approximation at time t n = n∆t and position x i of ρ(t n , x i ), with ∆t > 0 the time step. Moreover,

g n (x, v) ≈ g(t n , x, v), x n k ≈ x k (t n ), v n k ≈ v k (t n ) and w n k ≈ w k (t n ).
Let us remark that vk = 0, so that the velocities v k (t) are constant in time and we will note v n k = v 0 k =: v k for all n. Thus, the unknowns of the method are (ρ n i ), ∀i = 0, . . . , N x , n > 0, (x n k ) and (ω n k ), ∀k = 1, . . . , N p , n > 0. Our goal is then to extend the particle discretization of [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF] to diffusion scaling. Whereas the hydrodynamic scaling is considered in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF], we have here an additional scale of order 1/ε in front of the transport term. This scale is difficult to handle with a particle method, since it can not be stabilized by the collision term of order 1/ε 2 . Moreover, a specific treatment of the macro flux ∂ x vg is required in the diffusion scaling to ensure the AP property, as we will see later. To deal with the additional scale, we exploit the reformulation (3.5). As already said in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF], we have to use a splitting procedure between the transport part and the source part. Then, the (first-order) splitting writes

• start with an initial repartition of the N p particles (x 0 k , v 0 k ), with ω 0 k = g(t = 0, x 0 k , v 0 k )L x L v /N p ,
• solve the transport part

∂ t g = ε 1 -e -∆t/ε 2 ∆t v∂ x g, with the (non stiff) characteristics ẋk = ε 1 -e -∆t/ε 2 ∆t v k , (3.8)
• solve the source part

∂ t g = - 1 -e -∆t/ε 2 ∆t g -ε 1 -e -∆t/ε 2 ∆t [v∂ x ρ -∂ x vg ] ,
using the equation satisfied by the weights

ωk = - 1 -e -∆t/ε 2 ∆t ω k -ε 1 -e -∆t/ε 2 ∆t [v k ∂ x ρ(x k ) -∂ x vg (x k )] L x L v N p .
(3.9)

Remark 3.1 (Preservation of the micro-macro structure.). As detailed in subsection 4.2 in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF], we have to correct the particle weights in order to preserve at the numerical level the micro-macro structure. Indeed, the micro-macro decomposition technique uses the zero-mean property g = 0.

This property is preserved at the continuous level by the couple (3.8)-(3.9). But the splitting breaks the operator (I -• ) and so this property. Thus, nothing guarantees that this property is satisfied at the discrete level (on the weights ω k ). That is why we have to correct the weights, by applying a discrete approximation of the operator (I -• ) to the weights (ω k ) k , which is consistent with the continuous model. In the simpler case (that is by using a regularization of order ℓ = 0, see (3.12) below), it consists in the following correction:

∀k ∈ I i := k / x k ∈ [x i-1/2 , x i+1/2 ] , ω k ← ω k -∆x g (x i ) k∈I i p k p k ,
where

p k := ρ(x k )M (v k )L x L v /N p
is the weight associated to the Maxwellian and g (x i ) is computed in the same way as vg (x i ) in (3.11).

For more details of this correction (called projection step in following algorithms), we refer the reader to subsection 4.2 in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF]. Now, we detail the time discretization of the two steps. First, (3.8) is approximated by a simple forward Euler scheme

x n+1 k = x n k + ε(1 -e -∆t/ε 2 )v k . (3.10)
Second, we compute the last term in (3.9). The term vg is approximated on the spatial grid x i using

vg (x i ) ≈ Np k=1 ω n k B ℓ (x i -x n+1 k )v k , (3.11)
where B ℓ ≥ 0 is a B-spline function of order ℓ:

B ℓ (x) = (B 0 * B ℓ-1 )(x), with B 0 (x) = 1 ∆x if |x| < ∆x/2, 0 else. (3.12)
We then approximate the equation on the weights (3.9) using a first-order explicit integrator

ω n+1 k = e -∆t/ε 2 ω n k -ε(1 -e -∆t/ε 2 ) [α n k + β n k ] , (3.13)
with

α n k = v k ∂ x ρ n (x n+1 k ) L x L v N p and β n k = -∂ x vg (x n+1 k ) L x L v N p . (3.14)
To compute α n k (resp. β n k ), since ρ n (resp. vg ) is known on the spatial grid, we approximate ∂ x ρ n (resp. ∂ x vg ) by centered finite differences and evaluate at x n+1 k using an interpolation with B-spline functions, for example

∂ x ρ n (x n+1 k ) ≈ Nx i=1 ρ n i+1 -ρ n i-1 2∆x B ℓ (x i -x n+1 k ).
Let us remark that in the limit ε → 0, the particles do not move anymore (see (3.10)) and their weights ω k tend to zero (see (3.13)), as well as ωk .

The macro equation (3.4) is advanced through

ρ n+1 i = ρ n i - ∆t ε vg n+1 i+1 -vg n+1 i-1 2∆x , (3.15)
where vg n+1 i is computed using (3.11). However, this discretization of the macro equation does not produce a time diminishing AP scheme, since it is not accurate in the limits ε → 0 and N p → 0. Indeed, the error due to the particle approximation of g n+1 is amplified by the factor 1/ε. To ensure this time diminishing property, we perform a decomposition of (3.13) so that the macro flux becomes

vg n+1 i = Np k=1 ω n+1 k B ℓ (x i -x n+1 k )v k , = Np k=1 e -∆t/ε 2 ω n k -ε(1 -e -∆t/ε 2 )β n k B ℓ (x i -x n+1 k )v k -ε(1 -e -∆t/ε 2 ) Np k=1 α n k B ℓ (x i -x n+1 k )v k , = -ε(1 -e -∆t/ε 2 ) Np k=1 α n k B ℓ (x i -x n+1 k )v k + h n i , (3.16)
with B ℓ given by (3.12) and

h n i = e -∆t/ε 2 Np k=1 ω n k B ℓ (x i -x n+1 k )v k -ε(1 -e -∆t/ε 2 ) Np k=1 β n k B ℓ (x i -x n+1 k )v k . (3.17)
Since α n k is the weight of v∂ x ρ n (see (3.14)) and v 2 = 1/3, we can write

vg n+1 i ≈ -ε(1 -e -∆t/ε 2 ) 1 3 ∂ x ρ n i + h n i ,
so that the macro scheme becomes

ρ n+1 i = ρ n i + ∆t(1 -e -∆t/ε 2 ) 1 3 ρ n+1/2 i+1 -2ρ n+1/2 i + ρ n+1/2 i-1 ∆x 2 - ∆t ε h n i+1 -h n i-1 2∆x , (3.18)
where ρ n+1/2 can be chosen equal to ρ n or to ρ n+1 depending on the desired asymptotic scheme (explicit or implicit in time). Obviously, the choice of an implicit scheme (ρ n+1/2 = ρ n+1 in (3.18)) enables to get rid of the diffusion-type constraint: ∆t = O(∆x 2 ). We can write the following proposition.

Proposition 3.1. The scheme given by (3.10)-(3.13)-(3.18) enjoys the AP property, i.e. it satisfies the following properties

• for fixed ε > 0, the scheme is a first-order (in time) approximation of the original model (1.1),

• for fixed ∆t > 0, the scheme degenerates into an explicit or implicit first-order scheme of (1.2) (according to the choice of ρ n+1/2 ).

Proof. The consistency follows directly from standard approximation. For the asymptotic behavior, when ε goes to zero, we get immediately from (3.13) ω n+1 k = O(ε) ∀n ≥ 0. From (3.17), we deduce that h n i = O(ε 2 ) ∀n ≥ 1. The macro equation (3.18) then reduces when ε → 0 to a consistent discretization of (1.2).

The scheme is finally summarized in the following algorithm. Algorithm 3.1.

• Initialize (x 0 k , v 0 k ), ω 0 k and ρ 0 i . At each time step: • 1) Advance micro part:
advance the characteristics with (3.10),

-compute vg with (3.11),
advance the equation on the weights with (3.13).

• 2) Projection step: compute (I -• )g n+1 using [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF].

• 3) Advance macro part:

compute h n i with (3.17), -compute ρ n+1 with (3.18).

Second-order in time reformulation and its discretization

This section is devoted to the derivation of a second-order scheme for the micro-macro system (2.1). As for the first-order scheme, we will first reformulate the microscopic equation (2.2) in order to suppress stiff terms (see Subsection 3.1 for the first-order case), and then discretize the obtained micro-macro model to get an AP efficient numerical scheme (see Subsection 3.2 for the first-order case).

Second-order reformulation.

Let us start from the following (equivalent) reformulation of the micro part of (2.1)

∂ t e t/ε 2 g = - e t/ε 2 ε F (ρ(t), g(t)) ,
where F (ρ, g) is defined by (3.2). We now integrate with respect to t ∈ [t n , t n+1 ] and use a second-order mid-point quadrature

g(t n+1 ) = e -∆t/ε 2 g(t n ) - ∆te -∆t/2ε 2 ε F ρ(t n+1/2 ), g(t n+1/2 ) + O ∆t 3 .
To derive a continuous (in time) equation, we make appear a discrete time derivative on the lefthand side

g(t n+1 ) -g(t n ) ∆t = e -∆t/ε 2 -1 ∆t g(t n ) - e -∆t/2ε 2 ε F ρ(t n+1/2 ), g(t n+1/2 ) + O ∆t 2 .
We now look for a continuous (in time) equation for which the previous relation is a second-order numerical scheme. To do so, we perform Taylor expansions of the different terms at t n+1/2

∂ t g(t n+1/2 ) = e -∆t/ε 2 -1 ∆t g(t n+1/2 ) - ∆t 2 ∂ t g(t n+1/2 ) - e -∆t/2ε 2 ε F ρ(t n+1/2 ), g(t n+1/2 ) + O ∆t 2 . (4.1)
Finally, the microscopic equation of (2.1) is reformulated up to the second-order by

∂ t g = 2 ∆t e -∆t/ε 2 -1 e -∆t/ε 2 + 1 g - 2 ε e -∆t/2ε 2 e -∆t/ε 2 + 1 F (ρ, g) ,
and we can now consider the second-order reformulated micro-macro system

∂ t ρ + 1 ε ∂ x vg = 0, (4.2)
∂ t g = 2 ∆t e -∆t/ε 2 -1 e -∆t/ε 2 + 1 g - 2 ε e -∆t/2ε 2 e -∆t/ε 2 + 1 [v∂ x ρ + v∂ x g -∂ x vg ] . (4.3)

Time discretization.

We are now interested in the construction of an AP scheme for system (4.2)-( 4.3), based on a second-order splitting method for the time discretization and a Lagrangian method for the phase space discretization of the micro part.

In the sequel, we will use the same notations as in Subsection 3.2. The splitting method is based on a prediction step on ∆t/2 (first-order) and a correction step on ∆t. Then, a second-order (in time) scheme for (4.2)-(4.3) would read Prediction step on ∆t/2

g n+1/2 =g n + e -∆t/ε 2 -1 e -∆t/ε 2 + 1 g n - ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 F (ρ n , g n ) , (4.4)
ρ n+1/2 =ρ n - ∆t 2ε ∂ x vg n+1/2 , (4.5)
Correction step on ∆t

g n+1 =g n + 2 e -∆t/ε 2 -1 e -∆t/ε 2 + 1 g - 2∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 F ρ n+1/2 , g n+1/2 , (4.6)
ρ n+1 =ρ n - ∆t ε ∂ x vg n+1/2 . (4.7)
We still have to fix g in (4.6) in order to get a second-order scheme and to ensure the convergence of g n+1 to zero as ε goes to zero. It turns out the choice g = g n +g n+1

2 ensures the two conditions. Indeed, we get for the correction step of the micro part

g n+1 = e -∆t/ε 2 g n - ∆t ε e -∆t/2ε 2 F ρ n+1/2 , g n+1/2 .
Note that g = g n+1/2 ensures the second-order accuracy of the scheme but does not ensure the convergence of g n+1 to zero as ε goes to zero (it gives g n+1 = g n ).

Up to now, the micro part converges exponentially fast to zero (when ε goes to zero), so that the asymptotic behavior of the scheme is ρ n+1 = ρ n . Hence, the last but not the least step consists in modifying the macro flux ∂ x vg n+1/2 in (4.7) to capture the correct asymptotic limit.

As done in Section 3, we will modify the discretization of the macro flux in order to make appear the diffusion term directly in the macro part (4.7). This allows us to take the diffusion term implicit and thus to avoid a constraint of diffusion-type ∆t = O(∆x 2 ) in the limit ε → 0 (as done in Section 3). However, the modification has to be of order ∆t 3 for a fixed ε > 0 to not spoil the second-order accuracy of the scheme. The correction we propose consists in computing ρ n+1 i with (4.7), in which we add a diffusion term (discretized with a second-order Crank-Nicolson scheme) in front of which we added the coefficient ∆t(1 -e -∆t/ε 2 ) 2 . This coefficient is of order O(∆t 3 ) for fixed ε > 0 and degenerates to ∆t when ε goes to zero as required. We then obtain

ρ n+1 i = ρ n i - ∆t ε ∂ x vg n+1/2 i + ∆t(1 -e -∆t/ε 2 ) 2 1 3 ∂ xx ρ n+1 i + ρ n i 2 . (4.8)
Note that from (4.4), we see that g n+1/2 behaves like e -∆t/2ε 2 when ε → 0, and so the asymptotic model derived from (4.4)-(4.7) reduces to ρ n+1 = ρ n , which is not the desired diffusion limit.

Therefore we added the term ∆t(1 -e -∆t/ε 2 )

2 1 3 ∂ xx ρ n+1 i +ρ n i 2
in (4.8) in order to capture the right asymptotic regime. This additional term is of order ∆t 3 , so we keep the second-order accuracy for all fixed ε > 0. However when ε → 0, this term is no longer negligible and provides in this limit the right diffusion asymptotic limit. Of course other choices than the multiplicative factor (1-e -∆t/ε 2 ) 2 are possible. For example, we could have taken 1 -

1 1+∆t/ √ ε 2
, which gives second-order accuracy for all fixed ε > 0 and the right asymptotic limit too.

Lagrangian discretization.

We consider the same notations as in Subsection 3.2 and detail here the Lagrangian discretization of the micro-macro system (4.2)-(4.3). In the prediction step (4.4), we compute x n+1/2 k with a forward Euler integrator

x n+1/2 k = x n k + ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 v k , (4.9)
and advance the weights with

w n+1/2 k = 2e -∆t/ε 2 e -∆t/ε 2 + 1 w n k - ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 [v k ∂ x ρ n (x n k ) -∂ x v k g n (x n k ) ] L x L v N p . (4.10)
We end this prediction step by computing the flux vg n+1/2 with (3.11) to get the density

ρ n+1/2 i = ρ n i - ∆t 2ε ∂ x vg n+1/2 i . (4.11)
Now in the correction step, we compute the position at t n+1 with

x n+1 k = x n k + ∆t ε e -∆t/2ε 2 v k , (4.12)
then the weights are given by

w n+1 k = e -∆t/ε 2 w n k - ∆t ε e -∆t/2ε 2 v k ∂ x ρ n+1/2 (x n+1/2 k ) -∂ x v k g n+1/2 (x n+1/2 k ) L x L v N p . (4.13)
Now, using (4.8) in the last step, we have

ρ n+1 i = ρ n i - ∆t ε ∂ x vg n+1/2 i + ∆t(1 -e -∆t/ε 2 ) 2 1 3 ∂ xx ρ n+1 i + ρ n i 2 , (4.14)
where vg n+1/2 i is computed using (3.11). The previous scheme is a semi-discrete and second-order approximation in time of the original equation (2.1). To get a fully discretized scheme, the space derivatives are approximated exactly in the same way as in Subsection 3.2. We finally have the following result. • for fixed ε > 0, the scheme is a second-order (in time) approximation of the original model (1.1),

• for fixed ∆t > 0, the scheme degenerates into an implicit second-order (in time) scheme of (1.2).

Proof. When ε → 0, we get from (4.10) w n+1/2 k → 0 exponentially fast and then vg n+1/2 i → 0.

By injecting it in the macro equation (4.14), we have at the limit

ρ n+1 i = ρ n i + ∆t 3 ∂ xx ρ n+1 i +ρ n i 2
, which is a Crank-Nicolson discretization of the diffusion equation (1.2).

Remark 4.1. Let us emphasize that the moments • have to be computed with B-spline functions of order ℓ ≥ 1 in order to obtain a second-order in time scheme. Taking ℓ = 0 would lead to space discontinuities preventing the time scheme to be of second-order.

The scheme is finally summarized in the following algorithm. • Initialization of (x 0 k , v 0 k ), ω 0 k and ρ 0 i . At each time step: Prediction step: from t n to t n+1/2 .

• 1) Advance micro part:

advance the characteristics with (4.9), -compute vg with (3.11) and B-spline functions of order ℓ ≥ 1, -advance the equation on the weights with (4.10).

• 2) Projection step: compute (I -• )g n+1/2 using [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF].

• 3) Advance macro part:

compute vg n+1/2 with (3.11) and B-spline functions of order ℓ ≥ 1, -compute the density with (4.11).

Correction step: from t n to t n+1 .

• 4) Advance micro part:

advance the characteristics with (4.12), -compute vg with (3.11) and B-spline functions of order ℓ ≥ 1, -advance the equation on the weights with (4.13).

• 5) Advance macro part with (4.14).

Extension to the Vlasov-Poisson-BGK case

This section is devoted to the extension of our method to kinetic equation making appear an electric field in the velocity direction. We consider the Vlasov-Poisson-BGK system in the diffusion scaling

∂ t f + 1 ε v∂ x f + 1 ε E∂ v f = 1 ε 2 (ρM -f ), (5.1)
∂ x E = ρ -1, (5.2
)

Ω Edx = 0, ∀t ≥ 0, (5.3)
where

x ∈ Ω = [0, L x ] ⊂ R, ρ(t, x) = R f (t, x, v)dv and M (v) = 1 √ 2π exp -v 2 2
is the absolute Maxwellian. Let f 0 (x, v) = f (t = 0, x, v) the initial distribution function and let consider periodic boundary conditions in x:

f (t, 0, v) = f (t, L x , v), ∀ v ∈ V , E (t, 0) = E (t, L x ) ∀t ≥ 0.
We can extend our schemes to this problem by adapting the computations of Subsections 4.1 and 4.2. We do not give all the details of the computations but insist on difficulties coming from the electric field term and write the resulting schemes.

The second-order reformulated micro-macro system corresponding to (1.3) is

∂ t ρ + 1 ε ∂ x vg = 0, (5.4)
∂ t g = 2 ∆t e -∆t/ε 2 -1 e -∆t/ε 2 + 1 g - 2 ε e -∆t/2ε 2 e -∆t/ε 2 + 1 [vM ∂ x ρ + v∂ x g -∂ x vg M -vM Eρ + E∂ v g] . (5.5)
The limit model is here the drift-diffusion equation coupled to Poisson equation (2.4). Equipped with this system (5.4)-(5.5), we construct the following Lagrangian, second-order in time scheme for (5.1)-(5.2)-(5.3). In the prediction step, characteristics are solved through

x n+1/2 k = x n k + ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 v n k , v n+1/2 k = v n k + ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 E n (x n k ), (5.6)
the weights evolve with

w n+1/2 k = 2e -∆t/ε 2 e -∆t/ε 2 + 1 w n k - ∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 [v n k M (v n k )∂ x ρ n (x n k ) -∂ x v n k g n (x n k ) M (v n k ) -v n k M (v n k )E n (x n k )ρ n (x n k )] L x L v N p (5.7)
and the macro equation is advanced with

ρ n+1/2 i = ρ n i - ∆t 2ε ∂ x vg n+1/2 i + ∆t 2 (1 -e -∆t/ε 2 )∂ x ∂ x ρ n+1/2 i + ρ n i 2 -E n i ρ n i .
(5.8)

In the correction step, characteristics are solved through

x n+1 k = x n k + 2∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 v n+1/2 k , v n+1 k = v n k + 2∆t ε e -∆t/2ε 2 e -∆t/ε 2 + 1 E n+1/2 (x n+1/2 k ), (5.9)
the weights evolve with

w n+1 k = e -∆t/ε 2 w n k - ∆t ε e -∆t/2ε 2 v n+1/2 k M (v n+1/2 k )∂ x ρ n+1/2 (x n+1/2 k ) -∂ x v n+1/2 k g n+1/2 (x n+1/2 k ) M (v n+1/2 k ) -v n+1/2 k M (v n+1/2 k )E n+1/2 (x n+1/2 k )ρ n+1/2 (x n+1/2 k ) L x L v N p .
(5.10) and the macro equation is advanced with

ρ n+1 i = ρ n i - ∆t ε ∂ x vg n+1/2 i + ∆t(1 -e -∆t/ε 2 ) 2 ∂ x ∂ x ρ n+1 i + ρ n i 2 -E n+1/2 i ρ n+1/2 i
.

(5.11)

The limit has been directly written in the macroscopic equation and the diffusion term is managed by a Crank-Nicolson method, in the prediction as well as in the correction step. The previous scheme is a semi-discrete and second-order approximation in time of the original equation (2.3). To get a fully discretized scheme, the space derivatives are approximated exactly in the same way as in Subsection 3.2.

We have the following proposition. • for fixed ε > 0, the scheme is a second-order (in time) approximation of the original model (1.3),

• for fixed ∆t > 0, the scheme degenerates into a second-order (in time) scheme of (1.5).

The scheme is finally summarized in the following algorithm.

Algorithm 5.1.

• Initialize (x 0 k , v 0 k ), ω 0 k , and ρ 0 i .

• Compute E 0 i thanks to FFT or finite differences. At each time step: Prediction step: from t n to t n+1/2 .

• 1) Advance micro part:

advance the characteristics with (5.6), -compute vg with (3.11) and B-spline functions of order ℓ ≥ 1, -advance the equation on the weights with (5.7).

• 2) Projection step: compute (I -• )g n+1/2 using [6].

• 3) Advance macro part:

compute vg n+1/2 with (3.11) and B-spline functions of order ℓ ≥ 1, -compute ρ n+1/2 i with (5.8),

compute E n+1/2 i thanks to FFT or finite differences.

Correction step: from t n to t n+1 .

• 4) Advance micro part:

advance the characteristics with (5.9), -compute vg with (3.11) and B-spline functions of order ℓ ≥ 1, -advance the equation on the weights with (5.10).

• 5) Advance macro part:

compute ρ n+1 i with (5.11),

compute E n+1 i thanks to FFT or finite differences.

Remark 5.1. We propose to use an upwind discretization of the derivative

∂ x E n+1/2 i ρ n+1/2 i : ∂ x E n+1/2 i ρ n+1/2 i ≈ E n+1/2,+ i ρ n+1/2 i + E n+1/2,- i ρ n+1/2 i+1 ∆x - E n+1/2,+ i-1 ρ n+1/2 i-1 + E n+1/2,- i-1 ρ n+1/2 i ∆x ,
where the standard notations u + = max(u, 0) and u -= min(u, 0) are used. The same discretization is done for ∂ x (E n i ρ n i ) in the prediction step.

Numerical results

This section is devoted to some numerical experiments comparing the here designed micro-macro model with particles (denoted by MiMa-Part-1 for the first-order scheme and by MiMa-Part-2 for the second-order scheme) to the micro-macro Eulerian model (denoted by MiMa-Grid), the moment guided method (denoted by Moment G.) and (i) the Full PIC method in kinetic regimes (ε of order 1) or (ii) the limit scheme in diffusion regime (ε ≪ 1). MiMa-Part-1 corresponds to Proposition 3.1 and MiMa-Part-2 corresponds to Proposition 4.1. The micro-macro Eulerian scheme is presented in Appendix A. The moment guided method was first presented in [START_REF] Degond | The Moment Guided Monte Carlo Method[END_REF] and is adapted to our context in Appendix B. The Full PIC method [START_REF] Birdsall | Plasma Physics via Computer Simulation[END_REF] consists in applying the particle representation (3.6) to the whole function f (and not only to the perturbation g) and to solve the characteristics and the equations on weights coming from equation (1.1) or (1.3).

In the sequel, we consider three families of test cases: radiative transport equation (RTE) test cases with periodic boundary conditions (see equation (1.1)) in Subsection 6.1, Vlasov-Poisson-BGK test cases (see equations (1.3)-(1.4)) of Landau damping type in Subsection 6.2 and two-stream instability (TSI) test cases in Subsection 6.3.

RTE with periodic boundary conditions.

We consider the RTE test case given by the initial condition

f (t = 0, x, v) = 1 + cos 2π x + 1 2 , x ∈ [0, 1] , v ∈ [-1, 1] , (6.1)
with

M (v) = 1, ∀v ∈ [-1, 1]
, and periodic boundary conditions in x.

We propose here to verify numerically the convergence of MiMa-Part-2 (presented in Subsection 4.3), Lagrangian in space, of second-order in time and with an implicit treatment of the diffusion term (see Algorithm 4.1). In Figure 1, we plot the error in L ∞ norm of the density ρ at time T = 0.1 as a function of ∆t (from 10 -4 to 0.1) for the following parameters: N x = 16, N p = 1000. For ε ≥ 10 -3 , the reference solution is computed with MiMa-Part-2 using the same parameters but with ∆t = 10 -7 . Whereas for ε < 10 -3 , the reference is a numerical solution of the diffusion equation (computed on a space grid, with N x = 16 and ∆t = 10 -7 ). In Figure 2, the error in L ∞ norm is now represented as a function of ε for different values of ∆t: 10 -1 , 10 -2 , 10 -3 and 10 -4 . The reference solution is computed as previously. These plots confirm the fact that MiMa-Part-2 is second-order accurate in time for any fixed ε > 0, and also when ε → 0. However, for intermediate regimes (for instance ε = 0.1 and ε = 10 -3 ), order reduction is observed. This is a classical observation for AP schemes. Note that similar behaviour is obtained with L 2 norm.

In Figure 3, we verify the AP property of the MiMa-Part-2 scheme and plot the density ρ(T = 0.1, x) as a function of x for different values of ε: 1, 0.25, 10 -2 and 10 -6 . We take fixed parameters: N x = 64, ∆t = 10 -3 and N p = 10 4 . We compare the solutions obtained by MiMa-Part-2 to a numerical solution of the diffusion equation (computed on a space grid, with N x = 512 and ∆t = ∆x 2 ) and see that the ε-dependent solutions come closer to the diffusion one when ε decreases.

Moreover, we illustrate in Figure 4 the fact that the cost of our method is very small at the limit. For that, we plot the density ρ(T = 0.1, x) as a function of x for ε = 10 -6 (N x = 64 and ∆t = 10 -2 ) and see that N p = 100 is sufficient to represent in a good way (without noise) the density. The numerical cost is then very close to the one of the asymptotic model.

Landau damping.

In this subsection, we present Landau damping test cases in both regimes (kinetic when ε = O(1) and diffusive when ε → 0). This test is known to be relevant to check the accuracy of the numerical method. In particular, conventional PIC methods have difficulties to capture the long-time behaviour due to the statistical noise. The initial distribution function is given by

f (t = 0, x, v) = 1 √ 2π exp - v 2 2 (1 + α cos (kx)) , x ∈ 0, 2π k , v ∈ R, (6.2)
with the wave number k = 0.5 and α = 0.05. For the micro-macro model (1.3), the initial condition is ρ (t = 0, x) = 1 + α cos (kx) and g (t = 0, x, v) = 0. For the limit drift-diffusion equation (1.5), we have ρ (t = 0, x) = 1 + α cos (kx).

We first verify the order in time of the MiMa-Part-2 scheme detailed in Section 5 and plot in Figure 5 the error in L ∞ norm of the density ρ at time T = 0.1 as a function of ∆t (from 10 -4 to 0.1) for the following parameters: N x = 16, N p = 1000. For ε ≥ 10 -3 , the reference solution is computed with MiMa-Part-2 using the same parameters but with ∆t = 10 -7 . Whereas for ε = 10 -6 the reference is a numerical solution of the drift-diffusion equation (computed on a space grid, with N x = 16 and ∆t = 10 -7 ).

Results are similar to the RTE case: the second-order in time is preserved for big and small values of ε but not for intermediate regimes.

We are now interested in more qualitative tests by considering the time history of the electric energy E (t) =

Lx 0 E (t, x) 2 dx in semi-logarithmic scale for different values of ε. We compare the results obtained by MiMa-Part-2 (detailed in Algorithm 5.1) to other schemes: MiMa-part-1, Moment G. and Full PIC (for ε of order 1) or the scheme for the drift-diffusion model (for small values of ε).

We expect that the number of particles that is necessary to represent in a good way the perturbation g in MiMa-Part methods decreases when ε diminishes and consider N p = 10 5 if ε ≥ 0.5, N p = 10 4 if ε = 0.1 and N p = 10 2 if ε = 10 -4 . For comparison, we take the same N p for moment guided and Full PIC methods.

Results for ε = 10 are given in Figure 6. For the four particle methods, we take ∆t = 0.1, and N x = 128. With the same parameters, results for ε = 1 are presented in Figure 7. For ε = 0.5, we consider ∆t = 0.01 and N x 256 for the four particle methods. Results are given in Figure 8. For these three values of ε, the reference is given by MiMa-Grid with N x = N v = 512 and ∆t = ∆x 2 ≈ 6 × 10 -4 . First, we note that the behaviour of the electric energy is well described during time by micro-macro schemes (MiMa-Part and MiMa-Grid). As observed in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF], the Full PIC method suffers from numerical noise. This is due to the probabilistic character of particle methods (for instance the random initialization of particles). To reduce this noise, we should consider more particles, which would increase the numerical cost. As expected, the moment guided method gives better results than the Full PIC one, but suffers however also from this noise. In MiMa-Part schemes, only the perturbation g is represented by particles (not the whole distribution function f ), that is why for the same N p , the noise is lower, which enables to capture the reference solution for large times. In addition, MiMa-Part-2 is closer to the reference MiMa-Grid than the first-order version MiMa-Part-1.

For smaller values of ε, we compare the four AP schemes (MiMa-Part-1, MiMa-Part-2, MiMa-Grid and Moment G.) to the limit scheme. Results for ε = 0.1 are given in Figure 9. Parameters are the following: ∆t = 10 -3 and N x = 128 for particle methods, ∆t = 0.1∆x 2 ≈ 3.5 × 10 -3 and N x = N v = 64 for MiMa-Grid. We observe that MiMa-Part-2 is the best method since it almost coincides with the reference MiMa-Grid method. Finally, results for ε = 10 -4 are given in Figure 10, where we have ∆t = 10 -2 and N x = 128 for particle methods, ∆t = 0.1∆x 2 ≈ 3.5 × 10 -3 and N x = N v = 64 for MiMa-Grid. The asymptotic regime is well recovered by all these AP methods. As in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF], we remark that few particles are sufficient in the particle-micro-macro schemes to describe in a good way the solution when ε is small. The cost is then reduced at the limit.

In Figures 11 and12, we plot the spatial dependency of the densities (at T = 1) obtained by MiMa-Part-2 and Full PIC methods for different numbers of particles for ∆t = 0.1, N x = 128 and ε = 1. The statistical error observed for MiMa-Part-2 is lower than the one observed for Full PIC when N p = 10 5 is fixed. As expected, when N p increases (N p = 10 6 ), the noise decreases for both methods. Even with N p = 10 4 , MiMa-Part-2 gives rise to satisfactory results, whereas the corresponding Full PIC density would be not acceptable (not plotted on Figure 12).

Two stream instability.

We propose now a study in which the perturbation g is not zero initially and consider the Two-Stream Instability (TSI) test case in both regimes (kinetic and diffusive). The initial distribution function is given by

f (t = 0, x, v) = 1 √ 2π v 2 exp - v 2 2 (1 + α cos (kx)) , x ∈ 0, 2π k , v ∈ R, (6.3)
with the wave number k = 0.5 and α = 0.05. The initial condition for the micro-macro model (1.3) is ρ (t = 0, x) = 1 + α cos (kx) and g

(t = 0, x, v) = 1 √ 2π v 2 -1 exp -v 2 2
(1 + α cos (kx)). For the limit drift-diffusion equation (1.5), we have as in the Landau damping case ρ (t = 0, x) = 1 + α cos (kx).

We first verify the order in time of the MiMa-Part-2 scheme detailed in Section 5 and plot in Figure 13 the error in L ∞ norm of the density ρ at time T = 0.1 as a function of ∆t (from 10 -4 to 0.1) for the following parameters: N x = 16, N p = 1000. For ε ≥ 10 -3 , the reference solution is computed with MiMa-Part-2 using the same parameters but with ∆t = 10 -7 . Whereas for ε = 10 -6 , the reference is a numerical solution of the drift-diffusion equation (computed on a space grid, with N x = 16 and ∆t = 10 -7 ).

As for the RTE and the Landau damping cases, the second-order in time is preserved for big and small values of ε but not for intermediate regimes.

We are now interested in the time evolution of the electric energy E (t) = Lx 0 E (t, x) 2 dx in all regimes.

Results for ε = 10 are given in Figure 14. For the four particle methods, we take N p = 10 6 , ∆t = 0.1, and N x = 128. By taking N p = 10 5 , ∆t = 0.1, and N x = 128 for particle methods, we obtain results for ε = 1 presented in Figure 15. For ε = 0.5, we consider N p = 10 5 , ∆t = 0.01 and N x = 256 for the four particle methods. Results are given in Figure 16. For these three values of ε, the reference is given by MiMa-Grid with N x = N v = 512 and ∆t = ∆x 2 ≈ 6 × 10 -4 . The behaviour of the electric energy is well described during time by micro-macro schemes (MiMa-Part and MiMa-Grid). As previously, the Full PIC method, as well as moment guided method suffer from numerical noise. To illustrate the efficiency of the method, we plot f (T = 5, x, v) obtained by the reference MiMa-Grid, by MiMa-Part-2 and by Full PIC for ε = 10 on Figure 17 and for ε = 0.5 on Figure 18. For MiMa-Grid and MiMa-Part-2, f is reconstructed from g, ρ and M , whereas the approximation of f is directly given by the Full PIC scheme. The numerical parameters are the same as previously (see comments on Figures 14 and16). On Figure 17, we observe that the result obtained by MiMa-Part-2 and Full PIC are in good agreement with MiMa-Grid; however, some numerical noise can be distinguished on f obtained by the Full PIC method. On Figure 18 (ε = 0.5), we can see clearly that the level of the noise is higher for Full PIC, which prevents it from giving good results. On the contrary, MiMa-Part-2 produces good results compared to MiMa-Grid, since the noise only affects the micro part g, which is small in this regime. for MiMa-Grid.

For smaller values of ε, we compare the four AP schemes (MiMa-Part-1, MiMa-Part-2, MiMa-Grid and Moment G.) to the limit scheme. Results for ε = 10 -1 are given in Figure 19. Parameters are the following: ∆t = 10 -3 , N x = 128 and N p = 10 4 for particle methods, ∆t = 0.1∆x 2 ≈ 3.5 × 10 -3 and N x = N v = 64 for MiMa-Grid. As in the Landau damping case, MiMa-Part-2 gives the best result comparing to the reference MiMa-Grid. Finally, results for ε = 10 -4 are given in Figure 20, where we have ∆t = 10 -2 , N x = 128 and N p = 100 for particle methods, ∆t = 0.1∆x 2 ≈ 3.5 × 10 -3 and N x = N v = 64 for MiMa-Grid. The asymptotic regime is well recovered by all these AP methods.

Conclusion

In this paper, we have presented new micro-macro models for the kinetic radiative transport equation (RTE), as well as for the Vlasov-Poisson-BGK system, in the diffusion scaling with periodic boundary conditions. First-order in time and second-order in time models are derived, and their Lagrangian discretizations are detailed. The obtained schemes are proved to degenerate into implicit discretizations of the limit model (the diffusion equation in the RTE case and the drift-diffusion equation in the Vlasov-Poisson-BGK case) when ε → 0. This asymptotic property is shown in the numerical results too. Moreover, thanks to the use of particle methods for the microscopic equation, the numerical cost is reduced when ε diminishes. Finally, compared to a standard PIC method (where f is represented by particles, and not g), the numerical noise is reduced.

In future works, we would like to extend the present approach to high-dimensional Vlasov-Maxwell-BGK case on the one side. On the other side, we would like to combine our approach to Monte Carlo ones to handle diffusion or drift-diffusion limits, in the spirit of what has been proposed for the hydrodynamic limit of Vlasov-BGK in [START_REF] Crouseilles | Asymptotically stable and time diminishing schemes for rarefied gas dynamics[END_REF]. Indeed, this would enable us to adapt automatically the number of particles with respect to ε which is of great importance in applications (recall that in the present work, N p has to be fixed at the beginning). The main idea is to replace the equation of the weights ω k by an equation on the velocities v k to take into account the source part (see [START_REF] Crouseilles | Asymptotically stable and time diminishing schemes for rarefied gas dynamics[END_REF] for more details), as usual in Monte Carlo PIC methods.

A Time discretization for Eulerian schemes

We present the time discretization of (2.1) having in spirit a Eulerian discretization of the phase space. Obviously, the numerical scheme proposed in [START_REF] Bennoune | Uniformly stable numerical schemes for the Boltzmann equation preserving the compressible Navier-Stokes asymptotics[END_REF][START_REF] Crouseilles | An asymptotic preserving scheme based on a micro-macro decomposition for collisional Vlasov equations: diffusion and high-field scaling limits[END_REF][START_REF] Lemou | A new asymptotic preserving scheme based on micro-macro formulation for linear kinetic equations in the diffusion limit[END_REF] works well. Now, (3.5) also provides a numerical scheme that we will exploit in this appendix.

Let us consider staggered grids in the phase-space domain and adopt the following notations: x i = i∆x and x i+1/2 = i∆x + ∆x/2, i ∈ N, define two grids in space and v j = j∆v, j ∈ N, defines a grid in velocity, where ∆x (resp. ∆v) is the step in space (resp. in velocity). Time is also discretized with a time step ∆t and we note t n = n∆t, n ∈ N. The density ρ is discretized on the first space grid: ρ n i approximates ρ(t n , x i ), whereas the perturbation g is discretized on the second one: g n i+1/2,j approximates g(t n , x i+1/2 , v j ).

Let an approximation D of the spatial derivative, the numerical scheme we propose consists in computing g n+1 i+1/2,j with

g n+1 i+1/2,j = e -∆t/ε 2 g n+1 i+1/2,j -ε(1 -e -∆t/ε 2 ) v j ρ n i+1 -ρ n i ∆x +(I -• ) v + j (D - x g n ) i+1/2,j + v - j (D + x g n ) i+1/2,j , (A.1)
where h i+1/2,j = ( j h i+1/2,j ∆v), and then to compute ρ n+1 i with

ρ n+1 i = ρ n i - ∆t ε j v j g n+1 i+1/2,j -g n+1 i-1/2,j ∆x ∆v. (A.2)
Proposition A.1. The scheme given by (A.1)-(A.2) enjoys the AP property, i.e. it satisfies the following properties

• for fixed ε > 0, the scheme is a first-order (in time) approximation of the original model (1.1),

• for fixed ∆t > 0, the scheme degenerates into an explicit first-order (in time) scheme of (1.2).

Proof. We observe easily that when ε goes to zero, (A.1) gives g n+1 i+1/2,j = -εv j ρ n i+1 -ρ n i ∆x + O(ε 2 ), which, injected in the time discretization (A.2) for ρ, gives up to terms of order O(ε 2 )

ρ n+1 i = ρ n i + ∆t   j v 2 j ∆v   ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2 .
Since j v 2 j ∆v is an approximation of -∆t ε j v j h i+1/2,j -h i-1/2,j ∆x ∆v.

(A.3) Since h i+1/2,j = O(ε 2 ) as ε goes to zero after two iterations, the asymptotic preserving property is ensured. Moreover, the diffusion term can now be chosen as implicit, so that the macro equation becomes

ρ n+1 i = ρ n i + ∆t(1 -e -∆t/ε 2 ) j (v 2 j )∆v ρ n+1 i+1 -2ρ n+1 i + ρ n+1 i-1 ∆x 2
-∆t ε j v j h i+1/2,j -h i-1/2,j ∆x ∆v, (A.4) and the scheme is now free from the usual diffusion condition on the time step.

The algorithm finally writes

Algorithm A.1.

• Initialize g 0 i+1/2,j and ρ 0 i . At each time step:

• Advance micro part with (A.1).

• Advance macro part with (A.4).

And we have the following result.

Proposition A.2. The scheme given by (A.1)-(A.4) enjoys the AP property, i.e. it satisfies the following properties

• for fixed ε > 0, the scheme is a first-order (in time) approximation of the original model (1.1),

• for fixed ∆t > 0, the scheme degenerates into an implicit first-order (in time) scheme of (1.2).

We do not present here the extension to the Vlasov-Poisson-BGK case, but it is straightforward.

B Moment guided

In this section, we present the adaptation of the moment guided particle method proposed in [START_REF] Degond | The Moment Guided Monte Carlo Method[END_REF] to our context. For the sake of simplicity, we present it in the RTE case but note that these computations can also be extended to the Vlasov-Poisson-BGK case, without difficulty. The kinetic equation on f has to be reformulated to avoid the singularity linked to the transport term. To do that, we proceed as previously, but from (1.1). Indeed, we rewrite equation (1.1) as

∂ t (e t/ε 2 f ) = e t/ε 2 ε -v∂ x f + 1 ε ρ ,
and we integrate between t n and t n+1 to get f (t n+1 ) = e -∆t/ε 2 f (t n ) -e -t n+1 /ε 2 ε t n+1

t n e t/ε 2 v∂ x f -1 ε ρ dt.

We make the following approximation

f n+1 = e -∆t/ε 2 f n -ε(1 -e -∆t/ε 2 ) v∂ x f n - 1 ε ρ n ,
where f n ≈ f (t n ) and ρ n ≈ ρ(t n ), ∀n.

Making appear the discrete time derivative enables to write

f n+1 -f n ∆t = e -∆t/ε 2 -1 ∆t f n -ε 1 -e -∆t/ε 2 ∆t v∂ x f n - 1 ε ρ n ,
which we approximate by

∂ t f = e -∆t/ε 2 -1 ∆t f -ε 1 -e -∆t/ε 2 ∆t v∂ x f - 1 ε ρ . (B.1)
Following the spirit of the moment guided method (see [START_REF] Degond | The Moment Guided Monte Carlo Method[END_REF]), this equation is coupled with the macro one, that is

∂ t ρ + 1 ε ∂ x vf = 0, ∂ t f + ε 1 -e -∆t/ε 2 ∆t v∂ x f = e -∆t/ε 2 -1 ∆t f + 1 -e -∆t/ε 2 ∆t ρ.
To derive an AP scheme for this latter system satisfied by (ρ, f ), we adapt the strategy presented in [START_REF] Degond | The Moment Guided Monte Carlo Method[END_REF] to our diffusion framework. To do so, we first remark that vf = vg and using the expression of g obtained by (3.5), we get the following approximation for the macro flux (considered implicit in time) A Lagrangian method can be used to approximate the equation on f . As for the micro-macro scheme, we use a splitting procedure

• solve ∂ t f + ε 1 -e -∆t/ε 2 ∆t v∂ x f = 0 • solve ∂ t f = - 1 -e -∆t/ε 2 ∆t f + 1 -e -∆t/ε 2 ∆t ρ.
To do that, the transport part is solved with the (non stiff) characteristics

ẋk (t) = ε 1 -e -∆t/ε 2 ∆t v k (t). (B.3)
The source part is solved using the equation satisfied by the weights ωk (t) = -1 -e -∆t/ε 2 ∆t ω k (t) + 1 -e -∆t/ε 2 ∆t ρ(t, x k (t)). (B.4)

The last step consists in matching the moment of f n+1 obtained by the particle method with ρ n+1 obtained with (B.2). This can be done using the techniques proposed in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF]. Indeed, considering the function g = f -ρ, its weight can be written as

γ k = ω k -β k , with β k = ρ(x k ) L x L v N p .
Then, we apply the discrete version of (I -• ) to the weights γ k as in [START_REF] Crestetto | Kinetic/Fluid micro-macro numerical schemes for Vlasov-Poisson-BGK equations using particles[END_REF]]

ω new k = β k + (I -• )(ω k -β k ).

Proposition 4 . 1 .

 41 The scheme given by (4.9)-(4.10)-(4.11)-(4.12)-(4.13)-(4.14) enjoys the AP property, i.e. it satisfies the following properties

Proposition 5 . 1 .

 51 The scheme given by (5.6)-(5.7)-(5.8)-(5.9)-(5.10)-(5.11) enjoys the AP property, i.e. it satisfies the following properties

Figure 1 :Figure 2 :

 12 Figure 1: Error in L ∞ norm of ρ at time T = 0.1 as a function of ∆t for N x = 16, N p = 1000, ε = 1, 0.5, 0.1, 10 -3 and 10 -6 .

Figure 3 :

 3 Figure 3: AP property. Density ρ(T = 0.1, x) for ε = 1, 0.25, 10 -2 and 10 -6 . N x = 64, ∆t = 10 -3 and N p = 10 4 . Comparison with the diffusion solution.

Figure 4 :

 4 Figure 4: Cost at the limit. Density ρ(T = 0.1, x) for ε = 10 -6 , N x = 64, ∆t = 10 -2 and N p = 10 4 or N p = 100. Comparison with the diffusion solution.

Figure 5 :

 5 Figure 5: Error in L ∞ norm of ρ at time T = 0.1 as a function of ∆t for N x = 16, N p = 1000, ε = 1, 0.5, 0.1, 10 -3 and 10 -6 .

Figure 6 :

 6 Figure 6: Time history of the electric energy, ε = 10. ∆t = 0.1, N x = 128 and N p = 10 5 for the four particle methods.

Figure 7 :

 7 Figure 7: Time history of the electric energy, ε = 1. ∆t = 0.1, N x = 128 and N p = 10 5 for the four particle methods.

Figure 8 :

 8 Figure 8: Time history of the electric energy, ε = 0.5. ∆t = 0.01, N x = 256 and N p = 10 5 for the four particle methods.

Figure 9 :Figure 10 : 6 Figure 11 : 6 Figure 12 :

 910611612 Figure 9: Time history of the electric energy, ε = 0.1. ∆t = 10 -3 , N x = 128 and N p = 10 4 for the four particle methods.

Figure 13 :

 13 Figure 13: Error in L ∞ norm of ρ at time T = 0.1 as a function of ∆t for N x = 16, N p = 1000, ε = 1, 0.5, 0.1, 10 -3 and 10 -6 .

Figure 14 :

 14 Figure 14: Time history of the electric energy, ε = 10. ∆t = 0.1, N x = 128 and N p = 10 6 for the four particle methods.

Figure 15 :

 15 Figure 15: Time history of the electric energy, ε = 1. ∆t = 0.1, N x = 128 and N p = 10 5 for the four particle methods.

Figure 16 :

 16 Figure 16: Time history of the electric energy, ε = 0.5. ∆t = 0.01, N x = 256 and N p = 10 5 for the four particle methods.

Figure 17 :Figure 18 :

 1718 Figure 17: Representation of f (T = 5, x, v), ε = 10. Distribution function f reconstructed from MiMa-Grid on the left, from MiMa-Part-2 on the middle and obtained by Full PIC on the right. ∆t = 0.1, N x = 128 and N p = 10 6 for both particle methods. N x = N v = 512 and ∆t ≈ 6 × 10 -4 for MiMa-Grid.

Figure 19 :

 19 Figure 19: Time history of the electric energy, ε = 0.1. ∆t = 10 -3 , N x = 128 and N p = 10 4 for the four particle methods.

Figure 20 :

 20 Figure 20: Time history of the electric energy, ε = 10 -4 . ∆t = 0.01, N x = 128 and N p = 100 for the four particle methods.

1 - 1

 11 v 2 dv = 1/3, we obtain a consistent discretization of the diffusion equation.Proposition A.1 is of big interest for impliciting the diffusion term ∂ xx ρ. Indeed, let us rewrite (A.1) as followsg n+1 i+1/2,j = -ε(1 -e -∆t/ε 2 )v j ρ n i+1 -ρ n i ∆x + h i+1/2,j ,withh i+1/2,j = e -∆t/ε 2 g n+1 i+1/2,j -ε(1 -e -∆t/ε 2 ) (I -• ) v + j (D - x g n ) i+1/2,j +v - j (D + x g n ) i+1/2,j .Injecting this relation into the macro part, we getρ n+1 i = ρ n i + ∆t(1 -e -∆t/ε 2 ) j (v 2 j )∆v ρ n i+1 -2ρ n i + ρ n i-1 ∆x 2

  1 ε ∂ x vg n+1 = -(1 -e -∆t/ε 2 )∂ xx ρ n + 1 ε e -∆t/ε 2 ∂ x vg n .Then, we get the following scheme for ρρ n+1 = ρ n + ∆t(1 -e -∆t/ε 2 )∂ xx ρ n -∆t ε e -∆t/ε 2 ∂ x vf n . (B.2)

Acknowledgements

N. Crouseilles and M. Lemou are supported by the French ANR project MOONRISE ANR-14-CE23-0007-01 and by the Enabling Research EUROFusion project CfP-WP14-ER-01/IPP-03. A. Crestetto is supported by the French ANR project ACHYLLES ANR-14-CE25-0001 and by the French ANR project MoHyCon ANR-17-CE40-0027-01.


Emile Devillers

Jean Le Besnerais

Thierry Lubin

Michel Hecquet

Jean Philippe Lecointe

A review of subdomain modeling techniques in electrical machines: performances and applications

Keywords: Magnetic field, Electric machines, Analytical model, Harmonic analysis, Performance analysis, Analytical model, Reviews

This paper reviews the recent developments of semi-analytical subdomains modeling techniques to compute the flux density distribution in electrical machines by the exact solving of Maxwell equations.

It is shown that with an appropriate development methodology and numerical implementation, these harmonic models break the traditional compromise between accuracy and computation time that must be done using finite element or other analytical methods. Besides that, subdomains model development techniques have been improved to overcome its topological limitations This fact is demonstrated on three different subdomains models in comparison with finite element methods in terms of accuracy and processing time. The first one is a subdomains model of a surface permanent-magnet synchronous machine, the second one is for an inset permanentmagnet synchronous machine, and the third one is for a squirrelcage induction machine. Thanks to an efficient implementation method, a very low computation time is obtained. The robustness of the subdomains on the geometrical assumptions is also demonstrated.

I. INTRODUCTION

It is often necessary to estimate the rated power, efficiency, magnetic losses and even magnetic vibrations of electrical machines to optimize their design. This estimation relies on an accurate computation of the machine characteristics, such as electromagnetic torque, magnetic losses and air gap Maxwell forces. All these quantities can be computed if the magnetic field inside the machine is fully determined, meaning the space and time distribution for both radial and tangential components.

For this purpose, different methods have been developed and can be grouped in four main categories: analytical, semianalytical, numerical, and hybrid methods which result from the combination of the three first. Numerical methods are very flexible to various geometries, include non-linear and non-homogeneous materials, and enable coupling with other physics. Yet, this high level of complexity induces very timeconsuming simulations which slows down the design process. Simplified analytical models are consequently used for the first design steps as they are very fast and may give more physical insights, while FEM is more interesting for final validations.

Corresponding author: E.Devillers (email: emile.devillers@phd.ec-lille.fr) A semi-analytical method named "Subdomains Model" (SDM) has been recently developed to compute the magnetic field with great accuracy and much faster than numerical methods such as Finite Elements Method (FEM), representing an interesting alternative between analytical and numerical methods.

The first part of this paper presents the principles of the subdomain modeling technique and reviews the works related to SDM development and applications in electrical machines, showing how these works have enlarged the application range of SDMs. Then, a more detailed comparison with FEM is done to demonstrate the performance of the SDMs despite their modeling assumptions, especially in terms of accuracy, computing speed and robustness to geometry. Finally, some future applications of subdomain models are discussed.

II. STATE OF THE ART

A. Principles

SDM is a semi-analytical method that consists in dividing the problem into physical regions named subdomains in which Maxwell governing equations can be solved analytically. The main processing steps to obtain the magnetic field in each subdomain are presented in Fig. 1. For this purpose, the subdomains must fulfill specific conditions on geometry physics. In Fig. 2, the problem is composed of one air gap subdomain, stator slot subdomains and rotor slots subdomains, and limited by rotor and stator iron cores. Then, Maxwell equations are written in each subdomain for scalar potential or vector potential, from which the flux density is derived in step (1) in Fig. 1. The vector potential formulation can always be used whereas the scalar potential formulation may be used only if the magnetic field is irrotational, meaning there is no current density. Maxwell equations are finally solved using the Fourier method consisting in the separation of variables method (2). This analytical resolution gives in each subdomain the potential in Fourier series in function of several unknown constants. The analytical potential solutions may be then reformulated (3) in equivalent expressions for readability and numerical optimization.

By expressing the boundary conditions of potential and magnetic field continuity at each interface between two subdomains, a linear system of analytical independent equations linking all the unknown constants (4) can be obtained. If the potential and the magnetic field of both subdomains are expressed in a different Fourier basis, it is necessary to project one Fourier basis on the other.

Then, assuming a finite number of Fourier harmonics in each subdomain, the linear system of equations is put into matrix form and solved numerically [START_REF] Dubas | Two-dimensional analytical permanentmagnet eddy-current loss calculations in slotless PMSM equipped with surface-inset magnets[END_REF].

M * X = S (1)
The matrix M is often called "topological matrix", and S the "source vector". Solving the linear system gives the numerical value of the unknown constants X, which enables to compute the magnetic field's spatial and temporal distribution in each subdomain [START_REF] Gysen | General Formulation of the Electromagnetic Field Distribution in Machines and Devices Using Fourier Analysis[END_REF]. Because both analytical and numerical resolutions are successively accomplished, SDM may be classified as "semi-analytical" model. A similar methodology is also proposed in [START_REF] Teymoori | Two-Dimensional Analytical Magnetic Field Prediction for Consequent-Pole Permanent Magnet Synchronous Machines[END_REF].

It is important to mention that the formulation of the analytical solution (3) may strongly differ from an author to another. In this review, two main formulations are distinguished : the formulation A in [START_REF] Wu | An Improved Subdomain Model for Predicting Magnetic Field of Surface-Mounted Permanent Magnet Machines Accounting for Tooth-Tips[END_REF], [START_REF] Rahideh | Analytical magnetic field calculation of slotted brushless permanent-magnet machines with surface inset magnets[END_REF] and the formulation B in [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF], [START_REF] Dubas | Two-dimensional analytical permanentmagnet eddy-current loss calculations in slotless PMSM equipped with surface-inset magnets[END_REF]. Formulation B presents the interest of giving dimensionless expressions and a topological matrix M with only 1 on the diagonal and more 0 elsewhere. Specific algorithms are dedicated to optimize the inversion of such matrix.

Hence the choice of formulation is a real matter to design a SDM as it may compromise the numerical resolution during (5) if the topological matrix obtained after (4) is illconditioned. Such numerical problems are frequently pointed out in the SDM literacy, though only few articles such as [START_REF] Gysen | General Formulation of the Electromagnetic Field Distribution in Machines and Devices Using Fourier Analysis[END_REF] analyze the analytical formulation as regards on the numerical performances.

B. Development history

The main difficulties encountered by any analytical models are how to take into account slotting effect, as air gap length variations strongly influence the magnitude and shape of the magnetic field. The first methods based on the formal resolution of Maxwell equations were developed in the 1980's for both slotless Permanent Magnet Synchronous Machine (PMSM) [START_REF] Boules | Two-dimensional field analysis of cylindrical machines with permanent magnet excitation[END_REF] and Induction Machines (IM) [START_REF] Bolte | Three-dimensional analysis of linear motor with solid iron secondary[END_REF], and have been improved by several approaches to better account for slotting effect and radial and tangential air gap flux components .

In 1984, [START_REF] Boules | Two-dimensional field analysis of cylindrical machines with permanent magnet excitation[END_REF] used Carter's coefficients to transform a slotted stator into an equivalent slotless one. In 1993, [START_REF] Zhu | Instantaneous magnetic field distribution in brushless permanent\nmagnet DC motors. III. Effect of stator slotting[END_REF] introduced a relative permeance which modulates the radial air gap flux density previously computed without slotting effect. Another permeance model was developed in 1997 by [START_REF] Brudny | Modelisation de la denture des machines asynchrones : phenomenes de resonances[END_REF]. The relative permeance method was extended by [START_REF] Kim | Magnetic field calculation in permanent magnet motors with rotor eccentricity: Without slotting effect considered[END_REF] in 1998 to take into account both radial and tangential components. In 2003, [START_REF] Wang | Analytical calculation of air-gap magnetic field distribution and instantaneous characteristics of brushless DC motors[END_REF] used conformal transformation and more specifically Schwarz-Christoffel mapping to model slotting effect. This method was also adopted by [START_REF] Žarko | Analytical calculation of magnetic field distribution in the slotted air gap of a surface permanent-magnet motor using complex relative air-gap permeance[END_REF] in 2006, which applied the conformal transformation to the relative permeance model and deduced a complex permeance model, giving better accuracy for both components. One can refer to [START_REF] Gysen | General Formulation of the Electromagnetic Field Distribution in Machines and Devices Using Fourier Analysis[END_REF], [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motor taking into account slotting effect[END_REF], [START_REF] Zhu | An Accurate Subdomain Model forMagnetic Field Computation in Slotted Surface-Mounted Permanent-Magnet Machines[END_REF], [START_REF] Pfister | Slotless permanent-magnet machines: General analytical magnetic field calculation[END_REF] for their exhaustive history in PMSM analytical modeling, and to [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF] for IM modeling.

The first SDM for SPMSM were developed in 2008-2009 by [START_REF] Dubas | Analytical solution of the magnetic field in permanent-magnet motor taking into account slotting effect[END_REF], [START_REF] Liu | Accurate Prediction of Magnetic Field and Magnetic Forces in Permanent Magnet Motors Using an Analytical Solution[END_REF] and [START_REF] Bellara | Two-Dimensional exact analytical solution of armature reaction field in slotted surface mounted PM radial flux synchronous machines[END_REF], although the Fourier projection between subdomains was already used by [START_REF] Ackermann | Analytical modeling of the cogging torque in permanent magnet motors[END_REF] and [START_REF] Bolte | Analytical calculation of the two-dimensional field in the air gap and the slots of electrical machines[END_REF] a few decades ago. Also in 2008, [START_REF] Gysen | Semi-analytical calculation of the armature reaction in slotted tubular permanent magnet actuators[END_REF] used the same method for a linear actuator. The term "Subdomains model" appeared in 2010 in [START_REF] Zhu | An Accurate Subdomain Model forMagnetic Field Computation in Slotted Surface-Mounted Permanent-Magnet Machines[END_REF], [START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF] and was adopted by then in several contemporary major publications. This method is also referred as "exact analytical model" [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF], "semianalytical harmonic model" [START_REF] Sprangers | Electric circuit coupling of a slotted semi-analytical model for induction motors based on harmonic modeling[END_REF] or "Fourier-based Model" [START_REF] Hannon | 2-D Analytical Subdomain Model of a Slotted PMSM With Shielding Cylinder[END_REF]. Compared to the previous analytical models, SDMs provide both components of the magnetic field by exactly taking into account slotting effect and the influence between slots. In 2010, [START_REF] Lubin | Exact Analytical Method for Magnetic Field Computation in the Air Gap of Cylindrical Electrical Machines Considering Slotting Effects[END_REF] developed an elementary model to give a better understanding of slotting effect in SDMs.

C. Available topologies 1) Introduction:

Due to the PMSM popularity in the past decades, most of SDMs deal with them, at the detriment of IM. This can also be explained because of more complex physics. Besides, SDMs for other types of machines have been developed. A complete review on SDMs done until 2014 is presented in [START_REF] Pfister | Slotted Permanent-Magnet Machines: General Analytical Model of Magnetic Fields, Torque, Eddy Currents and Permanent Magnet Power Losses including the Diffusion Effect[END_REF]. Some examples of existing topologies are illustrated in Fig. 3.

2) Geometry aspects: The geometry is usually in twodimensional (2-D), but have already been extended to 3-D such as in [START_REF] De La Barrière | 3-D formal resolution of maxwell equations for the computation of the no-load flux in an axial flux permanent-magnet synchronous machine[END_REF], [START_REF] Meessen | General formulation of fringing fields in 3-D cylindrical structures using fourier analysis[END_REF]. In 2-D, axial and radial 3-D end-effects are neglected.

The problem is either expressed in polar or Cartesian coordinates. In case of polar coordinates, every subdomain geometry is approximated by a polar geometry. For example, the rectangular teeth are supposed to have radial edges with orthoradial tooth tips, as it is illustrated in [START_REF] Wu | An Improved Subdomain Model for Predicting Magnetic Field of Surface-Mounted Permanent Magnet Machines Accounting for Tooth-Tips[END_REF]. In case of Cartesian coordinates, the air gap is unrolled by considering an infinite radius of curvature, giving an equivalent rectangular topology such as in [START_REF] Tiegna | Analytical Modeling of the Open-Circuit Magnetic Field in Axial Flux Permanent-Magnet Machines With Semi-Closed Slots[END_REF]. The different topology approximations according to the chosen coordinate systems are gathered in the aforementioned methodology [START_REF] Gysen | General Formulation of the Electromagnetic Field Distribution in Machines and Devices Using Fourier Analysis[END_REF].

Besides, SDMs can be applied to internal or external rotor [START_REF] Zhu | An Accurate Subdomain Model forMagnetic Field Computation in Slotted Surface-Mounted Permanent-Magnet Machines[END_REF], [START_REF] Jahromi | 2-D Analytical Model for External Rotor Brushless PM Machines[END_REF]. Moreover, semi-closed slots may be used for a more realistic model [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF].

3) Physics aspects: As said previously, the physics is also approximated. The iron is considered to have infinite relative permeability, resulting in homogeneous boundary conditions at the interfaces between subdomains and the iron. The saturation is consequently neglected. For PM machines, magnets have an isotropic and homogeneous relative permeability and a linear B(H) curve. For induction machines (IM), the rotor bars are assumed to have a homogeneous electrical conductivity.

SDM can model both magnet and current sources. It accounts for any magnetization shapes such as radial, parallel or Hallbach magnetization. Concerning the armature reaction fields, the windings are usually designed by a connection matrix which enables to use (non)-overlapping single/double layer windings. In fact, the magnetic sources are expanded into Fourier series and injected in Maxwell equations at step [START_REF] Wu | An Improved Subdomain Model for Predicting Magnetic Field of Surface-Mounted Permanent Magnet Machines Accounting for Tooth-Tips[END_REF]. It results in a linear superposition of stator and rotor fields, hence the possibility to solve everything at once or separately.

4) Synchronous Machines (SM):

Several models exist for each topology of PMSM, depending on the chosen modeling level. For SPMSM with armature reaction field and semiclosed slots, one can refer to [START_REF] Wu | An Improved Subdomain Model for Predicting Magnetic Field of Surface-Mounted Permanent Magnet Machines Accounting for Tooth-Tips[END_REF], [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF], [START_REF] Jahromi | 2-D Analytical Model for External Rotor Brushless PM Machines[END_REF]. Inset PMSM (IPMSM) models with armature reaction field and semiclosed slots can be found in [START_REF] Teymoori | Two-Dimensional Analytical Magnetic Field Prediction for Consequent-Pole Permanent Magnet Synchronous Machines[END_REF], [START_REF] Rahideh | Analytical magnetic field calculation of slotted brushless permanent-magnet machines with surface inset magnets[END_REF], [START_REF] Lubin | Two-Dimensional Analytical Calculation of Magnetic Field and Electromagnetic Torque for Surface-Inset Permanent-Magnet Motors[END_REF].

Besides the above topologies, more singular SM have been modeled by the subdomain technique. For instance, SDMs exist for flux switching SM [START_REF] Boughrara | General Subdomain Model for Predicting Magnetic Field in Internal and External Rotor Multiphase Flux-Switching Machines Topologies[END_REF], double excitation SM [START_REF] Bali | Analytical Modeling of Open Circuit Magnetic Field in Wound Field and Series Double Excitation Synchronous Machines[END_REF], axial flux SM [START_REF] Tiegna | Analytical Modeling of the Open-Circuit Magnetic Field in Axial Flux Permanent-Magnet Machines With Semi-Closed Slots[END_REF], PMSM with noches [START_REF] Oner | Analytical On-Load Sub-domain Field Model of Permanent Magnet Vernier Machines[END_REF], pseudo direct drives SM [START_REF] Penzkofer | Analytical Modeling and Optimization of Pseudo-Direct Drive Permanent Magnet Machines for Large Wind Turbines[END_REF] and Switch Reluctance Motor (SRM) [START_REF] Sprangers | Semi-Analytical Framework for Synchronous Reluctance Motor Analysis including Finite Soft-Magnetic Material Permeability[END_REF]. The SRM model illustrates the difficulties to transform any geometry into a polar one.

5) Induction Machines (IM): IM SDMs have both common points and differences with SM ones. Assuming an internal rotor topology, stator slots subdomains and air gap subdomain remain the same as for PMSM. Though, the level complexity is increased because of the the induced current modeling in the rotor bars, and of the existence of two asynchronous frequencies as well as the space harmonics.

The former analytical models of IM were designed for laminated solid rotor, such as in [START_REF] Bolte | Three-dimensional analysis of linear motor with solid iron secondary[END_REF], [START_REF] Gerling | Three-dimensional analytical calculation of induction machines with multilayer rotor structure in cylindrical coordinates[END_REF], and extended by the SDM in [START_REF] Boughrara | 2-D Analytical Prediction of Eddy Currents, Circuit Model Parameters, and Steady-State Performances in Solid Rotor Induction Motors[END_REF].

The first model of Squirrel Cage IM (SCIM) in [START_REF] Lubin | Analytic Calculation of Eddy Currents in the Slots of Electrical Machines: Application to Cage Rotor Induction Motors[END_REF] accounts for rotor bars with induction and a current sheet at stator inner bore. A complete SCIM SDM was developed in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF], [START_REF] Sprangers | Electric circuit coupling of a slotted semi-analytical model for induction motors based on harmonic modeling[END_REF], using Electrical Equivalent Circuit (EEC) to drive the feeding current as a function of the slip value. 6) Linear or tubular PM machines: SDMS of linear actuators can be found in [START_REF] Gysen | Semi-analytical calculation of the armature reaction in slotted tubular permanent magnet actuators[END_REF], [START_REF] Shin | Analysis on the Pitching Moment in Permanent Magnet Linear Synchronous Motor for Linear Motion Stage Systems[END_REF]. These topologies are well suited to SDMs because the geometry can be directly implemented in Cartesian coordinates without simplifications.

7) Missing machines:

To the author's knowledge, no paper specifically deals with SDM of Direct-Current Machines (DCM) but the models used for synchronous machines can be applied with DC current sources. There is either no Buried PMSM (BPMSM) SDM, because of the saturation in the bridges between magnets and the difficulty to approximate interior magnets with polar geometry.

D. Accounting for time parameter

In SDMs, time is differently taken into account according to the machine type. Its influence on the topological matrix M and the source vector S must be carefully established. In fact, sources variations impact on S whereas reluctance variations modify M.

Particularly, for SPMSM with internal rotor and magnet permeability equal to 1, there is no reluctance variation with rotor rotation. Only the magnetization distribution changes with the rotation and this results in a source matrix whose columns are the source vector S for each time step. Consequently, a single linear system resolution can give the constants X for every time step.

For IPMSM, both reluctance and magnetization distribution change over time. In this case, it is the same principle as for FEM : M and S have to be evaluated at each time step and one linear system resolution may only give one time step.

It is a different method for IM. Every stator quantity, respectively rotor quantity, is first expressed as a phasor of pulsation ω, respectively slip * ω. Magnetic potentials are then solved in terms of complex amplitude and finally modulated with their respective pulsation. It implies that one system resolution gives every time step for one chosen slip value.

E. Current applications

Though it may have been restricted to torque and Back-ElectroMotive Force (BEMF) at the beginning, SDMs have been derived for many other applications.

1) Magnetic Forces and Flux linkage computation: Each local and global electromagnetic quantity may be derived from the magnetic potential using well-known physics laws. Electromagnetic forces as well as electromagnetic torque including torque ripple -such as cogging torque -may be directly computed from the magnetic field using the Maxwell stress tensor with a great accuracy. The flux linkage is computed by applying the Stokes theorem to the magnetic potential of each stator slot. Flux linkage knowledge enables to deduce the BEMF generated by the rotating field.

Thanks to the Fourier series formulation, these previous temporal integrations are converted in faster and more accurate summations on the magnetic field's harmonic components. The limit of any 2-D magnetic models is the fact that 3-D axial end-effects have been neglected so the torque and the BEMF may be overestimated.

2) Fault simulation: It is possible to simulate unbalanced magnetic sources distribution by injecting the corresponding Fourier series at the step (2), such as in [START_REF] Wu | An Analytical Model of Unbalanced Magnetic Force in Fractional-Slot Surface-Mounted Permanent Magnet Machines[END_REF], and deduce the Unbalanced Magnetic Force (UMF) -or Pull (UMP). UMF is also caused by rotor eccentricity. The effect of rotor eccentricity is modeled in [START_REF] Fu | Subdomain model for predicting magnetic field in slotted surface mounted permanent-magnet machines with rotor eccentricity[END_REF] by adding a first-order perturbation component to magnetic potentials, whereas [START_REF] Li | Analytical Synthesis of Air-gap Field Distribution in Permanent Magnet Machines with Rotor Eccentricity by Superposition Method[END_REF] introduces it with a superposition method. For SCIM, it is possible to simulate defective bars as in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF] by decreasing their conductivity.

3) Losses: Eddy-current losses may be computed using Helmholtz equation in the PM subdomains [START_REF] Dubas | Two-dimensional analytical permanentmagnet eddy-current loss calculations in slotless PMSM equipped with surface-inset magnets[END_REF], [START_REF] Wu | Simplified analytical model and investigation of open-circuit AC winding loss of permanent-magnet machines[END_REF] and in the windings [START_REF] Amara | Analytical Prediction of Eddy-Current Loss in Armature Windings of Permanent Magnet Brushless AC Machines[END_REF], [START_REF] Arumugam | Estimation of eddy current loss in semi-closed slot vertical conductor permanent magnet synchronous machines considering eddy current reaction effect[END_REF]. In [START_REF] Pfister | Slotted Permanent-Magnet Machines: General Analytical Model of Magnetic Fields, Torque, Eddy Currents and Permanent Magnet Power Losses including the Diffusion Effect[END_REF], a shielding cylinder is added at the surface of the PM to reduce the eddy-current losses in magnets.

4) EEC parameters:

The potential solution is used to compute leakage flux at stator slots [START_REF] Prieto | Analytical Calculation of the Slot Leakage Inductance in Fractional-Slot Concentrated-Winding Machines[END_REF] and due to endwindings [START_REF] Liu | Analytical Solution of the Magnetic Field and EMF Calculation in Ironless BLDC Motor[END_REF]. A method to compute primary and secondary impedances for SCIM is developed in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF], [START_REF] Sprangers | Electric circuit coupling of a slotted semi-analytical model for induction motors based on harmonic modeling[END_REF]. It is also shown that the estimated EEC parameters enable to compute and check previous quantities such as torque, back-EMF, power losses, etc. with another approach.

III. ADVANTAGES AND DRAWBACKS IN COMPARISON WITH FEM A. Introduction

In most of SDMs papers, the model is validated by comparing its accuracy with a parallel FEM analysis using the same approximations. For the same modeling level, SDM is naturally as accurate as FEM, since the former is an exact resolution of Maxwell equations. It is actually more exhaustive than FEM because the solution is continuously defined in each subdomain and not only at the mesh points. Tab. I shows the qualitative performances criteria in comparison with FEM. B. Comparison on the model limitations 1) Sensitivity to mesh in the FEM: Significant problems due to the meshing quality exist, as regards on computing derived magnetic quantities. This has been studied a lot for the evaluation of vibrations due to magnetic forces and cogging torque. For example, Fig 2 .38 in [START_REF]Etudes des vibrations d'origine électromagnétique d'une machine électrique : conception optimisée et variabilité du comportement vibratoire[END_REF] states a vibration variation up to 4 dB below 10 kHz between different meshing methods in Flux3D [START_REF] Flux3d | Electromagnetic and thermal finite element analysis[END_REF]. These problems may be solved by refining the mesh in the air gap, but it significantly slows down the computation.

2) Robustness to geometry in SDMs: One drawback of SDM is the geometry simplification. However by defining an equivalent polar geometry as in Fig. 4, the air gap flux computation is still accurate. In Fig. 5, a comparison has been done using MANATEE simulation environment [START_REF] Manatee | Magnetic acoustic noise analysis tool for electrical engineering[END_REF] on the SPMSM 6s/4p presented in [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF]. 3) Saturation in SDMs: Saturation is still a limitation that prevents from modeling several machine topologies. The first challenge is to compute the magnetic field assuming a finite iron permeability. This is done for example in [START_REF] Boughrara | 2-D Analytical Prediction of Eddy Currents, Circuit Model Parameters, and Steady-State Performances in Solid Rotor Induction Motors[END_REF] where the magnetic field is computed in the rotor iron core. In [START_REF] Sprangers | Magnetic Saturation in Semi-Analytical Harmonic Modeling for Electric Machine Analysis[END_REF], a method is proposed to account for linear soft-magnetic iron with finite permeability, but there exists no SDMs which account for non-linear materials to the authors' knowledge.

However, it is possible to couple SDM with Magnetic Equivalent Circuit (MEC) [START_REF] Gholizad | Dynamic analysis of highly saturated switched reluctance motors using coupled magnetic equivalent circuit and the analytical solution[END_REF], [START_REF] Ilhan | Analytical hybrid model for flux switching permanent magnet machines[END_REF], [START_REF] Laoubi | Hybrid analytical modeling of a flux switching permanent magnet machines[END_REF] or FEM [START_REF] Liu | A combined numerical and analytical approach for magnetic field analysis of permanent magnet machines[END_REF] in an iterative way to account for saturation.

C. Comparison on computing performance of SDMs 1) SDM numerical optimization: An efficient implementation is necessary to take advantage of SDM performances. To illustrate this, the SDM recently presented in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF], a SCIM 36s/4p at no-load condition has been developed and optimized. However, the analytical resolution at (3) is done with reformulation B to improve numerical performances, and the numerical implementation uses only vectorized operations.

Besides, the choice of the harmonics number in each subdomain at step (4) relies on a compromise between accuracy and computation time [START_REF] Tiegna | Overview of analytical models of permanent magnet electrical machines for analysis and design purposes[END_REF]. Fig. 6 shows the variation of the Mean Squared Error (MSE) between the magnetic field and the reference obtained with 1000 air gap harmonics, the computation time and the memory used by M in function of the air gap harmonics number. Harmonics numbers in the other subdomains are the same as in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF].

Fig. 6 illustrates the fast convergence of the Fourier series, as 150 harmonics gives a MSE of 0.2% compared with 1000 harmonics, within 0.34 seconds. In [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF], the SDM computation time is 3.24 minutes and very close to the FEM computation time. Furthermore, a few SDMs for SPMSM accounts for periodicity and one slot per pole approximation [15] [44] to reduce computation time.

The air gap harmonics number can be reduced by a previous analytical study of the magnetic field's harmonic content. For instance, in a SPMSM with Z s stator teeth and 2p poles, the fundamental at no load is at 2p and harmonics are linear combinations of 2p and Z s . Hence any other harmonics ranks can be ignored to improve performances.

2) Comparison with FEM: An efficient implementation of SDMs can lead to much lower computation time than FEM, as presented in [START_REF] Gilson | 2-D Analytical Subdomain Model for High-Speed Permanent-Magnet Machines[END_REF]. Tab. II shows a quantitative comparison on performances between FEM using Femm [START_REF] Femm | Finite element method magnetics[END_REF] and three particular SDMs : the SPMSM 6s/4p without armature load presented in [START_REF] Lubin | 2-D exact analytical model for surface-mounted permanent-magnet motors with semi-closed slots[END_REF], the IPMSM 15s/4p without armature load in [START_REF] Lubin | Two-Dimensional Analytical Calculation of Magnetic Field and Electromagnetic Torque for Surface-Inset Permanent-Magnet Motors[END_REF], and the previous SCIM 36s/4p in [START_REF] Boughrara | Analytical Analysis of Cage Rotor Induction Motors in Healthy, Defective and Broken Bars Conditions[END_REF].

For the SDM, the computation time includes building the linear system, solving it, and reconstructing the air gap magnetic field in the air gap subdomain. For the FEM, it includes building the topology, meshing it, solving the problem in the entire topology and extracting the air gap field. Both models are linear and without any symmetry or periodicity simplifications. When computing temporal quantities such as cogging torque, FEM computation time values in tab. II are multiplied by the number of time steps. As said previously, it is not the case for an SPMSM model whose computation time is very less sensitive to the number of time steps. For instance, considering the same topology of SPMSM 6s/4p, the computation time for the SDM is 0.139 seconds whereas it takes around 347 seconds for the FEM because of the 50 time steps, meaning 50 meshings and resolutions.

SDMs are implemented in the Matlab scientific environment -version R2014b. The computer is equipped with a CPU Intel Xeon E5 1620 v2 @ 3.70GHz and 24GB RAM DDR3 @ 797 MHz. All results in this paper have been repeated 10 times and then averaged for robustness purpose.

IV. CONCLUSION

The paper presents a review of the semi-analytical subdomains models applied to electrical machines. The methodology is first explained step by step, including both analytical and numerical aspects such as model assumptions and the choice of formulation. Then, a state of the art presents the main existing models and their numerous applications. Advantages and drawbacks of the method are finally illustrated in comparison with FEM for three machines topologies. The results validate the model assumptions and show great computation time performances especially when accounting for rotor rotation.

Furthermore, an efficient implementation of subdomains models strongly reduces the computation time and allows to increase the electromagnetic model complexity. This can be done to include:

• the effect of magnet shaping using harmonic superposition • the effect of winding space harmonics in induction machines using field superposition • 3-D effects such as fringing flux and skewing • strong coupling with electrical circuit (calculation of equivalent circuit parameters iteratively) • integration in complex multiphysics models, especially for vibroacoustic analysis [START_REF] Besnerais | Fast prediction of acoustic noise and vibrations due to magnetic forces in electrical machines[END_REF].

Fig. 1 .

 1 Fig. 1. Methodology graph of Subdomains Models in 7 steps.

Fig. 2 .

 2 Fig. 2. Division in polar subdomains for an induction machine.

Fig. 3 .

 3 Fig. 3. Three idealized topologies of existing SDMs. a) SPMSM ; b) IPMSM ; c) SCIM. The pictures are provided by MANATEE software.

Fig. 4 .

 4 Fig. 4. Polar approximation of constant tooth width and curved magnet.

Fig. 5 .

 5 Fig. 5. Magnetic field comparison between idealized and real geometries.

Fig. 6 .

 6 Fig. 6. MSE, computation time and memory used in function of the air gap harmonics number.

TABLE II COMPARISON

 II ON THE ELAPSED SIMULATION TIME FOR ONE TIME STEP.

		SPMSM IPMSM SCIM
	SDM
	Computation time [s]	0.109	0.199	0.7337
	Air gap harmonics	150	250	300
	Number of unknowns 1572	1909	2298
	FEM
	Computation time [s]	6.152	10.11	15.27
	Number of nodes	45888	64605	73003
	Number of elements	91054	128488	145284