Keywords: 2000 Mathematics Subject Classification. 35J20, 35Q31, 35Q60, 65M06, 65M08, 65M12, 65N20, 76N17, 76W05, 76X05 Fusion plasmas, Euler equations, Lorentz force, Large magnetic field, quasi-neutrality, Drift-fluid limit, Asymptotic-Preserving schemes, strongly anisotropic problems, Micro-macro decomposition

1.

Introduction. This paper is devoted to the construction of a numerical scheme for the simulation of a two-fluid Euler-Lorentz model: such a model represents the evolution of a mixture of an ion gas and an electron gas which are submitted to the Lorentz force. More precisely, we focus in this paper on a situation involving a strong Lorentz force and a low Mach number regime for both ion and electron fluids, i.e. we assume that pressure and Lorentz forces are of the same order as τ -1 where τ > 0 is the square of the ion Mach number and also represents the ratio between the ion gyro-period and the characteristic time scale of the experiment. When τ converges to 0, we reach an asymptotic regime which is referred to as the drift-fluid regime or gyro-fluid regime: in this limit regime, the pressure force for the ions and for the electrons is balanced by the Lorentz force. In return, the momentum equations within the two-fluid Euler system degenerate into a pair of equations in which the parallel components of the ion and electron velocities can be viewed as Lagrange multipliers of the zero total force equations in the direction of the magnetic field (see [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF][START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF]).

Such a model describes plasma physics experiments involving strong external magnetic fields, such as Magnetic Confinement Fusion (MCF) experiments in tokamak reactors. In such a case, the rescaled gyro-period of the confined particles, which is denoted by τ , can be close to 0. The assumption that τ is very small leads to a singularly perturbed Euler-Lorentz model. The limit τ → 0 is referred to as the gyro-fluid limit (see [START_REF] Beer | Toroidal gyrofluid equations for simulations of tokamak turbulence[END_REF][START_REF] Dorland | Gyrofluid turbulence models with kinetic effects[END_REF][START_REF] Hammett | Developments in the gyrofluid approach to tokamak turbulence simulations[END_REF]). It is also possible to consider a gyro-kinetic approach when a kinetic model is considered from the onset, instead of fluid equations (see [START_REF] Brizard | Nonlinear gyrokinetic tokamak physics[END_REF][START_REF] Brizard | Foundations of nonlinear gyrokinetic theory[END_REF][START_REF] Hitchon | Flux coordinates and magnetic field structure. A guide toa fundamental tool of plasma theory[END_REF][START_REF] Dubin | Nonlinear gyrokinetic equations[END_REF][START_REF] Frénod | Two-dimensional Finite Larmor Radius approximation in canonical gyrokinetic coordinates[END_REF][START_REF] Grandgirard | A drift-kinetic semilagrangian 4D code for ion turbulence simulation[END_REF][START_REF] Grandgirard | Computing ITG turbulence with a full-f semi-lagrangian code[END_REF][START_REF] Heikkinen | Full-f gyrokinetic method for particle simulation of tokamak transport[END_REF][START_REF] Lee | Gyrokinetic approach in particle simulation[END_REF][START_REF] Lee | Gyrokinetic particle simulation model[END_REF][START_REF] Littlejohn | A guiding center Hamiltonian : A new approach[END_REF]). For generalities on asymptotic regimes for fusion plasmas physics, we refer to [START_REF] Miyamoto | Controlled fusion and plasma physics[END_REF].

In many experimental cases, the value of τ is not uniform and can vary a lot between subdomains of the tokamak. Additionally, it may depend on the time variable: in most of MCF experiments, τ is very small in the plasma core whereas it can be of order 1 far from the plasma core. From a numerical point of view, the usual approach for simulating both cases together consists in a domain decomposition according to the local value of τ . More precisely, we choose to simulate the initial τ -dependent model in the regions where τ = O(1), and we choose the limit model in the regions where τ ≪ 1. Such an approach involves different numerical methods for solving either the Euler-Lorentz model or its drift-fluid limit according to the value of τ . Generally, the coupling of these methods is not straightforward and presents several drawbacks such as the treatment of the interface position (or cross-talk region): indeed, it can depend on the time variable and, in most cases, costly algorithms are required to simulate the motion of the interface and to couple it with the space mesh.

We choose a different method based on the resolution of the τ -dependent Euler-Lorentz model and on the design of a scheme which is able to handle both the cases τ = O(1) and τ ≪ 1. Then this so-called Asymptotic-Preserving (AP) method provides consistent approximation of the Euler-Lorentz model when τ = O(1) and of its limit regime when τ → 0, and does not require a τ -dependent stability condition. As a consequence, such a method can be used on the whole simulation domain for both the τ = O(1) and τ ≪ 1 regimes. AP schemes have been introduced by S. Jin [START_REF] Jin | Efficient Asymptotic-Preserving (AP) schemes for some multiscale kinetic equations[END_REF] and were applied on tranport models and their diffusive limits. Other applications can be found in plasma physics (see [START_REF] Belaouar | emphAn asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF][START_REF] Crispel | Quasi-neutral fluid models for current-carrying plasmas[END_REF][START_REF] Crispel | An asymptotic preserving scheme for the two-fluid Euler-Poisson model in the quasi-neutral limit[END_REF][START_REF] Crispel | A plasma expansion model based on the full Euler-Poisson system[END_REF][START_REF] Degond | Asymptotic-Preserving Particle-In-Cell method for the Vlasov-Poisson system near quasi-neutrality[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for strongly anisotropic elliptic problems[END_REF][START_REF] Degond | Numerical approximation of the Euler-Poisson-Boltzmann model in the quasi-neutral limit[END_REF][START_REF] Degond | Analysis of an Asymptotic-Preserving scheme for the Euler-Poisson system in the quasi-neutral limit[END_REF] for quasi-neutrality regimes and [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF][START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF] for strong magnetic fields regimes), low Mach number fluid dynamics (see [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF][START_REF] Mclarren | The effects of slope limiting on Asymptotic-Preserving schemes numerical methods for hyperbolic conservation laws[END_REF]), or other types of transport problems (see [START_REF] Buet | emphDiffusion limit of the Lorentz model: Asymptotic-Preserving schemes[END_REF][START_REF] Buet | Asymptotic-Preserving and positive schemes for radiation hydrodynamics[END_REF][START_REF] Carrillo | Simulation of fluid and particles flows: Asymptotic-Preserving schemes for bubbling and flowing regimes[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] Filbet | A class of Asymptotic-Preserving schemes for kinetic equations and related problems with stiff sources[END_REF][START_REF] Filbet | An Asymptotic-Preserving Scheme for the ES-BGK model of the Boltzmann equation[END_REF][START_REF] Klar | An Asymptotic-Preserving numerical scheme for kinetic equations in the low Mach number limit[END_REF][START_REF] Lemou | A new Asymptotic-Preserving scheme based on micro-macro decomposition for linear kinetic equations in the diffusion limit[END_REF]) or diffusion problems (see [START_REF] Degond | Duality-based asymptoticpreserving method for highly anisotropic diffusion equations[END_REF]).

The present paper has two main goals: the first one is to propose a new equivalent formulation of the two-fluid Euler-Lorentz model when τ > 0. In this new formulation, the parallel velocity equation for τ = 0 explicitly appears as the limit of the parallel velocity equation for τ > 0 by contrast to the original formulation. This reformulation is the building block for the Asymptotic-Preserving scheme for the two-fluid Euler-Lorentz model. For the scheme being simple, we choose an isothermal pressure law for both ion and electron fluids. We also assume that the fluids are coupled through a quasi-neutrality constraint which allows us to compute the self-consistent electric field within the Lorentz term, the magnetic field being external and given. This work is the last development of a program started in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF][START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF]: in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF] and [START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF], the one-fluid isentropic Euler-Lorentz model with given electric and magnetic fields has been investigated (in [START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF], under a uniform magnetic field and in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF], with any magnetic field and arbitrary coordinate system). Here, the specificity of this work is to consider a two-fluid Euler-Lorentz model with self-consistent electric field, computed through the quasi-neutrality hypothesis. This leads to a system of two coupled anisotropic diffusion equations, which brings some specific difficulties. In particular, it involves some singularity which will be treated through a regularization procedure. This work also bears relations with [START_REF] Degond | Duality-based asymptoticpreserving method for highly anisotropic diffusion equations[END_REF][START_REF] Degond | An Asymptotic-Preserving scheme for strongly anisotropic elliptic problems[END_REF] which are concerned with more general anisotropic diffusion equations (but not in the context of the Euler-Lorentz model). Again, [START_REF] Degond | Duality-based asymptoticpreserving method for highly anisotropic diffusion equations[END_REF] deals with a uniform anisotropy direction and [START_REF] Degond | An Asymptotic-Preserving scheme for strongly anisotropic elliptic problems[END_REF] with an arbitrary anisotropy direction and arbitrary coordinate systems compared to the direction of the anisotropy.

The present paper is organized as follows. In section 2, we present the isothermal two-fluid Euler-Lorentz model and the drift-fluid limit regime. The section 3 is devoted to the reformulation of the τ -dependent Euler-Lorentz model leading to a new equivalent formulation of these equations when τ > 0 which is equivalent to the drift-fluid limit when τ = 0. In section 4, we present a time semi-discretization of the Euler-Lorentz model which is also consistent with the new formulation of the model. Since this time semi-discrete scheme involves an ill-posed diffusion problem for the electric potential, we choose to recover the well-posedness of this problem by introducing a regularization of the quasi-neutrality constraint. This is the subject of the second part of section 4. In section 5, we present a fully-discrete finite volume scheme based on the AP scheme for the Euler-Lorentz model coupled with the perturbed quasi-neutrality constraint. Finally, in section 6, we present some numerical results which have been obtained with this scheme.

2. The isothermal two-fluid Euler-Lorentz model.

2.1.

Scaling. In this paragraph, we present the scaling of the Euler-Lorentz equations which leads to the following dimensionless model:

             ∂ t n τ + ∇ x • q τ α = 0 , (2.1a) ǫ α τ ∂ t q τ α + ∇ x • q τ α ⊗ q τ α n τ +T α ∇ x n τ = q α -n τ ∇ x φ τ + q τ α × B , (2.1b) α ∈ {i, e} , (2.1c)
where τ is the ratio between the ion gyro-period and the characteristic time scale but also the square value of the ion Mach number, T i = 1 and T e = cst. are the dimensionless ion and electron temperatures, and ǫ α and q α are defined by

ǫ α = 1 , if α = i, ǫ , if α = e, q α = 1 , if α = i, -1 , if α = e, (2.2)
where ǫ is the ratio between the unit electron mass and the unit ion mass. Finally, n τ , q τ i , q τ e , φ τ and B correspond to the dimensionless ion and electron density, the dimensionless ion momentum, the dimensionless electron momentum, the dimensionless electric potential and the external magnetic field and are functions of the position x ∈ R 3 and of the time t ≥ 0.

Our starting point is the two-fluid Euler-Lorentz model describing a mixture of an ion gas and an electron gas. This model writes

               ∂ t n α + ∇ x • q α = 0 , ∂ t q α + ∇ x • q α ⊗ q α n α + 1 m α ∇ x p α = q α e m α (n α E + q α × B) , ∂ t e α + ∇ x • e α + p α n α q α = q α e E • q α , α ∈ {i, e} , (2.3)
and the ion-electron coupling is ensured by the following quasi-neutrality constraint

n i = n e = n . (2.4)
In this two-fluid model, n i , q i , e i and p i (resp. n e , q e , e e and p e ) are respectively the density, the momentum, the total energy per mass unit and the pressure for the ion (resp. electron) gas. The physical constants m i , m e , and e stand for the unit ion mass, the unit electron mass and the absolute value of unit electron charge. The electric and magnetic fields are denoted with E and B and we assume that B is given whereas E is generated by the ions and the electrons through the constraint (2.4).

