Florian Bernard

Angelo Iollo

Gabriella Puppo

Simulation of Particle Dynamics for Rarefied Flows: Backflow in Thruster Plumes

A novel model and a new numerical method are presented for the transport of solid particles in rarefied flows. The model is based on a Vlasov type equation where the particles are represented by a distribution function. The rarefied flow is described by a BGK or ES-BGK approach. An accurate method is proposed to solve the particle transport equation in a fully Eulerian framework. Validations in 2D with respect to analytical solutions and a Lagrangian method are presented. The numerical model is then used to explain a peculiar particle dynamics observed in satellite thrusters.

Introduction

In many complex applications, fluid flows can contain solid particles, bubbles or droplets. These applications go from combustion in engines [START_REF] Doisneau | Eulerian multi-fluid models for the simulation of dynamics and coalescence of particles in solid propellant combustion[END_REF] to transport of particle pollutants in atmosphere or rivers. Numerical methods have been developed in the 70s and 80s to deal with this problem by coupling a macroscopic model for the fluid flow with a transport equation with possible source terms for the second phase usually called spray [START_REF] Dukowicz | A particle-fluid numerical model for liquid sprays[END_REF], [START_REF] Harlow | Numerical calculation of multiphase fluid flow[END_REF], [START_REF] Carrillo | Simulation of fluid and particles flows: Asymptotic preserving schemes for bubbling and flowing regimes[END_REF]. The simulation of particle transport in rarefied flows is also of interest but has been less addressed in the literature. In [START_REF] Charles | Modélisation mathématique et étude numérique d'un aérosol. Application à la simulation du transport de particules de poussière en cas d'accident de perte de vide dans ITER[END_REF] the author is interested in the transport of dust particles during a loss-of-vacuum accident in ITER. Other models have been developed in [START_REF] Østmo | Kinetic theory study of steady evaporation from a spherical condensed phase containing inert solid particles[END_REF] for particle transport in rarefied flows. We can also cite the work of Ferrari and Pareschi [START_REF] Ferrari | Modelling and numerical methods for the dynamics of impurities in a gas[END_REF] where the authors deal with diffusion of impurities in granular flows. Here, particle transport models are used in order to explain the peculiar phenomenon of contamination of optical devices carried by satellites due to incompletely burnt particles coming from thrusters (see figure 1). As it has been noticed by Dettleff et al. [START_REF] Dettleff | Basics of Plume Impingement Analysis for Small Chemical and Cold Gas Thrusters[END_REF], the firing of satellite thrusters in rarefied environment pollutes (or damages) a collar located in front of the nozzle (see figure 12 in [START_REF] Dettleff | Basics of Plume Impingement Analysis for Small Chemical and Cold Gas Thrusters[END_REF]). This means that incompletely burnt particles are ejected from the nozzle and turn back towards the collar. On a real set up, this collar represents mirrors, lenses or solar panels that become unusable. Thus, the comprehension of this particular phenomenon is of interest for the aerospace industry to protect optical devices from dust. Up to our knowledge, it has not yet been dealt from a numerical point of view.

For flows outside the planetary atmosphere, the continuum model is no longer valid. The molecular nature of the fluid becomes predominant.

In this case, the dynamics of each gas molecules has to be considered and a statistical approach is more suitable. In this sense, the Boltzmann equation [START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF] is used:

∂f ∂t (x, ξ, t) + ξ • ∇ x f (x, ξ, t) = Q(f, f ) (1)
where f is the mass density distribution function of the gas depending on velocity, space and time. x ∈ R D , x = (x, y, z) in 3D, D being the number of space dimensions, ξ ∈ R d is the microscopic velocity, ξ = (ξ u , ξ v , ξ w ) in 3D, d being the number of energy degrees of freedom), t is the time and the initial condition is f 0 = f (x, ξ, t = 0). Q is a bilinear operator called collision term that represents the interaction between the gas particles. The parameter that dictates whether or not the flow is rarefied is the Knudsen number Kn. It is the ratio of the mean free path between the particles λ and the characteristic length of the problem L. Different methods exist to solve the Boltzmann equation such as DSMC [START_REF] Bird | Molecular Gas Dynamics and the Direct Simulation of Gas Flows[END_REF] or deterministic schemes. The BGK model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] and the ES-BGK model [START_REF] Holway | Kinetic Theory of Shock Structure Using an Ellipsoidal Distribution Function[END_REF] are approximation of the Boltzmann equation viable for a large range of Knudsen numbers (< 1).

In the following these two models are used for the simulation of the rarefied flow. As they are costly to compute, a Cartesian grid is used to take advantage of its adequacy for massive parallel computation. It is also very convenient for the simulation of moving bodies since no remeshing step is required. Now that the model for the gas flow is determined, a model for the particle dynamic is added to the governing equations, to simulate the contamination around satellite thrusters by incompletely burned particles, with an interaction term between the gas flow and the particles.

Even if we consider that the solid particle flow (bubble, incompletely burned particles,...) is diluted, computing the motion of each solid particles would be computationally prohibitive. A statistical approach is then more suitable. Hence, the solid particle flow is described by a Vlasov type equation:

df p dt + ξ • ∇ x f p + ∇ ξ • Ff p = 0 (2)
where f p is the mass density distribution function of the solid particles depending on the space position x ∈ R D , the microscopic velocity ξ ∈ R d (ξ = (ξ u , ξ v , ξ w ) in 3D), the time t, with initial condition f p 0 = f p (x, ξ , t = 0). F is the total force acting on the solid particles. A similar model for the particle was employed in [START_REF] Bernard-Champmartin | A Semi-Lagrangian Approach for Dilute Non-Collisional Fluid-Particle Flows[END_REF] where the carrier flow is in the continuum regime. In [START_REF] Andrews | The multiphase particle-in-cell (MP-PIC) method for dense particulate flows[END_REF] and [START_REF] Snider | Sediment flow in inclined vessels calculated using a multiphase particle-in-cell model for dense particle flows[END_REF], the particle flow is also modelled with a Vlasov equation. However, the fluid equation is averaged taking into account the different phases with an interaction through a drag force. Non-linear models for this drag force are used in [START_REF] Patankar | Modeling and numerical simulation of particulate flows by the Eulerian-Lagrangian approach[END_REF]. In [START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF] another kind of coupling is performed based on moment closure equations. The authors present a method which seems to be efficient for non-equilibrium rarefied flows. The Vlasov equation 2 is usually solved with a particle method to avoid computations where there are no particles. Hence, ideally, computations are performed only where it is needed. Particle methods were initially introduced to simulate simple flows [START_REF] Rosenhead | The Formation of Vortices from a Surface of Discontinuity[END_REF]. Recently, we distinguish three main methods to solve particle motion all based on Lagrangian or semi-Lagrangian schemes. Particle-In-Cell method [START_REF] Evans | The Particle-in-cell Method for Hydrodynamic Calculatons[END_REF] is the most popular to solve the Vlasov equation. This method uses a grid and particles. The particles are moved according to the equation and the Eulerian field is recovered by representing each particle with an interpolation kernel. However, although the method is considered robust and a small number of particles yields satisfactory results, it generates significant numerical noise that could pollute the solution. Moreover, storing particles and a grid increases the memory requirement and could be prohibitive with a rarefied model.