In the present context, we assume that both ion and electron gases are isothermal, i.e. we assume that p i and p e are given by

p i = k B T i n i , p e = k B T e n e , (2.5)
with constant temperatures T i and T e , and we also assume that the electric field E derives from a potential φ, i.e. E = -∇ x φ. Consequently, the model (2.3)-(2.4) is reduced to

               ∂ t n + ∇ x • q α = 0 , ∂ t q α + ∇ x • q α ⊗ q α n + k B T α m α ∇ x n = q α e m α (-n ∇ x φ + q α × B) , α ∈ {i, e} .
(2.6)

We introduce the characteristic length x, time t, momentum q, ion temperature T , electric potential φ, and magnetic field B such that

x = x x ′ , t = t t ′ , T i = T , T e = T T ′ e , (2.7
)

n(x x ′ , t t ′ ) = n n ′ (x ′ , t ′ ) , φ(x x ′ , t t ′ ) = φ φ ′ (x ′ , t ′ ) , q α (x x ′ , t t ′ ) = q q ′ α (x ′ , t ′ ) , B(x x ′ , t t ′ ) = B B ′ (x ′ , t ′ ) .
(2.8)

We consider the natural ratio

q = x n t , (2.9)
and we also assume that the electric and magnetic forces are of the same order, which means in terms of characteristic scales that

qB = n φ x .
(2.10)

We define the characteristic sound speed c for the ions, the characteristic Mach number for the ions M and the characteristic cyclotron frequency ω for the ions by

c = k B T m i , M = q n c , ω = e B m i . (2.11)
Considering a low Mach number regime, we introduce

M = √ τ , (2.12)
with τ ≥ 0 small and assuming that the applied magnetic field is strong allows us to take

t ω = 1 τ . (2.13)
Finally, we denote the ratio m e /m i by ǫ and we assume that it is a dimensionless fixed constant. Then, removing the primed notations and adding the index τ , we finally obtain the rescaled isothermal two-fluid Euler-Lorentz model (2.1).

2.2.

The limit model. We denote by n 0 , q 0 i , q 0 e and φ 0 the limit of (n τ ) τ > 0 , (q τ i ) τ > 0 , (q τ e ) τ > 0 and (φ τ ) τ > 0 respectively when τ → 0. Then, taking the limit regime of (2.1) when τ → 0 indicates that n 0 , q 0 i , q 0 e and φ 0 satisfy the following model

     ∂ t n 0 + ∇ x • q 0 α = 0 , (2.14a)
T α ∇ x n 0 = q α -n 0 ∇ x φ 0 + q 0 α × B , (2.14b) α ∈ {i, e} , (2.14c)
in which the parallel part of q 0 i and q 0 e are implicit. Indeed, if we separate the parallel and perpendicular parts of (2.14b) for any α, the model (2.14) writes

             ∂ t n 0 + ∇ x • q 0 α = 0 , (2.15a)
T α b • ∇ x n 0 = -q α n 0 b • ∇ x φ 0 , (2.15b) (q 0 α ) ⊥ = 1 B b × (q α T α ∇ x n 0 + n 0 ∇ x φ 0 ) , (2.15c) α ∈ {i, e} , (2.15d)
where

q ⊥ = b × (q × b), b = B B and B 2 = B 2 x + B 2 y + B 2
z . We observe in this reformulated limit model that (q 0 i ) ⊥ and (q 0 e ) ⊥ can be algebraically computed from n 0 and φ 0 , but we do not get any explicit constraint for (q 0 i ) || and (q 0 e ) || . One way to answer to this difficulty is to couple the limit model (2.15) with the following equations:

           ∂ t (q 0 α ) || -∂ t (b ⊗ b) q 0 α + (b ⊗ b) ∇ x • q 0 α ⊗ q 0 α n 0 + lim τ → 0 1 ǫ α τ (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ = 0 , α ∈ {i, e} .
(2.16)

These equations are not more than the parallel part of (2.1b) for any α ∈ {i, e}, with τ → 0. However, such a coupling is permitted if

1 ǫ α τ (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ = O(1) , ∀ α ∈ {i, e} , (2.17)
when τ is close to 0. Indeed, this property implies that the limit term within equations (2.16) exists in R and allows to compute (q 0 α ) || . From a numerical point of view, the resolution of (2.15)-(2.16) requires the computation of the limit terms

lim τ → 0 1 ǫ α τ (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ , α ∈ {i, e} . (2.18)
Even if the existence of these limits is admitted, their numerical approximation is not straightforward since n τ and φ τ have to be known for any τ close to 0 and have to satisfy (2.17). As a consequence, a numerical scheme based on a direct resolution of (2.15)-(2.16) requires a numerical resolution of (2.1) for very small values of τ . For these reasons, we choose a numerical approach for computing some approximations of n 0 , q 0 i , q 0 e and φ 0 which differs from a direct resolution of the model (2.15)-(2.16). More precisely, this new approach uses the same tools as in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF][START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF]: we apply a well-chosen time semi-discretization on the initial model (2.1), we reformulate the obtained semi-discrete scheme such that it is a time semi-discrete equivalent to a certain model which is itself equivalent to (2.1). This is what we do in the following sections.

3. Reformulation of the τ -dependent model and of the limit model. In this section, we present a reformulation of the model (2.1) in which the (ǫ α τ ) -1 -order singularities have been regularized. This equivalent model is the following:

                                                 ∂ 2 t n τ - 1 ǫ α τ ∇ x • (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ = ∇ x • ∂ t (b ⊗ b) q τ α -(b ⊗ b) ∇ x • q τ α ⊗ q τ α n τ -∂ t ((q τ α ) ⊥ ) , (3.1a)
∂ t (q τ α ) || -∂ t (b ⊗ b) q τ α + (b ⊗ b) ∇ x • q τ α ⊗ q τ α n τ + 1 ǫ α τ (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ = 0 , (3.1b) (q τ α ) ⊥ = 1 B b × q α T α ∇ x n τ + n τ ∇ x φ τ + q α ǫ α τ B b × ∂ t q τ α + ∇ x • q τ α ⊗ q τ α n τ , (3.1c) α ∈ {i, e} , (3.1d)
These equations are obtained for (2.1) by a differentiation in time and position procedure and projections in the direction of b and perpendicularly to b. More details about these computations can be found in Appendix A. When τ → 0, we obtain

                           ∇ x • (b ⊗ b) T α ∇ x n 0 + q α n 0 ∇ x φ 0 = 0 , (3.2a) ∂ t (q 0 α ) || -∂ t (b ⊗ b) q 0 α + (b ⊗ b) ∇ x • q 0 α ⊗ q 0 α n 0 + lim τ → 0 1 ǫ α τ (b ⊗ b) T α ∇ x n τ + q α n τ ∇ x φ τ = 0 , (3.2b)
(q 0 α ) ⊥ = 1 B b × q α T α ∇ x n 0 + n 0 ∇ x φ 0 , (3.2c) α ∈ {i, e} , (3.2d)
which is equivalent to (2.15)-(2.16). Indeed, (3.2a) can be obtained directly from (2.15b) by applying the differential operator

∇ x • (b •).
4. Semi-discrete AP schemes. In this section, we propose a semi-discretization of (2.1) in time which, up to some space differentiations and combinations, is also a semi-discrete equivalent to the reformulated model (3.1): proceeding in such a way ensures us that the constructed numerical method will be a time semi-discrete equivalent to (3.2) when τ converges to 0 and then a time semi-discrete equivalent to the drift-fluid model (2.14). Another motivation of this methodology is that we want the numerical method to be conservative like finite volume schemes: this kind of property may be not ensured if we build our numerical method on a direct resolution of the reformulated model (3.1).

In practical terms, the semi-discretization we describe in the next lines is based on semi-implicit mass fluxes and fully implicit pressure and Lorentz forces. This strategy is motivated by the fact that we want to preserve the balance between the pressure gradient and the Lorentz term, i.e. to ensure that, at every time step t m ,

T α ∇ x n τ,m + q α n τ,m ∇ x φ τ,m -q τ,m α × B m = O(ǫ α τ ) , ∀ α ∈ {i, e} , (4.1)
where the notation θ m stands for an approximation of the function θ = θ(x, t) at the time step t = t m . This methodology differs from the AP semi-discretizations which were described in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF] and [START_REF] Degond | An asymptotic preserving scheme for the Euler equations in a strong magnetic field[END_REF]: indeed, in these papers, the authors considered a fully explicit mass flux, a fully implicit Lorentz term and a semi-implicit pressure gradient, and this leads to a non-conservative discretization of the velocity equation.

Firstly, we describe the semi-discretization of (2.1) and we reformulate the semidiscrete model by following the same approach as in Section 3. Then we discuss the difficulties which are brought by this reformulation and we introduce a regularization of the mass conservation equations (2.1a) which allows us to bypass these difficulties.

                           n τ,m+1 -n τ,m ∆t + ∇ x • (b m+1 ⊗ b m+1 ) q τ,m+1 α + ∇ x • (I -b m+1 ⊗ b m+1 ) q τ,m α = 0 , (4.2a) q τ,m+1 α -q τ,m α ∆t + ∇ x • q τ,m α ⊗ q τ,m α n τ,m + T α ǫ α τ ∇ x n τ,m+1 = q α ǫ α τ -n τ,m+1 ∇ x φ τ,m+1 + q τ,m+1 α × B m+1 , (4.2b) α ∈ {i, e} . (4.2c)
As it has been announced above, we chose to implicit the whole pressure and Lorentz terms in order to have (4.1) at every time step. The choice of the implicitation of the parallel part of the mass fluxes is motivated by the fact that we have to reformulate the model (4.2) by injecting the parallel part of (4.2b) in (4.2a). Such a procedure leads to the separate computation of (q τ,m+1 i

) m+1 || , (q τ,m+1 i ) m+1 ⊥ , (q τ,m+1 e ) m+1
|| and (q τ,m+1 e ) m+1

⊥ by applying projection operators in the parallel and orthogonal directions to b m+1 . This leads to

                                               (q τ,m+1 α ) m+1 ⊥ - q α ǫ α τ ∆t B m+1 b m+1 × (q τ,m+1 α ) m+1 ⊥ = 1 B m+1 b m+1 × q α T α ∇ x n τ,m+1 + n τ,m+1 ∇ x φ τ,m+1 + q α ǫ α τ B m+1 b m+1 × - q τ,m α ∆t + ∇ x • q τ,m α ⊗ q τ,m α n τ,m , (4.3a)
(q τ,m+1 α

) m+1 || = (b m+1 ⊗ b m+1 ) q τ,m α -∆t (b m+1 ⊗ b m+1 ) ∇ x • q τ,m α ⊗ q τ,m α n τ,m - ∆t ǫ α τ (b m+1 ⊗ b m+1 ) T α ∇ x n τ,m+1 + q α n τ,m+1 ∇ x φ τ,m+1 , (4.3b)
α ∈ {i, e} , (4.3c)
on one hand, and to a couple of anisotropic diffusion equations for n τ,m+1 and φ τ,m+1 on the other hand with an anisotropy carried by b m+1 (see the computations of Appendix B with C i = C e = 0). These diffusion equations are of the form

-∇ x • (b m+1 ⊗ b m+1 ) ∇ x n τ,m+1 + λ τ n τ,m+1 = τ R τ,m+1 , (4.4)
-∇ x • n τ,m+1 (b m+1 ⊗ b m+1 ) ∇ x φ τ,m+1 = τ S τ,m+1 , (4.5)
where λ only depends on ǫ, ∆t and T e , and where

R τ,m+1 = R ∆t, T e , ǫ, n τ,m , q τ,m i , q τ,m e , b m+1 , (4.6) S τ,m+1 = S ∆t, T e , ǫ, n τ,m+1 , n τ,m , q τ,m i , q τ,m e , b m+1 . (4.7)
Remark that the algebraic equations (4.3a) can be solved for any value of τ . Then the scheme (4.2) is equivalent to (4.3)-(4.7) and is Asymptotic-Preserving if and only if

b m+1 • (T α ∇ x n τ,m+1 + q α n τ,m+1 ∇ x φ τ,m+1 ) = O(ǫ α τ ) , ∀ α ∈ {i, e} . (4.8)
These hypotheses are validated if we take into account the following boundary conditions

(b m+1 • ∇ x n τ,m+1 ) (b m+1 • ν) = 0 , on ∂Ω, (4.9a) (b m+1 • ∇ x φ τ,m+1 ) (b m+1 • ν) = 0 , on ∂Ω, (4.9b)
alongwith the diffusion equations (4.4) and (4.5). Indeed, the solutions of the problems (4.4)-(4.9a) and (4.5)-(4.9b) satisfy

b m+1 • ∇ x n τ,m+1 = O(τ ) , b m+1 • ∇ x φ τ,m+1 = O(τ ) , (4.10)
which correspond to (4.8) up to some linear combinations. ) is well posed for any τ > 0. To be more precise, the limit of (4.4)-(4.9a) writes

-∇ x • (b m+1 ⊗ b m+1 ) ∇ x ñ0,m+1 = 0 , on Ω, (b m+1 • ∇ x ñ0,m+1 ) (b m+1 • ν) = 0 , on ∂Ω, (4.11)
and a solution ñ0,m+1 of (4.11) is defined up to a function c : Ω → R satisfying b m+1 • ∇ x c = 0. Since we want to compute the particular solution n 0,m+1 of (4.11) which is exactly the limit of (n τ,m+1 ) τ > 0 when τ → 0, we follow the same approach as in [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF] and we use the following theorem:

Theorem 4.1 (Brull, Degond, Deluzet [START_REF] Brull | Degenerate anisotropic elliptic problems and magnetized plasmas simulations[END_REF]). Let us consider the subspace K ⊂ L 2 (Ω) defined by

K = u ∈ L 2 (Ω) : b m+1 • ∇ x u = 0 , (4.12)
and the functional space W 0 defined by

W 0 = u ∈ L 2 (Ω) : ∇ x • (b m+1 u) ∈ L 2 (Ω) , (b m+1 • ν) u = 0 on ∂Ω , (4.13) provided with the norm u W0 = ∇ x • (b m+1 u) L 2 (Ω)
. Then we have the following properties:

1. K is a closed subset in L 2 (Ω), 2. W 0 is a Hilbert space and ∇ x • (b m+1 W 0 ) is a closed subset of L 2 (Ω), 3. L 2 (Ω) = K ⊕ K ⊥ with K ⊥ = ∇ x • (b m+1 W 0 ).
Having these results in hand and assuming that n τ,m+1 is in L 2 (Ω), we write

n τ,m+1 = π τ,m+1 + q τ,m+1 , (4.14)
with π τ,m+1 ∈ K and q τ,m+1 ∈ K ⊥ . Since the solution n τ,m+1 of (4.4)-(4.9a) is unique when τ > 0, the functions π τ,m+1 and q τ,m+1 are also unique as the projection of n τ,m+1 on K and K ⊥ respectively. Then the diffusion problem (4.4)-(4.9a) writes

     -∇ x • (b m+1 ⊗b m+1 ) ∇ x q τ,m+1 + λ τ (π τ,m+1 + q τ,m+1 ) = τ R τ,m+1 , on Ω, (b m+1 • ∇ x q τ,m+1 ) (b m+1 • ν) = 0 , on ∂Ω. (4.15)
If we consider the variational formulation of (4.15) over K, we find that (λ π τ,m+1 -R τ,m+1 ) ∈ K ⊥ , i.e. there exists h τ,m+1 such that

λ π τ,m+1 -R τ,m+1 = ∇ x • (b m+1 h τ,m+1 ) , on Ω, (b m+1 • ν) h τ,m+1 = 0 , on ∂Ω. (4.16)
By applying the operator b m+1 • ∇ x on this equation, we find that h τ,m+1 is the unique solution of

-b m+1 • ∇ x ∇ x • (b m+1 h τ,m+1 ) = b m+1 • ∇ x R τ,m+1 , on Ω, (b m+1 • ν) h τ,m+1 = 0 , on ∂Ω, (4.17)
which is a well-posed problem for any τ ≥ 0. Since q τ,m+1 ∈ K ⊥ , we claim that there exists l τ,m+1 ∈ L 2 (Ω) such that

q τ,m+1 = ∇ x • (b m+1 l τ,m+1 ) , on Ω, (b m+1 • ν) l τ,m+1 = 0 , on ∂Ω. (4.18)
If we consider now the variational formulation of (4.15) over K ⊥ , we find that l τ,m+1 is the solution of a fourth-order problem which can be written as two successive second-order problems of the form

     -b m+1 • ∇ x ∇ x • (b m+1 L τ,m+1 ) + τ λ L τ,m+1 = -τ b m+1 • ∇ x R τ,m+1 , on Ω, (b m+1 • ν) L τ,m+1 = 0 , on ∂Ω, (4.19)
-b m+1 • ∇ x ∇ x • (b m+1 l τ,m+1 ) = L τ,m+1 , on Ω, (b m+1 • ν) l τ,m+1 = 0 , on ∂Ω. (4.20)
Remark that these problems remain well-posed for any value of τ ≥ 0. Then, instead of solving the problem (4.4)-(4.9a) which becomes ill-posed when τ = 0, we solve the problems (4.17), (4. [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF]) and (4.20) for computing h 0,m+1 and l 0,m+1 , then we compute π 0,m+1 and q 0,m+1 by using (4.16) and (4.18) respectively, and we finally get n 0,m+1 as the sum of π 0,m+1 and q 0,m+1 .

Concerning the problem (4.5)-(4.9b) for the electric potential φ τ,m+1 , we remark that it is ill-posed for any value of τ ≥ 0. Indeed, assuming that this diffusion problem admits at least one solution φτ,m+1 , we can prove that this solution is not unique: for this purpose, we consider a function c

: Ω → R satisfying b m+1 • ∇ x c = 0 , on Ω. (4.21)
Then, it is straightforward to see that φτ,m+1 + c is also a solution of the problem (4.5)-(4.9b). Since this argument works for the problem (4.5)-(4.9b) for any value of τ ≥ 0, the decomposition of its solution does not provide unique projections on K and K ⊥ just as it is done for n τ,m+1 . Then, we have to find a way to restore the uniqueness of the solution of the diffusion problem for the electric potential. This is what we do in the next paragraph.

4.2.

Regularized two-fluid Euler-Lorentz model. As it is explained in the previous paragraph, the Asymptotic-Preserving scheme for the model (2.1) we have contructed leads to a diffusion problem for computing the electric potential at time step t m+1 which admits more than one solution. In order to bypass the difficulty and to restore the well-posedness of the problem in φ τ,m+1 , we choose to include a small regularization in the mass conservation equations (2.1a). That is why we introduce the terms C i ∂ t φ and C e ∂ t φ in such a way that this new model writes

           ∂ t n τ + C α ∂ t φ τ + ∇ x • q τ α = 0 , ǫ α τ ∂ t q τ α + ∇ x • q τ α ⊗ q τ α n τ + T α ∇ x n τ = q α -n τ ∇ x φ τ + q τ α × B , α ∈ {i, e} . (4.22)
Here C i , C e > 0 are two fixed small parameters which will be chosen later. Then we consider the same semi-discretization method as previously, i.e.

                                 n τ,m+1 -n τ,m ∆t + C α φ τ,m+1 -φ τ,m ∆t + ∇ x • (b m+1 ⊗ b m+1 ) q τ,m+1 α + ∇ x • (I -b m+1 ⊗ b m+1 ) q τ,m α = 0 , (4.23a) q τ,m+1 α -q τ,m α ∆t + ∇ x • q τ,m α ⊗ q τ,m α n τ,m + T α ǫ α τ ∇ x n τ,m+1 = q α ǫ α τ -n τ,m+1 ∇ x φ τ,m+1 + q τ,m+1 α × B m+1 , (4.23b)
α ∈ {i, e} . (4.23c)

By splitting parallel and perpendicular parts of (4.23b) for α = i and α = e, we get that (q τ,m+1 α ) m+1

⊥ and (q τ,m+1 α ) m+1

|| satisfy (4.3). Following the same procedure as in the previous paragraph (see Appendix B), we inject (4.3b) with α = i (resp. α = e) in (4.23a) with α = i (resp. α = e). Under the hypotheses

C i + ǫ C e = 0 and C i - ǫ T e C e = C , (4.24)
with C > 0 being given, we find that n τ,m+1 and φ τ,m+1 satisfy the following diffusion equations:

-∇ x • (b m+1 ⊗ b m+1 ) ∇ x n τ,m+1 + τ λ 1 n τ,m+1 = τ R τ,m+1 , (4.25)
-∇ x • n τ,m+1 (b m+1 ⊗ b m+1 ) ∇ x φ τ,m+1 + τ λ 2 φ τ,m+1 = τ S τ,m+1 , (4.26)
where λ 1 , λ 2 only depend on ǫ, ∆t, T e and C, and where

R τ,m+1 = R ∆t, T e , ǫ, n τ,m , q τ,m i , q τ,m e , b m+1 , (4.27)
S τ,m+1 = S ∆t, T e , ǫ, C, n τ,m+1 , n τ,m , φ τ,m , q τ,m i , q τ,m e , b m+1 . (4.28)

Remark that the constraints (4.24) are equivalent to

C i = T e C 1 + T e , C e = - T e C ǫ (1 + T e ) . (4.29)
As a consequence, it is necessary to take C > 0 small enough to ensure that C i and C e are close to 0. For that, we can take C = O(ǫ) provided that the ratio ǫ = m e /m i is small.

We couple (4.25) with the boundary condition given in (4.9a). This diffusion problem is well-posed for any τ > 0 and becomes ill-posed if τ = 0 because of a lack of uniqueness of the solution (see page 9). As a consequence, we can apply Theorem 4.1 and write n τ,m+1 under the following form:

n τ,m+1 = π τ,m+1 + q τ,m+1 , (4.30)
with π τ,m+1 ∈ K, q τ,m+1 ∈ K ⊥ defined by

π τ,m+1 = 1 λ 1 R τ,m+1 + ∇ x • (b m+1 h τ,m+1 ) , q τ,m+1 = ∇ x • (b m+1 l τ,m+1 ) , (4.31)
where h τ,m+1 and l τ,m+1 are the solutions of

-b m+1 • ∇ x ∇ x • (b m+1 h τ,m+1 ) = b m+1 • ∇ x R τ,m+1 , on Ω, (b m+1 • ν) h τ,m+1 = 0 , on ∂Ω, (4.32)
and

     -b m+1 • ∇ x ∇ x • (b m+1 L τ,m+1 ) + τ λ 1 L τ,m+1 = -τ b m+1 • ∇ x R τ,m+1 , on Ω, (b m+1 • ν) L τ,m+1 = 0 , on ∂Ω, (4.33)
-b m+1 • ∇ x ∇ x • (b m+1 l τ,m+1 ) = L τ,m+1 , on Ω, (b m+1 • ν) l τ,m+1 = 0 , on ∂Ω. (4.34)
Now, we couple (4.26) with the boundary condition given in (4.9b). This diffusion problem has the same properties as (4.25)-(4.9a) which has been discussed above. Then, we can write φ τ,m+1 under the following form:

φ τ,m+1 = πτ,m+1 + qτ,m+1 , (4.35)
with πτ,m+1 ∈ K and qτ,m+1 ∈ K ⊥ defined by

πτ,m+1 = 1 λ 2 S τ,m+1 + ∇ x • (b m+1 hτ,m+1 ) , qτ,m+1 = ∇ x • (b m+1 lτ,m+1 ) , (4.36)
where hτ,m+1 and lτ,m+1 are the solutions of

-b m+1 • ∇ x ∇ x • (b m+1 hτ,m+1 ) = b m+1 • ∇ x S τ,m+1 , on Ω, (b m+1 • ν) hτ,m+1 = 0 , on ∂Ω, (4.37)
and

     -b m+1 • ∇ x ∇ x • (b m+1 Lτ,m+1 ) + τ λ 2 Lτ,m+1 = -τ b m+1 • ∇ x S τ,m+1 , on Ω, (b m+1 • ν) Lτ,m+1 = 0 , on ∂Ω, (4.38)
-b m+1 • ∇ x ∇ x • (b m+1 lτ,m+1 ) = Lτ,m+1 , on Ω, (b m+1 • ν) lτ,m+1 = 0 , on ∂Ω. (4.39)
5. Fully-discrete scheme. In this section, we present the fully-discrete version of the Asymptotic-Preserving method for (4.22) we have presented in the previous paragraph. Before going further, we introduce some notations which will be used throughout this section. First, we consider a uniform mesh (x i , y j , z k ) = (i ∆x, j ∆y, k ∆z) on Ω and we define the following subsets of Z 3 :

I = (i, j, k) ∈ Z 3 : (x i , y j , z k ) ∈ Ω , I = {(i + α, j + β, k + γ) : (i, j, k) ∈ I, (α, β, γ) ∈ {-1, 0, 1} 3 . (5.1)
Then, we define the meshed domain Ω h by

Ω h = (i,j,k) ∈ I [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ] × [z k-1/2 , z k+1/2 ] . (5.2)
We also define the subsets I * and I * of Z 3 as follows:

I * = (i, j, k) ∈ Z 3 : (x i+1/2 , y j+1/2 , z k+1/2 ) ∈ Ω h , I * = (i, j, k) ∈ I * : (x i+1/2 , y j+1/2 , z k+1/2 ) / ∈ ∂Ω h . (5.3)
Finally, we assume that, for any K = (K x , K y , K z ) ∈ Z 3 and x K = (x Kx , y Ky , z Kz ), the notation "f |K " stands for an approximation of f (x K ) and that the notation "f |K * " stands for an approximation of f (x K * ). From now, we also denote the point x K = (x Kx , y Ky , z Kz ) by a cell center and the point x K * = (x Kx+1/2 , y Ky+1/2 , z Kz+1/2 ) by a node.

At this point, the choice of the methodology between "discretize first, then reformulate the discretized model" and "reformulate first, then discretize the reformulated model" is critical. Indeed, choosing the second approach induces that we can perform any usual numerical scheme for solving the diffusion equations (4.25)-(4.9a) and (4.26)-(4.9b) through the decomposition method suggested in paragraph 4.1.2 and 4.2, this method will typically give us

n τ,m+1 (x K ) = n τ,m+1 i,j,k + O(h s ) , φ τ,m+1 (x K ) = φ τ,m+1 i,j,k + O(h s ) , (5.4)
with some s > 0 and where h = max(∆x, ∆y, ∆z). Having these results in hands, we inject them into a fully-discrete version of (4.3), we observe that

∆t ǫ α τ (b m+1 (x K ) ⊗ b m+1 (x K )) T α ∇ x n τ,m+1 (x K ) + q α n τ,m+1 (x K ) ∇ x φ τ,m+1 (x K ) = ∆t ǫ α τ (b m+1 ⊗ b m+1 ) T α ∇ h n τ,m+1 + q α n τ,m+1 ∇ h φ τ,m+1 |K + O( ∆t h s ′ ǫ α τ ) , (5.5)
with s ′ > 0, and that we have to strongly refine the space mesh if a very small value of τ is considered. This is not the Asymptotic-Preserving property we expect for. As a consequence, we focus on the second approach, i.e. discretize first then reformulate the discretized model.

The next lines are structured as follows: firstly, we present the finite volume scheme based on the semi-discretization (4.23). Then, we reformulate the obtained fully-discrete scheme by following the same approach as in section 4 and we suggest a numerical method for solving the fully-discrete diffusion equations for n τ,m+1 and φ τ,m+1 . 5.1. Finite volume scheme. First, we introduce some notations for the explicit and implicit fluxes for the hydrodynamic part of (4.23):

f exp,τ,m α,a =     e a • (I -b m+1 ⊗ b m+1 ) q τ,m α q τ,m α,a q τ,m α n τ,m     , (5.6)
f imp,τ,m+1 α,a =     e a • (b m+1 ⊗ b m+1 ) q τ,m+1 α e a T α n τ,m+1 ǫ α τ     , (5.7)
where α ∈ {i, e}, a ∈ {x, y, z}, e a is a vector of the canonical basis of R 3 , q τ,m α = (q τ,m α,x , q τ,m α,y , q τ,m α,z ) and where I is the 3 × 3 identity matrix. Then, we consider different notations for divergence and gradient operator depending on whether they are applied on some component of the implicit or the explicit fluxes. More precisely, we define the operators ∇ F V h •, ∇ h • and ∇ h by linking them to the fluxes by the following relations:

  ∇ F V h • (I -b m+1 ⊗ b m+1 ) q τ,m α ∇ F V h • q τ,m α ⊗ q τ,m α n τ,m   = a ∈ {x,y,z} ∂ a f exp,τ,m α,a , (5.8)
  ∇ h • (b m+1 ⊗ b m+1 ) q τ,m+1 α T α ǫ α τ ∇ h n τ,m+1   = a ∈ {x,y,z} ∂ a f imp,τ,m+1 α,a , (5.9)
where α ∈ {i, e}. At this point, the divergence operators ∇ h • and ∇ F V h • are different: ∇ F V h • is obtained from an finite volume scheme applied on the explicit fluxes whereas ∇ h • is a divergence operator on the implicit fluxes and will be defined differently alongwith the gradient operator ∇ h .

First, we explicit the divergence operator ∇ F V h •. The fully-discrete model obtained from (4.23) is written as follows:

                                       n τ,m+1 |K -n τ,m |K ∆t + ∇ h • (b m+1 ⊗ b m+1 ) q τ,m+1 α |K + ∇ F V h • (I -b m+1 ⊗ b m+1 ) q τ,m α |K = -C α φ τ,m+1 |K -φ τ,m |K ∆t , (5.10a) q τ,m+1 α |K -q τ,m α |K ∆t + ∇ F V h • q τ,m α ⊗ q τ,m α n τ,m |K = - 1 ǫ α τ T α ∇ h n τ,m+1 + q α (n τ,m+1 ∇ h φ τ,m+1 -q τ,m+1 α × B m+1 ) |K , (5.10b)
α ∈ {i, e} .

(5.10c)

In finite volume terms, we have

(∂ a f exp,τ,m α,a ) |K = 1 ∆a F τ,m α,a | K+ea/2 -F τ,m α,a | K-ea/2 , (5.11)
with

F τ,m α,a | K+ea/2 = 1 2 f exp,τ,m α,a |K + f exp,τ,m α,a |K+e a - 1 2 D τ,m α,a | K+ea/2 W τ,m α |K+e a -W τ,m α |K .
(5.12)

where α ∈ {i, e} and a ∈ {x, y, z}. In these formulae, W τ,m α |K is

W τ,m α |K = n τ,m |K q τ,m α |K , (5.13)
and D τ,m α,a is the numerical viscosity matrix linked with the flux f exp,τ,m α,a

for any α ∈ {i, e}. For the sake of simplicity, we choose to compute the viscosity matrices with Rusanov's method (see [START_REF] Rusanov | The calculation of the iteraction of non-stationary shock waves and obstacles[END_REF] and [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF]). For this purpose, we denote the eigenvalues of the jacobian matrices Jac Wα (f exp,τ,m α,a

) by λ τ,m α,k,a (k = 1, 2, 3, 4). Then, D τ,m α,a is defined by

D τ,m α,a | K+ea/2 = I max k=1,2,3,4 max |λ τ,m α,k,a |K+e a | , |λ τ,m α,k,a |K | .
(5.14) 5.2. Reformulation of the fully-discrete scheme. Following the same approach as in sections 3 and 4, we reformulate the discretized model (5.10) by computing separately the perpendicular part and the parallel part of q τ,m+1 i |K and q τ,m+1 e |K and by solving two diffusion equations to find n τ,m+1

|K and φ τ,m+1 |K . More precisely, under the hypotheses (4.24) for C i and C e , we solve

-∇ h • (∇ h n τ,m+1 ) m+1 || |K + λ 1 τ n τ,m+1 |K = τ R τ,m+1 |K , (5.15)
for finding n τ,m+1 |K and

-∇ h • n τ,m+1 (∇ h φ τ,m+1 ) m+1 || |K + λ 2 τ φ τ,m+1 |K = τ S τ,m+1 |K , (5.16)
for finding φ τ,m+1 |K . These discrete diffusion equations are obtained by injecting the parallel part of (5.10b) according to b m+1 |K into (5.10a), then performing some linear combinations of the obtained equations (see Appendix C).

Having n τ,m+1 |K and φ τ,m+1 |K in hand, we can compute separately the parallel part and the perpendicular part of q τ,m+1 i |K and q τ,m+1 e |K by using the following formulae:

                                                 (q τ,m+1 α ) m+1 || |K = (q τ,m α ) m+1 || |K -∆t ∇ F V h • q τ,m α ⊗ q τ,m α n τ,m m+1 || |K - ∆t ǫ α τ T α ∇ h n τ,m+1 + q α n τ,m+1 ∇ h φ τ,m+1 m+1 || |K , (5.17a) (q τ,m+1 i ) m+1 ⊥ |K - q α ǫ α τ ∆t B m+1 |K b m+1 |K × (q τ,m+1 i ) m+1 ⊥ |K = 1 B m+1 b m+1 × (n τ,m+1 ∇ h φ τ,m+1 + q α T α ∇ h n τ,m+1 ) |K + q α ǫ α τ B m+1 |K b m+1 |K × - q τ,m α |K ∆t + ∇ F V h • q τ,m α ⊗ q τ,m α n τ,m |K , (5.17b)
α ∈ {i, e} .

(5.17c) 5.3. Three-point scheme. In this paragraph, we focus on the resolution of (5.15) and (5.16) provided with a discretization of the Neumann-like boundary conditions (4.9). For solving these diffusion equations, we follow the approach of Degond & Tang in [START_REF] Degond | All speed scheme for the low Mach number limit of the isentropic Euler equation[END_REF]: we choose a three-point scheme by replacing the equations (5.15) and (5.16) by

     -(∂ m+1 h, * ∂ m+1 h n τ,m+1 ) |K + τ λ 1 n τ,m+1 |K = τ R τ,m+1 |K , ∀ K ∈ I , (∂ m+1 h n τ,m+1 ) |K * = 0 , ∀ K ∈ I * \I * , (5.18)
and

     -∂ m+1 h, * (n τ,m+1 * ∂ m+1 h φ τ,m+1 ) |K +τ λ 2 φ τ,m+1 |K = τ S τ,m+1 |K , ∀ K ∈ I , (∂ m+1 h φ τ,m+1 ) |K * = 0 , ∀ K ∈ I * \I * , (5.19)
respectively. In these equations, n τ,m+1 * |K * stands for the average of n τ,m+1 on the node (x Kx+1/2 , y Ky+1/2 , z Kz+1/2 ) defined by

n τ,m+1 * |K * = 1 8 α,β,γ ∈ {0,1} n τ,m+1 | K+αex+βey +γez , (5.20)
and the operators ∂ m+1 h and ∂ m+1 h, * correspond to some approximations of b m+1 • ∇ x and ∇ x • (b m+1 •) respectively. These operators are defined by |K and φ τ,m+1 |K that satisfy

(∂ m+1 h p) |K * = b m+1 |K * •             β,γ ∈ {0,1} p | K+ex+βey +γez -p | K+βey+γez 4∆x α,γ ∈ {0,1} p |K+αe x+ey +γez -p |K+αe x+γez 4∆y α,β ∈ {0,1} p | K+αex+βey +ez -p | K+αex+βey 4∆z             , (
(T α ∂ m+1 h n τ,m+1 + q α n τ,m+1 * ∂ m+1 h φ τ,m+1 ) |K * = O(ǫ α τ ) , ∀ α ∈ {i, e} . (5.23)
We remark that we have the good Asymptotic-Preserving property on the nodes. However, we need it on the cell centers, i.e.

b m+1 • (T α ∇ h n τ,m+1 + q α n τ,m+1 ∇ h φ τ,m+1 ) |K = O(ǫ α τ ) , (5.24)
for any α ∈ {i, e}. To reach such a result, we introduce the discrete gradient ∇ h, * defined by

(∇ h, * p) |K * =             β,γ ∈ {0,1} p | K+ex+βey +γez -p | K+βey +γez 4∆x α,γ ∈ {0,1} p |K+αe x+ey +γez -p |K+αe x+γez 4∆y α,β ∈ {0,1} p | K+αex+βey +ez -p | K+αex+βey 4∆z             , (5.25)
and we use it to couple the three-point scheme we have presented with the formulae for (q τ,m+1 i

) m+1 ⊥ |K , (q τ,m+1 e ) m+1 ⊥ |K , (q τ,m+1 i ) m+1 || |K and (q τ,m+1 e ) m+1
|| |K : more precisely, we use the operator ∇ h, * for computing the following terms in (5.17):

                                                   T α ∇ h n τ,m+1 + q α n τ,m+1 ∇ h φ τ,m+1 m+1 || |K = 1 8 α,β,γ ∈ {0,1} (T α ∇ h, * n τ,m+1 + q α n τ,m+1 * ∇ h, * φ τ,m+1 ) m+1 || | K * -αex -βey -γez , (5.26a) 1 B m+1 b m+1 × (q α n τ,m+1 ∇ h φ τ,m+1 + T α ∇ h n τ,m+1 ) |K = 1 8 α,β,γ ∈ {0,1} 1 B m+1 b m+1 × n τ,m+1 * ∇ h, * φ τ,m+1 + ∇ h, * n τ,m+1 | K * -αex-βey -γez , (5.26b) α ∈ {i, e} . (5.26c)
As a consequence, we obtain the properties (5.15) and (5.16) on cell centers and we can compute the parallel part and the perpendicular part of q τ,m+1 i and q τ,m+1 e by using separately the formulae (5.17).

6. Numerical results. In this last section, we present some 2D numerical results which have been obtained with the AP scheme we have presented in sections 4.2 and 5 for the perturbed two-fluid Euler-Lorentz model (4.22).

6.1. Validation of the three-point scheme for the diffusion problems. Since the AP scheme relies on the truthfulness of the properties (4.8), we first present some numerical results from the three-point scheme used for the resolution of the diffusion problems for n τ,m+1 and φ τ,m+1 (see paragraph 5.3). More precisely, the main goal of the first test sequence is to solve the diffusion problems (4.25)-(4.9a) and (4.26)-(4.9b) for any value of τ ≥ 0 and to ensure that

∀ K , n τ,m+1 |K → n 0,m+1 |K , φ τ,m+1 |K → φ 0,m+1 |K , (6.1)
as τ → 0. To perform this validation, we apply our method to the following diffusion problem:

-∇ x • H τ (b ⊗ b) ∇ x p τ ) + τ λ p τ = τ f τ , on Ω, H τ (b ⊗ b) • ∇ x p τ • ν = 0 , on ∂Ω, (6.2)
where λ > 0,

f τ : Ω → R, H τ : Ω → R * + and b : Ω → R 3 are given.
Let us consider a function sequence (p τ ) τ ≥ 0 defined by

p τ = p 0 + τ p τ 1 , (6.3)
with p 0 and p τ 1 satisfying b • ∇ x p 0 = 0 , on Ω, (6.4) and H τ (b ⊗ b) ∇ x p τ 1 • ν = 0 , on ∂Ω. (6.5) We assume from now that

f τ = λ p τ -∇ x • H τ (b ⊗ b) ∇ x p τ 1 , (6.6)
which implies that p τ is the analytic solution of the problem (6.2).

In Figures 1 and2, we plot the evolution of the error between p τ and its approximation (denoted with p τ app ) as a function of the space step h. In these results which are presented in decimal logarithmic scale, we have chosen τ = 10 -2 and τ = 10 -9 . The domain Ω is set to [START_REF] Beer | Toroidal gyrofluid equations for simulations of tokamak turbulence[END_REF][START_REF] Belaouar | emphAn asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF]

× [1, 2] ⊂ R 2 , p 0 and p τ 1 to p 0 (x, y) = 2 , p 1 (x, y) = (x -1)(2 -x)(y -1)(2 -y) 3 , (6.7)
and λ, H τ , and b are chosen as λ = 1 , (6.8) H τ (x, y) = 1 + sin 2 (x) sin 2 (y) , (6.9) b = (sin θ, -cos θ) , with θ(x, y) = arctan(y/x).

(6.10) As we can remark in theses figures, the solver for the diffusion problem (6.2) based on the micro-macro decomposition presented in Section 4.1 and on the discrete differential operators ∂ h and ∂ h, * is second order accurate in h since we observe that the error p τ -p τ app L p (p = 1, 2, ∞) is linearly decreasing in log 10 scale when h → 0, with a slope which is equal to 2. This is due to the fact that, according to the definitions (5.21) and (5.22), the operators ∂ h and ∂ h, * are themselves second -

Error in log scale

Space step in log scale

Error on p in log scale

Max norm L1 norm L2 norm Figure 2. L 1 , L 2
and L ∞ norms of the error between p τ and its approximation p τ app as functions of h: case with τ = 10 -9 .

order accurate. Furthermore, we have this property for τ = 10 -2 and τ = 10 -9 , so we can conclude that the second order accuracy of the solver is not penalized by the smallness of τ . Together with this convergence results in h, we plot in Figure 3 the error between p τ app provided by the solver and p 0 = p 0 as a function of τ with a 100 × 100 uniform mesh, and we take the same values of p 0 , p τ 1 , λ, H τ , and b as above. By definition of the analytic solution of p τ and p 0 (see (6.3)), we except this error to be of the same order of τ . This is confirmed by Figure 3: indeed, the error p τ app -p 0 L p (p = 1, 2, ∞) is linearly decreasing in log 10 scale when τ converges to 0 with a slope which is equal to 1.

From these two results, we can claim that lim

h → 0 lim τ → 0 p τ app = p 0 , (6.11)
which is exactly to say that the numerical solver for the diffusion problem (6.2) is Asymptotic-Preserving when τ → 0. Then, we can use it for solving the diffusion problems (4.25)-(4.9a) and (4.26)-(4.9b) for n τ,m+1 and φ τ,m+1 .

6.2. Validation of the AP scheme for the two-fluid Euler-Lorentz model near the drift-fluid limit. In order to validate the AP scheme we have developed for the perturbed Euler-Lorentz problem (4.22), we compare the results which are computed by the AP scheme to those which can be produced with a fully explicit finite volume method. From now, we denote with classical method a finite volume scheme which is based on the following time semi-discretization:

                           n τ,m+1 -n τ,m ∆t + C α φ τ,m+1 -φ τ,m ∆t + ∇ x • q τ,m α = 0 , q τ,m+1 α -q τ,m α ∆t + ∇ x • q τ,m α ⊗ q τ,m α n τ,m + T α ǫ α τ ∇ x n τ,m = q α - 1 ǫ α τ n τ,m ∇ x φ τ,m + q τ,m+1 α × B m+1 , α ∈ {i, e} . (6.12)
We easily remark that the stability condition of such a method strongly depends on τ : as in a low Mach number numerical experiment, the smaller τ is, the smaller the time step ∆t must be in order to ensure that a method based on (6.12) is stable.

As a example, if we solve the hydrodynamic part of (6.12) with Rusanov' scheme, we must have ∆t = O(h τ 1/2 ) at least, where h = min(∆x, ∆y, ∆z).

In the next lines, we distinguish two opposite situations:

• The resolved case: The time step is small enough in order to ensure that both classical and AP methods capture the fast time variations within the solution, • The under-resolved case: The time step does not allow the capture of fast time variations but ensures at least the stability of the AP scheme. The test case we present here is based on the perturbation of the following stationary case:

• The magnetic field is uniform and reads B = (sin α, -cos α, 0) with α ∈ R fixed, • n τ,0 (x, y) = n 0 , φ τ,0 (x, y) = φ 0 , q τ,0 i = q τ,0 e = B with some constants n 0 and φ 0 , • The whole system (4.22) does not depend on the variable z. q τ i,x (AP scheme) q τ i,x (classical scheme) q τ i,y (AP scheme) q τ i,y (classical scheme)

Figure 4. Resolved case at time t = 6×10 -6 : x and y components of the ion momentum q τ i as functions of (x, y) computed with the AP scheme (left) and the classical scheme (right).

Remarking that both classical and AP schemes compute the exact solution provided with these initial datas, we choose to introduce a small perturbation at the initial time step. More precisely, we choose to replace n τ,0 (x, y) = n 0 by n τ,0 (x, y) = n 0 + τ max 0, 1 -η (x -x 0 ) 2 -η (y -y 0 ) 2 , (6.13) with η ≥ 0 and (x 0 , y 0 ) ∈ Ω.

q τ e,x (AP scheme) q τ e,x (classical scheme)

q τ e,y (AP scheme) q τ e,y (classical scheme)

Figure 5. Resolved case at time t = 6×10 -6 : x and y components of the electron momentum q τ e as functions of (x, y) computed with the AP scheme (left) and the classical scheme (right).

We also assume that the physical domain Ω is [START_REF] Beer | Toroidal gyrofluid equations for simulations of tokamak turbulence[END_REF][START_REF] Belaouar | emphAn asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF] × [START_REF] Beer | Toroidal gyrofluid equations for simulations of tokamak turbulence[END_REF][START_REF] Belaouar | emphAn asymptotically stable semi-lagrangian scheme in the quasi-neutral limit[END_REF] and is meshed by a 100 × 100 uniform mesh. We also precise the initial datas by taking τ = 10 -8 , ǫ = 1, T e = 3, C = 10 -2 , α = 2π 3 , η = 80, (x 0 , y 0 ) = ( 3 2 , 3 2 ), n 0 = 1 and φ 0 = 0. In Figures 45, we present some results in the resolved situation for both classical and AP schemes and ∆t = 5 × 10 -9 is taken as time step. As we can see in these figures, the results which are produced by the AP scheme are very close to the classical method's ones.