Another kind of particle method also considers regularization of the particles with an interpolation kernel. The pioneering idea has been developed almost simultaneously by Gingold and Monaghan in [START_REF] Gingold | Smoothed particle hydrodynamics -Theory and application to non-spherical stars[END_REF] and by Lucy in [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF]. They are called Smooth Particle Hydrodynamic (SPH) methods. A more recent review is found in [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF]. Instead of PIC method, no grids are considered and usually, the interpolation kernels used are more accurate. They are very efficient for front tracking and free surface motion [START_REF] Cottet | A level set method for fluidstructure interactions with immersed surfaces[END_REF]. However, recovering macroscopic quantities induces errors due to interpolation kernels. In addition, the use of high order interpolation kernels could lead to negative values of the distribution function which are non physical.

Finally, we mention here a class of particle methods for Euler or Navier-Stokes equations. These methods are based on a vorticity formulation of the equation [START_REF] Cottet | Vortex Methods: Theory and Practice[END_REF] and they are particularly effective for incompressible flows. The drawback of these methods is that singularities in the flow can lead to non accurate or non physical solutions [START_REF] Beale | Vortex methods. I. Convergence in three dimensions[END_REF], [START_REF] Cottet | Vortex Methods: Theory and Practice[END_REF].

Remeshing techniques [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] have been developed for vortex methods to improve the approximation when vortices overlap or get too close [START_REF] Cottet | Particle methods revisited: a class of high order finite-difference methods[END_REF] by redistributing them on a grid. In particle flows, remeshing becomes also necessary in very smooth flows since the particles distribution may become non uniform in space. An additional constraint in our case, is the constraint on the preservation of positivity of the distribution function (as in PIC methods). We also want the function to remain on a Cartesian mesh in phase space. This allows also an easy integration of the distribution function to recover the number and density of the particles.

In the following, after a brief presentation of the rarefied models used, the particle dynamics model based on Vlasov equation is introduced. We then present a classical approach based on a fully Lagrangian scheme and a new method based on a remeshing step to solve this model. Two dimensional test cases are presented to validate the method. We finally investigate a realistic test case of a nozzle plume ejecting particles (previously described) which represents the main motivation of this work. Up to our knowledge it is the first attempt to simulate this phenomenon.

Governing Equations

In this section the two rarefied models used to simulate rarefied flows are briefly presented for mono-atomic and mono-species gases. Then, we detail the model describing the particle transport on which we will focus later for nozzle plumes.

BGK and ES-BGK models

We consider two rarefied models to simulate rarefied gas flows. For moderate Knudsen numbers (< 10 -2 ), the BGK model [START_REF] Bhatnagar | A Model for Collision Processes in Gases. I. Small Amplitude Processes in Charged and Neutral One-Component Systems[END_REF] is a good compromise between accuracy and computational cost. In dimensionless form, the BGK equation reads as:

∂f ∂t + ξ • ∇ x f = 1 τ M f -f (3)
with M the equilibrium distribution function in dimensionless form computed from macroscopic quantities. It is expressed as:

M f (x, ξ, t) = ρ(x, t) (2πT (x, t)) 3/2 exp - |ξ -U(x, t)| 2 2T (x, t)
The relaxation time τ depends on the local variables, the reference viscosity µ 0 at the reference temperature T 0 , the reference density ρ 0 , the specific gas constant R and the characteristic length of the problem L:

1 τ = √ RT 0 ρ 0 L µ 0 ρT 1-δ = 1 Kn ∞ ρT 1-δ ( 4
)
where δ is the exponent of the viscosity law of the gas. Kn ∞ is the Knudsen number in reference conditions. The macroscopic quantities (ρ the density, U the velocity and E the total energy) characterizing the flow can be recovered from the moments of f :

                 ρ = R 3 f dξ ρU = R 3 ξf dξ E = R 3 |ξ| 2 2 f dξ (5)
One of the main drawbacks of this model is that the transport coefficients are not correct. In particular, the Chapmann-Enskog expansion for the BGK model gives a Prandtl number of 1 instead of 2/3 for a monoatomic gas. A popular model that fixes this problem is the ES-BGK model [START_REF] Holway | Kinetic Theory of Shock Structure Using an Ellipsoidal Distribution Function[END_REF] which corrects the stress tensor to fix the Prandtl number [START_REF] Andries | Numerical comparison between the Boltzmann and ES-BGK models for rarefied gases[END_REF]. It is quite similar to the BGK model, only the equilibrium function differs:

∂f ∂t + ξ • ∇ x f = 1 τ G f -f (6)
The equilibrium distribution function is a Gaussian calculated as follows:

G f (x, ξ, t) = ρ(x, t) det(2πT (x, t))) exp - (ξ -U(x, t))T -1 (ξ -U(x, t)) T 2
As for the BGK model, density, velocity and energy can be recovered from eq.( 5). In the case of the ES-BGK model, we also need to get the pressure tensor Θ and heat flux q:

       ρΘ = R 3 c ⊗ cf dξ q = R 3 1 2 c|c| 2 f dξ (7)
where c = ξ -U is the peculiar velocity. The symmetric tensor T can be defined as:

T (x, t) = 1 P r T (x, t)I + (1 - 1 P r )Θ(x, t)
where I is the identity matrix. Then we can prove that:

ρT (x, t) = R 3 c ⊗ cGdξ
The relaxation time for the ES-BGK model can be expressed similarly to the one of the BGK model:

1 τ = √ RT 0 ρ 0 L µ 0 P rρT 1-δ = P r Kn ∞ ρT 1-δ (8)
with P r = µcp κ the Prandtl number, µ the viscosity, c p the specific heat and κ the thermal conductivity. In the following, we will always consider δ = 1 for both models.

Particle transport

The idea is to introduce a model representing a spray of particles in the diluted gas flow. Let's consider a set of particles represented by the distribution function f p with microscopic velocities ξ passively transported by the fluid. We assume that the flow of these particles is so diluted that they do not collide between each other. Each particle moves with its own velocity. This velocity is modified by the drag force due to the gas flow. The particle flow can be then modelled by the Vlasov type equation:

∂f p ∂t + ∇ x • ξ f p + ∇ ξ • af p = 0 ( 9
)
The acceleration is due to the drag force [START_REF] Julia | The effects of linear and quadratic drag on falling spheres: an undergraduate laboratory[END_REF]. Here we choose a = D(U f (x)ξ ). In [START_REF] Boudin | Fluid-particles flows: a thin spray model with energy exchanges[END_REF], an analogy to granular flows validated experimentally this approach.

However, the following work does not depend on this specific expression of the drag force. It can be extended to more complex models. We do not consider other forces acting on the particles, such as gravity for example. U f (x) is the velocity field given by the solution of the rarefied model (BGK or ES-BGK model) for the gas:

U f (x) = R 3 ξf dξ R 3 f dξ
The number of particles N p in a domain Ω can be easily recovered:

N p = Ω ρ p dx = Ω R 3 f p dξ dx ( 10
)
where ρ p is the particle density.

Note that in this model we disregard particle collisions. We also suppose that the particle flow is so diluted that it has negligible impact on the gas flow. Thus, no feedback on the rarefied equation for the gas is present in this model that can be considered as a first step towards a fully coupled model.In principle, one could add a drag feedback on a macroscopic model like Euler or Navier-Stokes equation in a relatively straightforward manner (see for example [START_REF] Bernard-Champmartin | A Semi-Lagrangian Approach for Dilute Non-Collisional Fluid-Particle Flows[END_REF]). On the other hand, we consider here a driving flow that is rarefied and we should then model how the feedback from the particles to the rarefied flow is distributed in velocity. This issue is left to future investigations.

3 Numerical schemes

Discretization of the BGK and ES-BGK models

This section is devoted to the numerical schemes used to solve both rarefied models. The space discretization is first presented in the case of the BGK model. Then, we will introduce the time discretization. More details can be found in [START_REF] Bernard | Accurate Asymptotic Preserving Boundary Conditions for Kinetic Equations on Cartesian Grids[END_REF] since we are extending that scheme to the present case.

Space discretization

The space discretization is performed on a Cartesian grid. On this type of grid, numerical schemes are simple and can be easily implemented, with massive parallel computation which is convenient as kinetic models are computationally heavy due to the large number of independent variables. The discretization is the same for both rarefied models presented in section 2.1 and is presented in the case of the BGK model in 2D. Let us consider a domain Ω x :

Ω x = i=1..n j=1..m Ω i,j
x where Ω i,j

x represents the cell (i, j) of the computational domain and such that (x i , y j ) are the coordinates of the center of the cell (i, j) and (x i+1/2 , y j ) are the coordinates of the center of the interface between cells (i, j) and (i + 1, j). We have also that ∆x = x i+1/2,j -x i-1/2,j = ∆y = y i,j+1/2 -y i,j-1/2 . On a space cell Ω i,j x , eq.( 3) is integrated with a finite volume method:

∂f i,j ∂t + ξ • ∂Ω i,j x f n ∂Ω i,j x dσ = 1 τ i,j (M fi,j -f i,j ) (11)
where f i,j = 1

|Ω i,j x | Ω i,j x f dxdy and M fi,j = 1 |Ω i,j x | Ω i,j
x M f dxdy. The scheme can be simply written as:

∂f i,j ∂t + 1 ∆x (F i+ 1 2 ,j -F i-1 2 ,j + F i,j+ 1 2 -F i,j-1 2 ) = 1 τ i,j (M fi,j -f i,j ) (12)
The flux F i+ 1 2 ,j is now expressed as (with a similar notation for the other fluxes):

F i+ 1 2 ,j = max(0, ξ u )f i,j + min(0, ξ u )f i+1,j (13)
where ξ u is the first component of the microscopic velocity and has to be replaced by ξ v , the second component of the microscopic velocity, to compute the fluxes along the second direction y. This is a first order scheme in space but can be easily extended to second order with slope reconstruction and MinMod limiters for example.

Time discretization

The time discretization can be performed for all terms explicitly. But in this case, the time step will be determined by the space discretization (∆x), the maximum velocity of the velocity grid and the relaxation time τ . For small Knudsen numbers, the relaxation part becomes very stiff (τ very small) and imposes a very strong restriction on the time step. Asher et al. [START_REF] Ascher | Implicitexplicit Runge-Kutta methods for time-dependent partial differential equations[END_REF] first presented IMEX schemes to cure this issue. Here, the IMEX scheme [START_REF] Christopher | Additive Runge-Kutta schemes for convection-diffusion-reaction equations[END_REF], [START_REF] Pareschi | Implicit-Explicit Runge-Kutta Schemes and Applications to Hyperbolic Systems with Relaxation[END_REF] is chosen. The relaxation term is treated implicitly while the convective part is non stiff and linear which means that an explicit scheme is more efficient.

The time integration for a ν-stages IMEX Runge-Kutta scheme applied to the BGK model reads as follows [START_REF] Hairer | The Numerical Solution of Differential-Algebraic Systems by Runge-Kutta Methods[END_REF]:

f n+1 i,j = f n i,j -∆t ν k=1 ωk ξ∇ x f (k) i,j + ∆t τ i,j ν k=1 ω k (M (k) fi,j -f (k) i,j ) f (k) i,j = f n i,j -∆t k-1 l=1 Ãk,l ξ∇ x f (l) i,j + ∆t τ i,j k l=1 A k,l (M (l) fi,j -f (l) i,j ) f (1) i,j = f n i,j + ∆t τ i,j A 1,1 (M (1) fi,j -f (1) i,j ) (14)
where A and à are ν × ν matrices, with Ãi,s = 0 if s ≥ i and A i,s = 0 if s > i, ω and ω are vectors of size ν and are the weights for the Runge-Kutta scheme between the solutions at different intermediate steps. These coefficients are derived from a double Butcher's tableaux (see [START_REF] Pieraccini | Implicit-Explicit schemes for BGK kinetic equations[END_REF] for the coefficient values):

à ωT A ω T
All the quantities until stage k -1 are known so the equation for stage k becomes:

f (k) i,j = τ A k,k ∆t + τ f n i,j -∆t k-1 l=1 Ãk,l ξ∇ x f (l) i,j + ∆t τ k-1 l=1 A k,l (M (l) fi,j -f (l) i,j )+ A k,k ∆t τ M (k) fi,j (15)
We are interested in a first order positivity-preserving IMEX scheme given by the following tableaux:

0 1 1 1
We need to distinguish the case of the BGK model and the case of ES-BGK model. In the first case, f (k) i,j can be computed explicitly since all the right hand side is known. Indeed, since the moments of the relaxation term are zero, the macroscopic variables at stage k can be computed integrating in velocity space the second equation of [START_REF] Andrews | The multiphase particle-in-cell (MP-PIC) method for dense particulate flows[END_REF], see [START_REF] Pieraccini | Implicit-Explicit schemes for BGK kinetic equations[END_REF]. Hence, the Maxwellian M fi,j and the distribution function are not strictly the same. In particular, the third moment does not give the same tensor. The trick used for the BGK model cannot be applied for the ES-BGK model. However, Filbet et al. [START_REF] Filbet | An Asymptotic Preserving Scheme for the ES-BGK Model of the Boltzmann Equation[END_REF], Alaia [START_REF] Alaia | A Time Dependent Domain Decomposition Method for a Multiscale Hydrodynamic-Kinetic System of Equations[END_REF] showed that the IMEX scheme can still be applied to the ES-BGK model. If we consider the second equation of ( 14), the three first moments ρ

(k) i,j , U (k) i,j , T (k) i,j can still be obtained explicitly. But to define G (k) fi,j , one also needs Θ (k) i,j . Let us define the tensor Σ (k) i,j : Σ (k) i,j = Ω ξ ξ ⊗ ξf (k) i,j dξ = ρ (k) i,j (Θ (k) i,j + U (k) i,j ⊗ U (k) i,j ) (16)
If [START_REF] Snider | Sediment flow in inclined vessels calculated using a multiphase particle-in-cell model for dense particle flows[END_REF] in the case of the ES-BGK model is multiplied by ξ ⊗ξ and integrated we get:

Σ (k) i,j = τ i,j P r A k,k ∆t + τ i,j P r Σ n i,j -∆t k-1 l=1 Ω ξ ξ ⊗ ξ Ãk,l ξ • ∇ x f (l) i,j + 1 τ i,j A k,l (G (l) fi,j -f (l) i,j ) dξ + A k,k ∆t τ i,j P r + A k,k ∆t ρ (k) i,j (T (k) i,j I + U (k) i,j ⊗ U (k) i,j ) Σ (k)
i,j can be then calculated explicitly and Θ

(k) i,j is deduced. Finally, G (k)
fi,j can be computed and so f (k) i,j .

Resolution of the particle transport

In this section, two methods to solve the Vlasov equation for the particle transport are presented. The first way to treat the transport of the particles is fully Lagrangian that will be used in the numerical results section to be compared with. In a second part we present an other method consisting in remapping the particles on a fixed mesh.

Lagrangian scheme

The first approach considers a Lagrangian scheme to solve [START_REF] Cercignani | The Boltzmann Equation and Its Applications[END_REF]. Each particle has its own trajectory, depends on the other particles and is considered as a Dirac:

f p (x, ξ , t) = p m p δ x (x -x p (t))δ ξ (ξ -ξ p (t)) ( 17
)
where m p is the mass of particle p and δ x (respectively δ ξ ) is the Dirac function in the physical space (respectively the velocity space). This method is particularly efficient for passive transport and is easily parallelizable. A splitting is performed to solve first the transport in physical space and then the transport in velocity space. Eq.( 9) becomes:

     ∂f p ∂t + ∇ x • ξ f p = 0 ∂f p ∂t + ∇ ξ • af p = 0 (18)
These equations are solved by tracking the position of the set of particles initially defined by the distribution function in phase space:

             f p (x, ξ , t) = f p (x 0 , ξ 0 , t = 0) dx dt = ξ dξ dt = D(U f (x) -ξ ) (19a) (19b) (19c)
If the velocity U f (x) is not known analytically (which is usually the case if it comes from the resolution of another equation), it has to be interpolated at the particle position.

One drawback of this model is that no information is exchanged between the particles. Then particles can overlap and lead to a degradation of accuracy or non-physical phenomena [START_REF] Krasny | A study of singularity formation in a vortex sheet by the point-vortex approximation[END_REF]. Moreover, the structure of the grid is lost and it is almost impossible to imagine a feedback from the particles to the rarefied models. Indeed, this feedback would be imposed through the rarefied equation in each cell. A grid structure (identical as in the rarefied equation if possible) is then required for the particle transport to recover the Eulerian field.

In the following we propose another technique based on methods widely used in fluid dynamics.

A particle method with remeshing

The underlying idea of the particle method using a remeshing step is to keep all flows information on the initial mesh, in our case, the Cartesian grid on which the rarefied model for the gas is solved. The scheme is the following. A Lagrangian step is performed at each time step. Then the particles are redistributed in phase space, in each Cartesian cell, according to an interpolation kernel. New equivalent particles are created in the center of each cell while the old ones are suppressed. Hence, at each time step, all data is known in the initial mesh, in the Cartesian cells. In other word, at the beginning of each time step, x i = x p and ξ i = ξ p where i is the cell index and p denotes the particle itself. To avoid a remeshing step in 6D (3D+3D) a splitting is performed between physical space and velocity space. Particles are first transported in physical space at velocity ξ and remeshed in the same space. Then, the transport in velocity space is realized according to the acceleration term. Finally, particles are redistributed in velocity space. At each time step, the equations to solve are:

                 dx p dt = ξ i fp (x i , ξ i ) = R x (f n p (x p , ξ i )) dξ p dt = a p f n+1 p (x i , ξ i ) = R ξ ( fp (x i , ξ p )) (20)
where R x and R ξ are respectively the remeshing operator in space and in velocity. It applies the remeshing kernel on all particles. Figure 2 shows the two different cases for a 1D configuration of the transport of a particle p initially in x = x i . The particle moves to the position x = x p and is remeshed on the two closest grid points (in the case of a two point interpolation kernel). α p and β p are the weights associated to these two grid points for the remeshing of the particle p.

x i-1 x i x i+1 x p α p β p (a) Case ξ u > 0 x i-1 x i x i+1 x p α p β p (b) Case ξ u < 0
For the remeshing step, different kernels are present in the literature preserving the moments of the distribution function up to a certain order (see for example [START_REF] Lagaert | Particle method: an efficient tool for direct numerical simulations of a high Schmidt number passive scalar in turbulent flow[END_REF]). However, in our case positivity and diffusivity properties of the kernel are a real issue. In particular, we need the distribution function of the particles to stay positive after the remeshing step. Moreover, the remeshing stencil has to be as compact as possible to avoid the spreading of the particles due to a wide stencil. Finally, if the position of the particle does not change, the weight associated to the collocated grid point has to be one such that it does not induce numerical diffusion. The best compromise is found with the Λ 1 kernel which is applied in space and in velocity:

Λ 1 = α = 1 -y β = y ( 21
)
where y is defined from figure 2 as:

     y = x p -x i ∆x if ξ u > 0 y = x p -x i-1 ∆x if ξ u ≤ 0
Then, the distribution function after remeshing (still in a general 1D case) is:

f p i = n ω ni f p n ( 22
)
where f p i is the particle distribution function after the remeshing step in x = x i , f p n is the particle distribution function after the transport in x = x n and ω ni are the remeshing weights in x = x i for the particle n. They are defined as:

ω ni = α if x p > x i β if x p < x i
One can note that the remeshing preserves the positivity since the weights ω n , defined from α and β, are always positives.

In 2D, the weights are computed with a tensor product. If α 1 and α 2 (respectively β 1 and β 2 ) are the 1D weights (as in ( 21)) in the first direction (respectively in the second direction), the weights associated to the 2 × 2 stencil (see Figure 3) are computed as:

ω n i,j = α i β j (23)
x i,j ω n i,j

x i,j+1

ω n i,j+1
x i+1,j ω n i+1,j

x i+1,j+1

ω n i+1,j+1

x p The particle distribution function in cell (i, j) is then computed similarly to the 1D case as:

f p i,j = n ω n i,j f p n (24)
One advantage of this method is that computations are done only were particles are actually present and have to be remeshed. This is a significant advantage in computational requirements with respect to classical finite volume schemes where computations would be necessary everywhere even when the particles are present only in one part of the phase space.

This kernel may lead to local inconsistencies where the acceleration term changes sign. However, when the velocity field of the carrier fluid is not constant, this inconsistency occurs only on a few grid points, where the information contained in the distribution function is negligible (see convergence results in section 4). A three point scheme would regularize the discrete acceleration discontinuity across 0, but it would represent an significant additional cost in parallel computations since it requires to communicate a larger layer at the boundaries between processes.

Time integration

In both methods, two Lagrangian steps are performed. The first one is a transport in physical space:

dx p dt = ξ ( 25
)
The integration is done with a first order Euler scheme.

The second equation is:

dξ dt = a = D(U f (x) -ξ ) (26)
Here the acceleration depends on the particle velocity. The solution would benefit from a high order integration scheme. But since the splitting is of first order and the first equation is solved with a first order Euler scheme, we use, here also, a first order Euler scheme.

One can note that in equation ( 25) ξ depends on the position of the particle through [START_REF] Koumoutsakos | High-resolution simulations of the flow around an impulsively started cylinder using vortex methods[END_REF] and on the velocity field U f (x). In this case, a higher order splitting method (for example a Strang splitting, [START_REF] Strang | On the construction and comparison of difference schemes[END_REF]) can be used for eq. [START_REF] Beale | Vortex methods. I. Convergence in three dimensions[END_REF]. In [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF], a second order remeshing method is presented where there is no additional computational cost with respect to a first order method.

For particle methods, the choice of the CFL number is crucial. We choose it such that a particle should not go over more than half a cell in order to have a sufficient condition to avoid crossing particles effect. We choose the time step such that:

∆t = CFL 2 min ∆x max(|ξ u |) , ∆y max(|ξ v |) , ∆ξ u max(D|U f u -ξ u |) , ∆ξ v max(D|U f v -ξ v |) (27)
In practice we will choose a CFL equal to 0.8 defined with the previous equation, although for the Lagrangian step, we could use a larger time step. Moreover, since U f (x) is also a function of time in a general case, we solve the Vlasov equation after the BGK (or ES-BGK) equation such that U f (x) is known at time t n+1 .

Numerical results

In the first three test cases, the particle methods are tested with a given velocity field, where the analytical solution can be computed. The domain is [-5,5]×[-5,5] in space and velocity. These test cases are used to validate the numerical methods (Lagrangian scheme and remeshing).

The last example is the passive transport of particles in a nozzle plume. The velocity field is given from the resolution of a rarefied model. In all test cases, the mass of the particles is constant.

Test 1: Zero velocity field

The domain is [-5,5]×[-5,5] in space and velocity. We consider a zero velocity field, constant in time. All the particles are initially concentrated in (0,0) and have a gaussian distribution in velocity space. The analytical solution can be easily computed in this case for a given mesh in space and velocity:

     x(t, ξ u0 ) = ξ u0 D (1 -exp(-Dt)) y(t, ξ v 0 ) = ξ v 0 D (1 -exp(-Dt)) (28)
where subscript 0 means at t=0. The larger is the drag coefficient D, the closer to the initial condition the solution is.

The steady state solution for different drag coefficients is shown in Figure 4 and5.

The solution is first computed with the Lagrangian method. At steady state, in each cell, the density of particles at the cell center is interpolated. The error with respect to the analytical solution does not depend on the space discretization. For the numerical test, the space grid is kept constant (51×51) and the velocity grid is refined from 11×11 to 401×401.

For the remeshing method, the comparison is also done on the number of particles in each cell. The initial distribution function is given by:

f p (x = 0, ξ , t = 0) = 100 2π|Ω x0 | exp(- |ξ -1| 2 2 )
where Ω x0 is the cell where all the particles are initially situated. In this case, the center of this cell is (x 0 , y 0 ) = (0, 0). The total number of particles initially in this cell (and therefore in the field) is N p0 = Ωx 0 R 2 f p (0, ξ , 0)dξ dx = 100. The analytical solution can be easily computed in this case with the error function. In a cell Ω i,j = [x i-1/2 , x i+1/2 ] × [y j-1/2 , y j+1/2 ]:

N i,j = Ωi,j R 2 f p (x, ξ , t)dξ dx (29)
The analytical solution of the problem gives the number of particles in a position x and it is determined by the initial microscopic velocity corresponding to that position. In other words, for a set of particles to be in x = (x, y) at steady state with U = 0, one needs (see [START_REF] Andries | Numerical comparison between the Boltzmann and ES-BGK models for rarefied gases[END_REF]):

ξ u (t = 0) = Dx ξ v (t = 0) = Dy Then, N i,j = D 2 |Ω x0 | Dx i+1/2 Dx i-1/2 Dy j+1/2 Dy j-1/2 f p (0, ξ , 0)dξ (30)
By introducing the error function erf in the previous equation:

N i,j = 25D 2 erf(Dx i+ 1 2 -1)-erf(Dx i-1 2 -1) erf(Dy j+ 1 2 -1)-erf(Dy j-1 2 -1) (31)
The simulations are run on different grids going from 11×11 to 151×151 in space and velocity.