In Figures 67, we present some simulations which are obtained in the underresolved case, i.e. where the time step is taken much larger than the time step which is required to ensure the stability of the classical scheme. In the present case, we have chosen ∆t = 10 -6 , which is 200 times larger than the time step which has been used for the resolved case above. As we can see in these figures, the AP scheme remains stable and produces the same results as in the resolved case. However, the classical scheme blows up after a small number of time iterations, which is not surprising because the time step we have chosen is too large for satisfying the stability condition of this scheme.

From this numerical experiment, we can conclude that the AP scheme we have developed for the perturbed two-fluid Euler-Lorentz model (4.22) allows us to take a time step which does not satisfy the stability condition required for capturing the fast time variations of the solution. Furthermore, such a time step choice does not penalize the quality of the results which are obtained with the AP scheme. q τ i,x (AP scheme) q τ i,x (classical scheme) q τ i,y (AP scheme) q τ i,y (classical scheme) Figure 6. Under-resolved case at time t = 6 × 10 -6 : x and y components of the ion momentum q τ i as functions of (x, y) computed with the AP scheme (left) and the classical scheme (right).

q τ e,x (AP scheme) q τ e,x (classical scheme)

q τ e,y (AP scheme) q τ e,y (classical scheme)

Figure 7. Under-resolved case at time t = 6 × 10 -6 : x and y components of the electron momentum q τ e as functions of (x, y) computed with the AP scheme (left) and the classical scheme (right). 6.3. Impact of the perturbation parameter C. Since the AP scheme we have built for the resolution of the perturbed Euler-Lorentz model (4.22) has been validated in terms of quality of results, the impact of the value of C needs to be investigated. Indeed, we recall that this parameter is linked with the constants C i and C e by the relations

C i = T e C 1 + T e , C e = - T e C ǫ (1 + T e ) , (6.14)
and these constants are introduced to recover the uniqueness of the solution for the diffusion problem (4.5)-(4.9b) (see Section 4). Since C i and C e are introduced in (2.1) through a perturbation of the mass conservation equations, these constants are assumed to be as close to 0 as possible. However, this is equivalent to assume that C is close to 0, so the diffusion problem In the last test sequence, we run the AP method with the initial datas which have been used in the previous paragraph, i.e.

• We take a 100 × 100 uniform mesh over

Ω = [1, 2] × [1, 2],
• The magnetic field is uniform and is defined as B = (sin α, -cos α, 0) with α = 2π 3 , • The initial electric potential is φ τ,0 = φ 0 with φ 0 = 0, • The initial ion momentum and the initial electron momentum are defined by q τ,0 i = q τ,0 e = B, • The initial density is n τ,0 defined by

n τ,0 (x, y) = n 0 + τ max 0, 1 -η (x -x 0 ) 2 -η (y -y 0 ) 2 , (6.15)
with n 0 = 1 constant, η = 80 and (x 0 , y 0 ) = ( 3 2 , 3 2 ), • We choose T e = 3, ǫ = 1 and τ = 10 -8 . In Figure 8, we plot the x-component of the ion momentum q τ i which is obtained with the AP scheme (4.23) at time t = 6 × 10 -6 with C = 10 -2 and a time step ∆t which is equal to 10 -6 , 10 -7 and 10 -8 . We can remark that the AP method is stable for ∆t ≤ 10 -6 since all the results are similar. Now, we do again this numerical experiment with C = 10 -3 instead of C = 10 -2 .