The error with respect to the number of velocity grid points in each direction is shown in figure 6a in L 1 and L ∞ norm for the Lagrangian method while the error of the remeshing method is shown on figure 6b. This first test shows that both methods converge towards the analytical solution but with different behaviour. The Lagrangian method convergence is noisier (see figure 7a), in the sense that the error does not decrease regularly. This is due to the way in which the local number of particles is calculated. In each cell the number of particles is interpolated to find the value at the cell center with the Λ 1 interpolation kernel. Then, the results closely depends on how much information is present in the cell. However, the problem tends to disappear as the number of velocity grid points increases, in particular, when ∆x ∆ξ (see figure 7b). It explains why the error decreases suddenly around 51 points in each direction in velocity space. This is a very well known problem of Lagrangian method when one wants to recover an Eulerian field and can be prohibitive in our case if a fine grid is used in space.

For the remeshing method, the convergence results are more stable in the sense that the error decreases regularly. This method gives a first order convergence. The transport in physical space and velocity space is done exactly so the error of the numerical solution is only due to the remeshing process. A second order interpolation kernel is used to redistribute the particles which is equivalent, for a CFL smaller than one to a first order upwind scheme (in the sense of finite differences [START_REF] Magni | Accurate, non-oscillatory, remeshing schemes for particle methods[END_REF], [START_REF] Cottet | High order semi-Lagrangian particle methods for transport equations: numerical analysis and implementation issues[END_REF]). As expected, the convergence rates observed on figure 6b show a first order convergence. The error is higher compared to the Lagrangian scheme (about one order of magnitude) due to the particular type of advection field that is constant. Hence, in the Lagrangian method the error is mostly due to the remeshing process that is performed only once at steady state. In contrast with the Lagrangian method, the solution is not noisy even for large ∆ξ since the particles are redistributed at each time step. The error of the remeshing method is even lower for a small number of grid points which is a great advantage when the computational cost is important (like with the full system with a rarefied model). However, the stencil being compact (2 points in each directions), for large ∆ξ the particles are remeshed in priority along the grid axis and a bias along the coordinate axis may appear. One can also note that the L ∞ tends to stop decreasing on fine grids. This is due to the particular velocity field which is constant and makes the distribution function converge towards a Dirac mass in zero. In this specific case, an inconsistency is observed.

Test 2: Translational velocity field

The same test as Test 1 is performed with a constant velocity field that is not zero. In particular, we consider U f (x) = (1, 1). The solutions are compared at t = 2 for D = 1. The initial distribution is now:

f p (x = 0, ξ , t = 0) = 100 2π|Ω x0 | exp(- |ξ | 2 2 )
As expected, the same conclusions hold for this velocity field (see figure 8). The remeshing method gives a first order convergence. The convergence of the Lagrangian method shows a sudden decreasing when the order of the space discretization is the same in space and velocity. It is not the case for the remeshing method. One can also remark that the error is smaller with the Lagrangian method (about one order of magnitude). However, since the error is calculated on the Eulerian mesh, the Lagrangian solution has to be projected on this mesh. This is done with the same interpolation kernel as in the remeshing. Since this kernel is a second order kernel, we obtain a second order convergence.

As in the previous case, on very fine meshes, the order of convergence of the remeshing method decreases. In this case, the distribution function tends to a Dirac mass in the velocity space which generates an inconsistency of the remeshing method.

Test 3: Rotating velocity field

The data are the same as in test 1 except that a non zero velocity field is imposed: U f (x) = (-y, x) with D = 5. The trajectory of a single particle is computed with the Lagrangian scheme. The analytical solution is computed solving the system of four first order differential equations:

                       dx dt = ξ u dy dt = ξ v dξ u dt = -D(y + ξ u ) dξ v dt = D(x -ξ v ) (32)
This velocity field is not constant in space any more. The solution computed with the Lagrangian method follows correctly the analytical trajectory for a particle initially in (0,0) with initial velocity (0.2,0.2) (see figure 9) and also for a particle initially in (1,0) with initial velocity (0,0) (see figure 10).

We also consider a cluster of particles that are initially in the cell containing the point (1,0). For these particles, the problem is solved with the remeshing method, and compared to the analytical solution. Figures 11 shows the density computed with the remeshing method for a grid 121×121 in space and velocity. At t = 5, one can observe that the position of non zero particle density is not clearly concentrated in a point. This is due to the initial particle distribution that allows a diffusion of the particles in the neighbouring cell before their velocities converge towards the velocity field. This diffusion is also biased by the interpolation kernel used during the remeshing step and tends to spread the particles even when their microscopic velocities have converged to the velocity field.

Figure 12a shows the time evolution of the density at (x, y) = (0, 0) and for D = 1 for different grids with respect to the analytical solution. One can observe that refining the grid, the profiles converge towards the analytical one. Figure 12b shows the error with respect to the analytical solution for the remeshing method at t = 2 and ((x, y) = (0, 0) and for D = 1. These values have been chosen to simplify the computation of the analytical solution. Also in this case, a first order convergence is observed. Moreover, in this case, the convergence is smooth and does not degrade on fine meshes. Since the advection field is not constant, the phenomenon observed in the first two test cases does not occur. The distribution function does not tend to a Dirac mass either in space or in velocity. The remeshing method remains consistent. Now that the two methods have been validated on several test cases where the analytical solution was known, we focus on an actual application where the velocity field is computed from the rarefied model. The previous tests also show the drawbacks of the different methods. In particular, the necessity to have a sufficient number of particles (sufficient number of velocity grid points) for the Lagrangian method becomes prohibitive for a real test case. Moreover, in a real test case, the velocity field computed from the rarefied model is usually not constant avoiding the consistency issue of the remeshing method. Thus, in the following, we will focus on the ability to simulate the one-way coupled system solving the particle transport with the remeshing method.

A nozzle plume

The coupling between the rarefied model and the particle model is investigated through the ejection of particles in a nozzle plume. This phenomenon has been observed experimentally in [START_REF] Dettleff | Basics of Plume Impingement Analysis for Small Chemical and Cold Gas Thrusters[END_REF]. This is a challenging test case that implies several physical phenomena and since it requires the use of rarefied models it computationally intensive. Although no quantitative data are available for this kind of problem, it is worth studying it in a first step qualitatively. It is then important to quantify how many particles come and stick on the collar (represented by Γ in figure 13) and to evaluate the resulting opacity after the firing of the thruster. This evaluation requires the ability of simulating a nozzle plume in highly rarefied environment (space) and the transport of particles in this flow. This simulation is done in 2D with the geometry and initial conditions presented in figure 13. The geometry is represented by the zero-isoline of the so-called level set function defined as the signed distance between a grid point and the nozzle. It is prolonged from the nozzle outlet to the right part of the domain. As the jet expands, the moving part of the level set (see figure 13) represents the contact discontinuity between the gas coming from the nozzle and the surrounding gas initially at rest. The velocity of the moving part of the level set is the one imposed as a boundary condition for the gas. It is calculated solving a Riemann problem between the fluid state and the surrounding pressure P atm considering the Riemann invariant. This velocity is then used to transport the level set function with a WENO5 scheme. When a cell initially in the surrounding domain is reached by the jet, it is initialized with a Maxwellian distribution function computed with the corresponding boundary condition given by the Riemann problem. On the nozzle, the Euler asymptotic preserving boundary condition presented in [START_REF] Bernard | Accurate Asymptotic Preserving Boundary Conditions for Kinetic Equations on Cartesian Grids[END_REF] is used. It enforces an equilibrium distribution function on the boundary computed with a zero normal velocity and the continuity of the temperature. The density is imposed considering a zero mass flux through the boundary. Full details on the method used to describe the level set function and its transport can be found in [START_REF] Bernard | Accurate Asymptotic Preserving Boundary Conditions for Kinetic Equations on Cartesian Grids[END_REF]. As initial condition, we take T tot = 0.6, P tot = 1. The pressure at the outlet can be approximated with quasi-1D relationship supposing an isentropic flow. If we consider that M = 1 is reached at the throat of the nozzle (which is true for P atm low enough), the outlet pressure only depends on the areas of the throat and the outlet. Hence, it is a constant for a fixed geometry. Such pressure is called adaptation pressure P c . In the following, the surrounding pressure P atm is imposed through the pressure ratio r = P c /P atm . The Knudsen number Kn ∞ is set to 10 -5 corresponding to the conditions in the nozzle where the hydrodynamic regime is expected.

When a particle hits the nozzle, we assume that it sticks to the solid boundary (its microscopic velocity is set to 0). Outside the nozzle, if a particle reaches the jet boundary its velocity is set to the minimum between its own velocity and the jet boundary velocity. Thus, it cannot cross the zero-isoline of the levelset and go in the surrounding domain. On the boundary of the domain, free flow boundary conditions are imposed for both the particles and the gas except at the bottom of the domain where symmetry condition is enforced. In the following, the space domain is [-1,5]×[0,4] and is discretized with 150×100 cells. The velocity space [-10,10]×[-10,10] is discretized with 151 points in each directions. It corresponds to a number of degrees of freedom of 10 9 .

The particle flux is calculated through the surface Γ (see figure 13). It has been observed that even for large pressure ratios (until 10 6 ), the steady state does not allow particles to go through Γ. Indeed the dynamics of the gas flow has already been studied in [START_REF] Bernard | Accurate Asymptotic Preserving Boundary Conditions for Kinetic Equations on Cartesian Grids[END_REF] and in [START_REF] Bernard | Simulation of Diluted Flow Regimes in Presence of Unsteady Boundaries[END_REF] and is shown, as an example, in figure 14 for P c /P atm = 200000.

Constant velocity field

The gas flow is taken at steady state for a pressure ratio P c /P atm = 200000 (it corresponds to the velocity field in the last picture of figure 14). Particles of mass m = 2.56 × 10 -6 are injected at the inlet of the nozzle at t = 0 through the same distribution function as in test case 1.

Figure 15 shows the density field computed with the remeshing method. We mention that the lower value of the density field is 0 but for convenience, we set it to 0.0001 to see the field in logarithmic scale. The remeshing method overestimates the density near the symmetric axe of the nozzle. This is due to the remeshing process that partly remeshes the particles that arrives on the symmetry axis. This overestimation is concentrated on the first cell and vanishes as ∆x goes to zero. No particles are turning back at the outlet of the nozzle even if the angle of the jet does. This is mostly due to the inertia of the particles when they arrive at the outlet but also to the numerical discontinuity of the velocity at the outlet boundary of the nozzle (in the first cell outside the nozzle, the horizontal velocity is positive while in the jet immediately above the nozzle, it is negative). A very fine grid in space and velocity is required to observe the particle turning back, something that we cannot afford. However, using a larger stencil for the remeshing step could be a solution to make the particles go through this discontinuity. But such stencils introduce too much numerical diffusion for our needs. Another solution is to add a diffusion term in the velocity divergence that physically corresponds to the Brownian motion of the particles due to the temperature for example [START_REF] Boudet | Shock Front Width and Structure in Supersonic Granular Flows[END_REF], [START_REF] Boudin | Fluid-particles flows: a thin spray model with energy exchanges[END_REF], [START_REF] Einstein | On the movement of small particles suspended in stationary liquids required by the molecular-kinetic theory of heat[END_REF]. Here, we choose to add a perturbation in the velocity space. After the transport in velocity, a white noise is added to each microscopic velocity such that: with rand(-1, 1) a random number between -1 and 1 and ∆ξ the velocity grid spacing.

ξ n+1 = ξ n + ∆tD(U f (x) -ξ n ) + rand(-1, 1)∆ξ

Time dependant particles flow

We are now interested in a realistic test case of the ejection of particles during the firing of the thrusters. The firing of thrusters in rarefied conditions induces a pollution of optical devices (mirror, lenses) that are usually present on satellites especially during the transitional state. We want to quantify how many particles go towards the optical devices placed on a satellite. These devices are usually located above the inlet of the nozzle perpendicularly to its symmetric axis. In initial conditions, the nozzle is filled with gas and particles. At the inlet, particles are injected continuously. This boundary condition is imposed through a constant distribution function in the ghost cell equal to the initial state. We recall the initial conditions on figure 13. We specified here a surface Γ that goes from above the inlet of the nozzle and until the upper boundary of the domain. It represents the position of optical devices placed on the satellites. In the following, we compute the particle flux through Γ, to understand whether the particle flow can contaminate (or damage) the devices. The gas flow field is solved with the rarefied models while the particles are injected as the gas flow evolves. In particular, for a pressure ratio of 200000, the velocity fields at different times is shown in figure 14 for the BGK model.

Figure 16 shows that the particles go out of the domain (through the surface Γ on figure 13) sooner for higher pressure ratios. The jet reaches faster the inlet of the nozzle because the velocity of the jet boundary increases with the pressure ratio. On the same figure, one can note that even for large pressure ratios, the particles flux through Γ goes to zero at steady state. Indeed, at steady state, the expanded jet turns back but not until the left boundary of the domain. Thus, no particles are ejected in this direction at steady state. On figure 16 one can remark that there is a peak in the flux when the particles start to go out of the domain. This bump is actually higher for P c /P atm = 10000 than for P c /P atm = 10 6 . It is due to the particles that have initially a higher microscopic velocity than the boundary of the jet. Since the particles are not allowed to cross the boundary (due to the boundary condition we chose) they artificially stick to the jet contour. As the pressure ratio increases, the jet boundary velocity increases and less particles have a microscopic velocity higher than the jet boundary velocity. Thus, less particles stick to the jet boundary making the peak lower.

Comparing figures 16a and 16b one can deduce that the difference between the BGK model and the ES-BGK model increases with the pressure ratio r = P c /P atm . A more precise comparison is shown on figure 17a for two other pressure ratios (r = 10 4 , r = 4.10 5 ). The general behaviour is the same, the particles start to go out of the domain at the same time and a peak is observed. Also, we can see that the fluxes for the ES-BGK are lower. As the pressure ratio increases, the local Knudsen number in the jet increases too (order Kn ∞ r) and explains the differences between the two models.

From these data, we can recover the number of particles that left the domain through Γ. Figure 17b shows this number for different pressure ratios normalized with the number of particles initially in the nozzle.

For pressure ratios lower than 10 3 , the expanded jet does not go back until the inlet of the nozzle so no particles are ejected in front of it.

For very high pressure ratios, the number of particles going through Γ tends to stabilize because no particles are going through Γ during steady state. The ejection of particles on optical devices is a purely transitional phenomenon.

This method allowed us to divide by three the number of particle used with respect to the number of grid points in space. It is also efficient in parallel since the load balancing can be done anticipating the cost in each cell and thus, keeping a reasonable scalability.

Comparing the results given by the two models, we can note that the BGK model seems to overestimate the number of ejected particles with respect to the ES-BGK model. For very high pressure ratios, the Knudsen number becomes of

Conclusion

In this work, we proposed a numerical method to solve applied problems dealing with particle dynamics in rarefied flows. A Vlasov equation with a drag force has been used to simulate the particle flow while the gas flow has been modelled with models of the Boltzmann equation, the BGK and ES-BGK models. We successfully applied a Lagrangian scheme with a remeshing technique to the Vlasov equation. The advantage of this method with respect to finite volume or finite difference schemes is that the equation is solved only where there are particles. It is also easy to recover an Eulerian field with this method. The method has been validated on different 2D test cases with respect to analytical solutions. We also added a Brownian motion in the particle transport model through a perturbation in the velocity space. A nozzle flow has been presented where the two rarefied models have been solved to compute the particle dynamics interacting with the gas velocity field. On this test case, we successfully reproduced the contamination phenomenon observed experimentally [START_REF] Dettleff | Basics of Plume Impingement Analysis for Small Chemical and Cold Gas Thrusters[END_REF]. Moreover, based on the results, it was possible to show that the pollution is due to the transitional regime. At steady state, no particles are ejected in front of the nozzle. However, the lack of data on such phenomenon in the literature does not allow a quantitative validation. In future work, the diffusive term including the Brownian motion could be expressed as in [START_REF] Carrillo | Simulation of fluid and particles flows: Asymptotic preserving schemes for bubbling and flowing regimes[END_REF] as a divergence in velocity space and then directly integrated in the model. Moreover, as in the case of thick sprays, a feedback from the particle motion to the gas dynamic could be included in the model as a source term in the BGK or ES-BGK equation.

Figure 1 :

 1 Figure 1: Backflow phenomena and position of optical devices around the satellite thruster.

  is known. The case of the ES-BGK model is slightly more complicated. Computing the distribution function at stage k requires the Gaussian distribution G (k) fi,j . For the ES-BGK model, the moments of the Gaussian G (k)

Figure 2 :

 2 Figure 2: 1D configuration after a transport step in physical space (same behaviour in velocity space with ξ u = a u )

Figure 3 :

 3 Figure 3: Stencil with interpolation weights associated to the cells for the remeshing of a particle in x p .

Figure 4 :

 4 Figure 4: Normalized density distribution at steady state for D=0.5,1

Figure 5 :

 5 Figure 5: Normalized density distribution at steady state for D=5,10

Figure 6 :

 6 Figure 6: Test 1: Errors for the Lagrangian scheme and the remeshing method with respect to the analytical solution (D = 1).

  Grid 51×51 in space, 41×41 in velocity.

  Grid 51×51 in space, 51×51 in velocity.

Figure 7 :

 7 Figure 7: Test 1: Steady state solution by the Lagrangian method with two different velocity grids (D = 1).

Figure 8 :

 8 Figure 8: Test 2: L 1 , L 2 and L ∞ norm of the error.

Figure 9 :

 9 Figure 9: Test 3: Particle position and velocity over time solved with the Lascheme. Analytical solution in solid line. Particle initially in (0,0) with (u,v)=(0.2,0.2).

  particle

Figure 10 :

 10 Figure 10: Test 3: Particle position and velocity over time solved with the Lagrangian scheme. Analytical solution in solid line. Particle in (1,0) with (u,v)=(0,0) at t = 0.

Figure 11 :

 11 Figure 11: Test 3: Particle density with the remeshing method.

  Density solution in time at (x, y) = (0, 0).

  Convergence order of the remeshing method

Figure 12 :

 12 Figure 12: Test 3: density profile at (x, y) = (0, 0) and convergence order with the remeshing method with D = 1 at t = 2.

Figure 13 :

 13 Figure 13: Initial configuration of the computational domain for the simulation of the particle flow in a nozzle plume.

  (a) t = 0.64 (b) t = 3.22 (c) t = 6.44 (d) t = 9.8 (e) t = 14.09 (f) Steady state

Figure 14 :

 14 Figure 14: Mach number and velocity vectors for P c /P atm = 200000 with BGK model.

Figure 15 :

 15 Figure 15: Density field computed with the remeshing method, D = 1.

Figure 16 :

 16 Figure 16: Particles flux for three pressure ratios (1200, 10000, 10 6 ).

  Particle flux for the BGK model and the ES-BGK model for two different pressure ratios (10 4 , 4.10 5 ).

  Number of particles exiting the domain through Γ.

Figure 17 :

 17 Figure 17: Particles fluxes and number of particle leaving the domain for different pressure ratios and different rarefied models.

Acknowledgement

Experiments presented in this paper were carried out using the PlaFRIM experimental testbed, being developed under the Inria PlaFRIM development action with support from LABRI and IMB and other entities: Conseil Régional d'Aquitaine, FeDER, Université de Bordeaux and CNRS (see https://plafrim.bordeaux.inria.fr/).

This study has been carried out with financial support from the French State, managed by the French National Research Agency (ANR) in the frame of the "Investments for the future" Programme IdEx Bordeaux (ANR-10-IDEX-03-02), Cluster of excellence CPU.

The authors thank Dr. Lisl Weynans for fruitful discussions.