As we can remark in Figure 9 where q τ i,x is plotted at time t = 4 × 10 -6 , the AP method is stable with ∆t = 10 -7 and ∆t = 10 -8 . However, we observe some important boundary effects in the ∆t = 10 -6 case.

Finally, we perform this numerical experiment with C = 10 -4 . In Figure 10, we plot q τ i,x at time t = 2 × 10 -6 with ∆t = 10 -6 , 10 -7 , 10 -8 . No numerical artifacts are present in the ∆t = 10 -8 case whereas we observe some boundary effects if ∆t = 10 -6 and even the blowing up of the method when ∆t = 10 -7 is considered. As a consequence, we can say that the AP scheme we have developed for the perturbed two-fluid Euler-Lorentz model (4.22) is Asymptotic-Preserving when C > 0 is fixed and when τ → 0. However, it is not Asymptotic-Preserving when C → 0 and τ > 0 is fixed.

7. Conclusions and perspectives. In this paper, we studied the isothermal twofluid Euler-Lorentz system coupled with a quasi-neutrality constraint in a low Mach number regime and a strong magnetic field regime. After having presented the model and its limit regime, we proposed a reformulation of this model which is compatible with the construction of an Asymptotic-Preserving scheme. Then we scheme is recalled here:

                                 n τ,m+1 -n τ,m ∆t + C α φ τ,m+1 -φ τ,m ∆t + ∇ x • (b m+1 ⊗ b m+1 ) q τ,m+1 α + ∇ x • (I -b m+1 ⊗ b m+1 ) q τ,m α = 0 , (B.1a) q τ,m+1 α -q τ,m α ∆t + ∇ x • q τ,m α ⊗ q τ,m α n τ,m + T α ǫ α τ ∇ x n τ,m+1 = q α ǫ α τ -n τ,m+1 ∇ x φ τ,m+1 + q τ,m+1 α × B m+1 , (B.1b) α ∈ {i, e} . (B.1c)
First, we separate the parallel and the perpendicular parts of (B.1b) according to b m+1 . More precisely, we obtain the equations (4.3a) by performing the vector product of b m+1 by (B.1b). Concerning the equations (4.3b), we obtain them by multiplying the tensor (b m+1 ⊗ b m+1 ) by (B.1b).

In order to obtain the diffusion equations (4.25) and (4.26), we put (4.3b) in (B.1a). We obtain a system of two diffusion equations satisfied by (n τ,m+1 , φ τ,m+1 ):

                           n τ,m+1 -n τ,m ∆t + C α φ τ,m+1 -φ τ,m ∆t + ∇ x • q τ,m α -∆t ∇ x • (b m+1 ⊗ b m+1 ) ∇ x • q τ,m α ⊗ q τ,m α n τ,m - ∆t ǫ α τ ∇ x • (b m+1 ⊗ b m+1 ) T α ∇ x n τ,m+1
+ q α n τ,m+1 ∇ x φ τ,m+1 = 0 , α ∈ {i, e} . ) .

(B.8)

Appendix C. Reformulation of the fully-discrete problem. In this paragraph, we develop the reformulation procedure for the finite volume scheme which is detailed in Section 5. This scheme writes

                               n τ,m+1 |K -n τ,m |K ∆t + ∇ h • (b m+1 ⊗ b m+1 ) q τ,m+1 α |K + ∇ F V h • (I -b m+1 ⊗ b m+1 ) q τ,m α |K = -C α φ τ,m+1 |K -φ τ,m |K ∆t , (C.1a) q τ,m+1 α |K -q τ,m α |K ∆t + ∇ F V h • q τ,m α ⊗ q τ,m α n τ,m |K = - 1 ǫ α τ
T α ∇ h n τ,m+1 + q α (n τ,m+1 ∇ h φ τ,m+1q τ,m+1 |K and (C.1b) and we obtain (5.17b). In order to obtain the discrete diffusion equations (5.15) and (5.16), we follow the same procedure as in the semi-discrete case (see Appendix B): we replace (q τ,m+1 i ) m+1 || (resp. (q τ,m+1 e ) m+1

||

) by its expression given by (5.17a) with α = i (resp. α = e) and we obtain two diffusion equation of the form

                           -∇ h • (T α ∇ h n τ,m+1 + q α n τ,m+1 ∇ h φ τ,m+1 ) m+1 || |K + ǫ α τ ∆t 2 (n τ,m+1 + C α φ τ,m+1 ) = ǫ α τ ∇ h • - 1 ∆t q τ,m α + ∇ F V h • ( q τ,m α ⊗ q τ,m α n τ,m ) m+1 || + 1 ∆t 2 n τ,m + C α ∆t 2 φ τ,m - 1 ∆t ∇ F V h • (q τ,m α ) m+1 ⊥ |K , α ∈ {i, e} . (C.2)
Finally, we consider the constraints (4.24) to make these 2 equations uncoupled and, up to some linear combinations, they can be rewritten under the form

+ ∇ F V h • ( q τ,m i ⊗ q τ,m i n τ,m ) + ǫ ∇ F V h • ( q τ,m e ⊗ q τ,m e n τ,m ) m+1 || - 1 ∆t ∇ F V h • (q τ,m i ) m+1 ⊥ + ǫ ∇ F V h • (q τ,m e ) m+1
+ ∇ h • - 1 ∆t (q τ,m i - ǫ T e q τ,m e ) + ∇ F V h • ( q τ,m i ⊗ q τ,m i n τ,m ) - ǫ T e ∇ F V h • ( q τ,m e ⊗ q τ,m e n τ,m ) m+1 || - 1 ∆t ∇ F V h • (q τ,m i ) m+1 ⊥ - ǫ T e ∇ F V h • (q τ,m e ) m+1
⊥ |K .

(C.5)

4. 1 .

 1 Time semi-discretization.

4. 1 . 1 .

 11 Asymptotic-Preserving property. The considered time semi-discretization of (2.1) is the following:

  K * -βey -γez -(b m+1 x p) | K * -ex -βey -γez K * -αex -βey -(b m+1 z p) | K * -αex-βey -ez 4∆z . (5.22) Then, replacing b m+1 •∇ x and ∇ x •(b m+1 •) by ∂ m+1 h and ∂ m+1 h, * in the decomposition procedure (4.30)-(4.39), we compute n τ,m+1

Figure 1 .

 1 Figure 1. L 1 , L 2 and L ∞ norms of the error between p τ and its approximation p τ app as functions of h: case with τ = 10 -2 .

20 S

 20 . BRULL, P. DEGOND, F. DELUZET AND A. MOUTON -

Figure 3 .

 3 Figure 3. L 1 , L 2 and L ∞ norms of the error between p τ and p 0 = p 0 as functions of τ : case with a 100 × 100 uniform mesh.

7 ∆t = 10 - 8 Figure 8

 71088 Figure 8. q τ i,x at time t = 6 × 10 -6 with C = 10 -2 and ∆t = 10 -6 , 10 -7 , 10 -8 .

∆t = 10 - 6 ∆t = 10 - 7 ∆t = 10 - 8 Figure 9

 1061071089 Figure 9. q τ i,x at time t = 4 × 10 -6 with C = 10 -3 and ∆t = 10 -6 , 10 -7 , 10 -8 . The results which are provided by the whole test sequence above indicate that we have a stability condition for the AP scheme which clearly depends on C and which would be of the form ∆t = O(C) . (6.16)

(B. 2 )-ǫ ∆t 2 4 )e ∆t 2

 2242 By doing some linear combinations, these diffusion equations write-∇ x • (b m+1 ⊗ b m+1 ) ∇ x n τ,m+1 + τ 1 + ǫ ∆t 2 (1 + T e ) n τ,m+1 + τ C i + ǫ C e ∆t 2 (1 + T e ) ∆t 2 n τ,m + C i + ǫ C e ∆t 2 φ τ,m + ∇ x • (b m+1 ⊗ b m+1 ) ∇ x • ( q τ,m i ⊗ q τ,m i n τ,m ) + ǫ ∇ x • ( q τ,m e ⊗ q τ,m e n τ,m) ,(B.3) and -∇ x • n τ,m+1 (b m+1 ⊗ b m+1 ) ∇ x φ τ,m+1 + τ T e C i -ǫ C e ∆t 2 (T e -1) φ τ,m+1 + τ T e -ǫ ∆t 2 (T e -1) T e n τ,m + T ǫ C i -ǫ C e ∆t 2 T e φ τ,m + ∇ x • (b m+1 ⊗ b m+1 ) ∇ x • (If we consider the constraints (4.24) for C i and C e , we decouple the diffusion equations. Firstly, we compute n τ,m+1 by solving-∇ x • (b m+1 ⊗ b m+1 ) ∇ x n τ,m+1 + τ λ 1 n τ,m+1 = τ R τ,m+1 , ∆t 2 n τ,m + ∇ x • (b m+1 ⊗ b m+1 ) ∇ x • ( it to compute φ τ,m+1 by solving -∇ x • n τ,m+1 (b m+1 ⊗ b m+1 ) ∇ x φ τ,m+1 + τ λ 2 φ τ,m+1 = τ S τ,m+1 , T e (n τ,m+1 -n τ,m ) + C ∆t 2 φ τ,m + ∇ x • (b m+1 ⊗ b m+1 ) ∇ x • ( q τ,m i ⊗ q τ,m i n τ,m ) -ǫ T e ∇ x • ( q τ,m e ⊗ q τ,m e n τ,m

i×

  B m+1 ) |K , (C.1b) α ∈ {i, e} . (C.1c) As in the previous appendices, we separate the parallel part and the perpendicular part of (C.1b) according to b m+1 |K . By multiplying the tensor (b m+1 |K ⊗ b m+1 |K ) by (C.1b), we obtain (5.17a) for each α. Concerning (q τ,m+1 α ) m+1 ⊥ |K , we compute the vector product of b m+1

-

  ∇ h • (∇ h n τ,m+1 ) m+1 || |K + λ 1 τ n τ,m+1 |K = τ R τ,m+1 |K , (C.3)and-∇ h • n τ,m+1 (∇ h φ τ,m+1 ) m+1 || |K + λ 2 τ φ τ,m+1 |K = τ S τ,m+1 |K , (C.4) with λ 1 , λ 2 , R τ,m+1|K and S τ,m+1 |K defined byλ 1 = 1 + ǫ ∆t 2 (1 + T e ), λ 2 = T e C ∆t 2 (T e -

  e -1 ǫ -T e ∆t 2 T e (n τ,m+1 -n τ,m ) + C ∆t 2 φ τ,m

Acknowledgments. This work has been supported by the french magnetic fusion programme FR-FCM, by the CEA-Cadarache in the frame of the contract 'AP-PLA' (# V3629.001 av. 2), by the INRIA large-scale initiative 'FUSION', by the University Paul Sabatier in the frame of the contract 'MOSITER' and by the Fondation 'Sciences et Technologies pour l'Aéronautique et l'Espace', in the frame of the project 'Plasmax' (contract # RTRA-STAE/2007/PF/002). The authors would like to thank the referees for their careful reading and helpful remarks.

supported by NSF grant xx-xxxx.

presented the time semi-discrete AP scheme itself and its reformulation leading to the resolution of some anisotropic diffusion equations for n τ,m+1 and φ τ,m+1 . The equation for φ τ,m+1 being ill-posed, we restored the uniqueness of the solution of this equation by introducing a small perturbation in the mass conservation equations. Finally, we performed some numerical tests of this scheme by comparing the results from the AP scheme to those from a fully explicit method and we tested the influence of the perturbation parameter C on the behaviour of the AP scheme. At this point, several work pathes can be investigated. The first one is to go back to the time semi-discretization described in section 4.1 by invoking boundary conditions for the computation of the electric potential which are different from Neumann conditions. The second one is to generalize the present work to some two-fluid Euler-Lorentz models involving other pressure laws for p i and p e .

Appendix A. Reformulation of the Euler-Lorentz model. In this paragraph, we detail the reformulation procedure of the Euler-Lorentz model

leading to the model (3.1). First, we separate the parallel and perpendicular parts of (A.1b) for each value of α: we obtain

In order to obtain (3.1a) for any α ∈ {i, e}, we compute the divergence in space of (A.2b) on one hand and the derivative in time of (A.2a) on the other hand. We obtain

and

for any α ∈ {i, e}. By doing some linear combinations of these equations, we obtain a system of 2 equations for n τ and φ τ of the form

which are exactly the equations (3.1a). Then the Euler-Lorentz model (A.1) is equivalent to the combination of (A.2b), (A.2c) and (A.5).

Appendix B. Reformulation of the semi-discrete problem. This paragraph is devoted to the reformulation of the semi-discrete problems (4.2) and (4.22). Since the model (4.2) is not more than (4.22) with C i = C e = 0, we present the reformulation of the semi-discrete scheme for the perturbed Euler-Lorentz model. This


Alice Baudot

Joseph Moysan

Cédric Payan

Noëllie Ylla

Jocelyne Galy

Blandine Verneret

André Baillard

Improving adhesion strength analysis by the combination of ultrasonic and mechanical tests on single lap joints

Keywords: Ultrasonics, lap shear, epoxy, adhesion

The aim of this paper is to improve knowledge about adhesion strength in single-lap joints by comparison of data from both Non Destructive Testing (NDT) and standard mechanical tests. study presents a high-frequency ultrasonic imaging of aluminum-epoxy-aluminum joints that carefully chemically prepared using two different surface treatments of the aluminum substrate: silanization and chromic-sulfuric etching. The surface energy of the treated substrate is complete set of ultrasonic data are obtained using high-frequency focused transducers. After inspection, single-lap shear tests were performed, and the mechanical strength of each joint is determined. The good cor relations among the surface measurements, the mechanical tests and ultrasonic results are presented and discussed. The demonstration of a close correlation between ultrasonic images obtained at high frequency and the fracture surface is promising.

Introduction

Structural bonding for the automotive industry, the naval industry and aeronautics was one of the major industrial changes in the field of material assembly [1]. This major change was accompanied by the necessity of developing methodologies to assess the strength of a bonded assembly. Assessing this strength remained a real technological and scientific challenge for many years because of regular changes in the materials that were used, as well as changing standards.

This challenge led to many interdisciplinary studies, as it involved both material science and mechanical science. The question of how an assembly should be mechanically tested has produced dozens of different mechanical tests. Concurrently, finite-element modeling offered new possibilities for understanding the effective mechanical strength of a structure and the impacts of many parameters (e.g., thickness, roughness, and stress concentration) and for attempting to predict cohesive or adhesive rupture [2]. Structural bonding created specific challenges for NDT science, as well. As bonding assemblies most of the time the result of an innovative process, material parts often possessed new properties and required the adaptation of NDT methods or the development of new ones [3]. The of adhesive bonding in primary structures has been limited by the lack of non-destructive testing procedures capable of guaranteeing the reliability of the joint. To date, no reliable NDT test for adhesion strength of a bond has been developed. At present, this problem is addressed by strict of the adherend-surface-preparation procedures [3]. In comparison with classical NDT methods previously developed for weld testing, the NDT testing of bonded assemblies presents new as material parts can have very different thicknesses due to the finer thickness of the bond layer.

Ultrasonic and acoustic methods are the approaches that have been most commonly pursued early stages of structural bonding applications [3].

The field of NDT applications is continuously expanding; they are used for assessing the initial manufacturing quality of the bonding process, for in-service testing and for assessing the quality of bonding subject to aging degradation. Three classical defects are investigated [3,4]: a) complete voids, disbonds or porosity b) poor adhesion (a weak bond between the adhesive and one or both adherends) c) poor cohesive strength (a weak adhesive layer)

This study aims to propose advances in NDT techniques for the poor-adhesion case. Poor adhesion is a reduction in the strength of the bond between the adhesive and the adherend. It is extremely difficult to detect using classical ultrasonic techniques because the bond strength is governed by a layer thinner than a conventional ultrasonic wavelength. By preference, a highfrequency transducer should be used for signal analysis [5]. Severe cases in which the poor adhesion tends to zero are commonly termed kissing bonds or, alternatively, zero-volume disbonds [6]. This form of poor adhesion stands between the poor-adhesion case and the wholly disbonded case. Numerous ultrasonic testing solutions, have been developed for this specific case including non-linear acoustics [5,[7][8][9] but also guided waves [10][11]. Alternative solutions are also in development, such as Digital Image Correlation [12]. One of the obstacles to the development and demonstration of an acoustic method capable of quantifying the adhesive strength of a bonded joint is the absence of perfect reference samples full knowledge of their assembly. Multiple solutions have been proposed in the past, and still are currently to create reference samples for NDT testing. Defects can be artificially created by foreign part at the interface between two parts [13] or by modifying the bond itself by release agent [14] or by incorporating particles [15][16]. To study the adhesion strength, the most appropriate approach would be to vary the surface preparation or the surface treatment [17]. parameters can modify the bond strength, so the true difficulty is to exactly reproduce the same preparation steps in the same environment [18]. More rarely, some studies use a chemical obtain imperfections similar to those created by pollution during industrial processes [19]. way to produce variable adhesion is to compress a bond; this case is representative of drykissing bonds [7].

The single-lap shear joint is the most widely studied type of adhesive joint in the literature. The manufacture of this joint is relatively easy, and the corresponding test requires a classic tensile testing machine. Many corresponding standards have been developed in industry, such as the ASTM D1002, the popular standard test method for evaluating the apparent shear strength of single-lap-joint adhesively bonded metal specimens by tension loading, or the European Standard EN 2243-1. Adhesive-lap-joint shear-strength tests do not create perfect shear solicitation, and other physical factors are uncertain, such as the overlap length and the bondline thickness [17]. Stress concentrations should also be taken into account for accurate mechanical analysis. The form of the adhesive edges greatly influences the shear-stress concentrations. Appropriate geometries were proposed to improve mechanical analysis by reducing edge effects [20]. This joint geometry is chosen for our study because it also enables a simple NDT evaluation using through-transmission ultrasonic testing. In the following, the experiments are described, and the results are presented and discussed in detail. Conclusions are proposed regarding the evaluation of variable adhesion on single-lap shear joints using high-frequency ultrasonic analysis.

Experimental methods

Materials

The adhesive formulation is based on a two-component epoxy system. Diglycidyl ether of bisphenol A (DGEBA, DER 332, Dow Chemicals, Paris Saint Denis, France) with an epoxy equivalent weight of 171-175 g/eq is cross-linked with a cycloaliphatic amine hardener, isophorone diamine (IPD, Aldrich, Saint-Quentin Fallavier, France). The epoxy and the amine are mixed at room temperature, at a stoichiometric ratio aminohydrogene/epoxy equal to 1; the mixture is degassed to avoid air-bubble formation in the adhesive bond.

glycidoxypropyltrimethoxysilane ( -GPS) is obtained from Aldrich (Saint-Quentin Fallavier, France).

The chemical structures of both components are given in Table 1.

The metallic substrate is an aluminum alloy, A2024, in the form of 1.6 mm thick sheets.

Sample preparation

The aluminum sheets are cut to 100 mm x 25 mm, and their surfaces are cleaned and treated before bonding to eliminate surface contamination and to promote different adhesion strengths.

The following two surface treatments are employed prior to the adhesive bonding of the lapshear test samples:

* Silane treatment (Si). There are four steps in this treatment: i) solvent degreasing with a tissue soaked in isopropanol, ii) alkaline cleaning by immersing the specimens for 10 min in a Sococlean (Socomore, Vannes, France) A3431 solution (c = 10 %) heated to 45 °C followed by rinsing in water, iii) alkaline etching by immersing the specimens for 2 min in a solution of 40 g/l NaOH + 10 g/l sodium glucoheptonate heated to 40 °C followed by rinsing in water and iv) silanization by immersing the specimens for 10 min in a 1 wt% solution of GPS in distilled water; prior to immersion, the pH of the solution is adjusted to 5 using acetic acid, and the solution is stirred at room temperature for 1 hour for the hydrolysis of the methoxy groups. The silane-treated specimens are heated for 1 hour at 93 °C in an oven to allow condensation reactions [21]. * Chromic-sulfuric acid treatment (CS). There are three steps in this treatment i) degreasing the specimens with acetone, ii) immersing them in a water solution (~300 g/l H2SO4, ~30 g/l CrO3, ~0.5 g Cu, ~10 g Al) heated at 65 °C for 15 min, and iii) rinsing them in distilled water and drying them in an oven for 30 min at 60 °C.

After the surface treatment, the aluminum pieces are assembled into single-lap shear joints with overlap length of 12.5 mm. The lap-shear geometry is illustrated in Fig. 1. A specific metallic designed to prepare the samples with a given thickness. This mold allows the simultaneous of seven specimens. The mold is placed under a press, and a contact pressure of 2.2 bars is then the samples are cured for 1 hour at 80 °C, followed by curing for 2 hours at 180 °C. This cycle fully cured the DGEBA-IPD adhesive.

Surface energy

The surface energies of the aluminum substrate with and without surface treatment are determined from contact angles (CA) using a GBX Digidrop contact-angle goniometer following the Owens-Wendt method [22]. The average CAs are obtained by measuring the contact angles of at least six drops each for two probe liquids, water and diiodomethane.

Mechanical testing

The samples are tested following EN 2243-1, using a universal testing machine (MTS,Eden Prairie, Minessota,USA) at a constant crosshead speed of 0.1 mm/min. The joint adhesive strength is calculated using the following formula:

= (1)
where is the adhesive strength in MPa, P max is the maximum load at fracture in Newton, and A is the average cross-sectional area in mm 2 . Lap-shear-strength data are obtained, and the corresponding standard deviations (std) are calculated. The full bandwidth is very large (from 0 to 30 MHz), and as a consequence, we obtain short signals with a duration of two periods. This enables the separation of the two surface echoes, but the epoxy layer is too thin to create separated echoes (Fig. 4 and Fig. 6).

NDT characterization

A complete set of ultrasonic signals with both transmitted and reflected signals is obtained. The C-scan images (X, Z cartography) obtained in the transmission mode are the most informative for the discussion presented in section 4. The displacement steps along both the X and Z axes for the C-scan are 0.2 mm, allowing the complete recovery of the focal width.

Results and discussion

3.1. Effect of the surface treatment on the surface energy of the substrate

The contact angles and the calculated surface energies are reported in Table 2 for the asreceived substrates as well as the silane-treated and sulfochromic-treated aluminum substrates.

The surface energy of the as-received A2024 has been found to be equal to 50 mJ/m 2 . This value is in agreement with the data published in the literature [23]. The surface energy of the GPSsilanized aluminum surface is similar to that of the as-received aluminum. Moreover, the polar and non-polar contributions are the same. However, the surface is modified because the protocol used in this work has proven to be effective in causing the condensation of silanol groups on the Al surface [24]. Finally, the sulfochromic treatment leads to a higher surface energy, equal to 72 mJ/m 2 ; this increase is mainly attributed to the increase of the polar component.

Adhesive properties

Protocols for the two surface treatments and the preparation of the lap-shear specimens have validated using preliminary tests. The optimization of the physical-chemistry protocol has been achieved over 28 samples. The present study focuses on the samples (four samples from the silanization set (Si) and three samples from the chromic-sulphuric set (CS)) that have been both ultrasonic C-scan and mechanical tests. The results from the mechanical measurements are summarized in Table 3 for the two adherend treatments, silanization and chromic-sulfuric best adhesive strength has been obtained for the chromic-sulfuric surface treatment. At the end test, joint failure has been observed. In all cases, the failure mode is adhesive, sometimes with more jumps of the fracture from one interface to the other or with some adhesive remaining on the aluminum substrates.

Ultrasonic testing results

Each adhesive joint sample tested using the transmission method produces a C-scan of the bonded area. Our first goal is to verify the quality of the joint. Indeed, before making a prediction regarding the adhesion strength, it is important to assess the joint homogeneity. Cscan images of the samples Si-1 and CS-2 are presented in Fig. 3. and Fig. 5.

A first qualitative observation can be done as images textures are very different. For the sample Si-1 (Fig. 3) large homogeneous areas with little differences in amplitudes are observed. In these large areas some localized defects are also present probably due to large gas entrapment as our bond is too liquid. Gas entrapment is a true difficulty for manufacturing adhesive joints and resulting voids should be taken into account for modes of failures [25]. The epoxy is degassed during 210 seconds using UIP1000hd device (Hielscher, Teltow, Germany) [26]. For the CS-2 sample the image texture, corresponding to ultrasonic amplitude distribution, is more heterogeneous and the spatial scale of the distribution is lower. The texture aspect could can be created by small voids distributed on all the joint surface. The adherend surface roughness contributes also to the texture aspect of C-scan images.

As the transmission method is used, if the joint is properly bonded, the received signal should high amplitude; in the opposite case, if there is a defect involving the presence of air, the the signal should be very low. A-scan data has been recorded at various locations (X, Z) for each sample. The obtained results are presented in Fig. 4 and Fig. 6, which correspond to samples one the silane set (Si-1) and the second of the chromic-sulfuric set (CS-2), respectively. The first echo corresponds to the direct travel of the ultrasound through the first aluminum plate, the joint and second aluminum plate. The second echo corresponds to the same travel path plus a round trip aluminum.

All images have been obtained with the same transducer settings. In area 1 of sample Si-1, the amplitude is very low (dark blue), which allows us to define the amplitude threshold below which the joint is considered not bonded. The corresponding threshold is 0.45 mV; thus, we can calculate the size of the bonded area and the size of the disbonded one, if any. The C-scan analyses have been conducted with MATLAB. The surface size of the adhesive joint is first visually defined. This surface is represented in Fig. 4 and 6 by the bold rectangle. Next, the size of the bonded area is automatically calculated depending on the threshold value. Thus, for sample Si-1, the surface of the joint is 234.6 mm², the bonded area is found to be 211.3 mm², and the disbonded zone is 23.2 mm². For sample CS-2, the surface of the adhesive joint is 273.7 mm², and there is no disbonded surface. Considering the geometry of the specimens (Fig. 1), the theoretical bonded area should be equal to 312.5 mm 2 ; however, the epoxy formulation used has a very low viscosity, and some leaks may occur in the joint. These leaks are the origin of the smaller-than-expected measured bonded area and of the observed defects.

Correlation between mechanical tests and NDT

This section presents the correlation between the mechanical tests and the ultrasonic tests. To for a link between the mechanical strength and the bonded surface area, a plot representing the surface area as a function of the failure load is shown in Fig. 7. Two different populations are set CS has higher failure strengths and greater bonded surface areas than set Si. Not simple linear relation exists between the bonded area and the mechanical strength [3]. It is in a single-lap shear test, the eccentricity of the tensile axis results in peel and transverse stress distribution is highly non-uniform, so the strength is much more sensitive to the integrity some areas of the joint than to others. Therefore, measurement of the bond area, without taking localization of this bonded area into account, does not necessarily provide a linear correlation strength. Literature indicates that the presence of large defects may not modify the ultimate there are localized in the middle of the joint. This behavior depends also of the ductility of the adherend [27].

With the knowledge of the real bonded area, we can evaluate a more accurate joint adhesive strength using formula (1). Table 4 presents the corrected strengths. It reveals that the silane treatment may yield the same adhesion strength as the chromic-sulfuric acid treatment, as the average strengths are very similar. It also reveals that in this study, the silane treatment produces less reproducible results. Such information is not evident from only the results obtained with the mechanical tests (see Table 3). Using strength results, no significant difference arises between the two surface treatments. It is due to the presence of voids identified through the ultrasonic data set. Moreover, additional interesting correlations between C-scan images and the fracture surface are described below.

The mechanical tests continue until failure of the samples. Therefore, at the end of these tests, the fracture surfaces can be optically observed (Fig. 8 and9). A strong similarity in shape can be observed between the visible features of the fracture surfaces and the ultrasonic C-scans. Traces of epoxy are indicated with the letter E on the fracture photographs (Fig. The adhesive failure mode is the most commonly observed in Fig. 8 and9. It is the expected failure mode, as the surface treatment is the variable parameter in this study. Also as expected, the cohesive failure mode is not observed. The third case, the presence of one or more defects, is apparent in zone 1 (defined in Fig. 3 above) of sample Si-1, which is consistent with our analysis of signal 1: there is no ultrasound transmission. This case is not expected and is attributed to an excessively fluid epoxy adhesive.

The position of the defect in zone 1, where there is no stress concentration, could explain why this sample demonstrates a similar strength to the other samples during the shear test despite the presence of a large defect.

The pictures of the fracture surfaces and the C-scans corresponding to samples Si-1 and CS-2 are presented in Fig. 8 and Fig. 9. The bonded area is more homogeneous for sample CS-2. The small, circular, weakly adhesive areas that are visible on the C-scans (darker areas of the image) are also apparent on the fracture surfaces.

Conclusion and prospects

This study shows how the combination of mechanical and ultrasonic NDT methods can provide additional information in the view of a better understanding of adhesion strength. For the material aspect of this challenge, the use of different surface treatment is the preferred solution.

Using a high-frequency C-scan, it is possible to obtain a corrected value for the bonded surface area, allowing the proper evaluation of the adhesive strength. It is also shown that these ultrasonic images exhibit a notably good correlation with the fracture surfaces obtained after failure. Efforts in progress aim at improving the C-scan method using higher frequencies as well as the epoxy adhesive which will be more viscous. It should allow separating echoes in ultrasonic data, increasing the information available such as speed of sound and attenuation. A complementary approach using C-scan images associated with a detailed analysis of the stress concentrations would greatly help multidisciplinary studies concerning adhesion strength.

  Each sample is tested by the immersion technique using two focused Panametrics V376 NDT, Waltham, MA, USA) transducers with a diameter of 12.7 mm and a focal distance of The experiment is performed in a water tank with two motorized arms. Each arm has five axes. The experimental set-up is piloted by the software UTWIN (Mistras, Paris, France). Fig.2presents the experimental set-up for the normal-incidence transmission mode. The transducers' theoretical central frequency is 30 MHz. They have been characterized with a hydrophone to the real acoustic field. The real -3 dB bandwidth ranges from 14 MHz to 22 MHz. The measured diameter of the focal zone is 0.63 mm.

  8 and 9), while apparent aluminum surfaces are indicated with Al. Various indications of likely failure modes are thus obtained: -adhesive failure when the photographs give Al-E or E-Al indications -cohesive failure when the photographs give E-E indications -defect presence when the photographs give Al-Al indications

Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :

 1234 Fig. 1 : Geometry and dimensions of lap joint

Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :

 56789 Fig. 5 : C-scan of sample CS-2 (amplitude in milivolts)

Table 1 : Chemical structure of monomers

 1

	Name	Chemical structure	Molar mass (g/mol)
	DGEBA		348
	IPD		170
	-GPS		236

Table 2 : Surface energy of the substrates

 2

			Mean
	Final Sets	Surface		Number of samples
			(MPa) / std (MPa)
		treatment
	Si	Silanization	9.2 / 0.7	4
	CS	Chromic-	12.6 / 0.7	3
		sulfuric

Table 3 :

 3 Mechanical results

			Mean	Mean	Number of
	Final Sets	Surface treatment
			Bonded area	(MPa) / std (MPa)	samples
	Si	Silanization	(mm²) 222.3 / 51.5	13 / 3.1	4
	CS	Chromic-sulfuric	281.5 / 8.20	14.0 / 1	3

Table 4 : Mechanical results corrected for the bonded area measured with the UT through- transmission method

 4

Acknowledgments

The authors wish to thank Erwin Taviot (Aix-Marseille Université) for his contribution to the ultrasonic experiments, and Professor R. Adams for useful discussions and comments about this

The collaboration with Aircelle and Safran Composites is gratefully acknowledged. This work was supported under the PRC Composites, French research project funded by DGAC, involving Group, ONERA and CNRS.