Asymptotic Preserving scheme for a kinetic model describing incompressible fluids

The kinetic theory of fluid turbulence modeling developed by Degond and Lemou in 7] is considered for further study, analysis and simulation. Starting with the Boltzmann like equation representation for turbulence modeling, a relaxation type collision term is introduced for isotropic turbulence. In order to describe some important turbulence phenomenology, the relaxation time incorporates a dependency on the turbulent microscopic energy and this makes difficult the construction of efficient numerical methods. To investigate this problem, we focus here on a multi-dimensional prototype model and first propose an appropriate change of frame that makes the numerical study simpler. Then, a numerical strategy to tackle the stiff relaxation source term is introduced in the spirit of Asymptotic Preserving Schemes. Numerical tests are performed in a one-dimensional framework on the basis of the developed strategy to confirm its efficiency.

various turbulence models. The concepts from Kinetic Theory of Gases have been of significant use in some of the models for turbulence. Extending this approach further, Degond and Lemou [7] developed a model fully based on Boltzmann type equation and its expanded analogy for a description of turbulence. Chen et al. [4] also followed a similar approach in developing a different model. This work is a first step, following the approach of [7] for study, analysis and simulation of turbulence from kinetic theory.
In [7] it is proposed to describe a turbulent incompressible fluid flow through a probability distribution function f (t, x, v) of fluid elements (or structures) depending not only on time and on the position x of these elements but also on their velocity. In this description, the fluid elements are assumed to interact in order to bring the fluid to an isotropic distribution in velocity (isotropic turbulence) and this is described by a relaxation collision operator. On the other hand, the pressure in the fluid acts in order to maintain the probability nature of the distribution function in such a way that the incompressibility (divergence-free) condition is satisfied. In the relaxation term (collision kernel), the relaxation time can be tuned in a such way that the model incorporates some turbulence phenomenology and is able to describe two important regimes: the molecular viscosity regime and the turbulence (the so-called inertial range) regime which follows the well-known Kolmogorov law on the energy spectrum of the fluid structures. In particular, this relaxation time should depend on the turbulent velocity and, as we shall see, this makes more difficult the construction of efficient numerical schemes solving the corresponding kinetic model. Let us mention that this kind of relaxation operator has been studied in different frameworks such as rarefied gas dynamics in [2], [23], [20] or to describe wave-particle collision in plasma physics [9], [10], or for cometary flows modeling [11,14].
We first give the basic kinetic equation introduced in [7]. The distribution function f (t, x, v), with v ∈ R d , x ∈ Ω ⊂ R d , d = 1, 2, 3, t > 0, satisfies the following equation where η > 0 is a constant parameter. We shall now define the different terms in this equation. First the operator Q diss is a dissipative operator which allows to take into account the dissipation of turbulence in a fluid, that is its tendency to become laminar. We shall omit it here in our preliminary numerical study and refer to [7] for a detailed description and different possible choices of such operators. Secondly, the operator Q u f is the relaxation collision operator that brings the distribution function to an isotropic state and is given by where u f is defined from the following (nonlinear) equation and the relaxation time τ is a function of the microscopic energy |v−u f | 2 2 which will be made precise later on. The operator Π is the projection on to the space of isotropic functions around u f where S d−1 is the sphere in R d (|S d−1 | denotes its measure). The distribution function f is constrained to satisfy and this allows to see the pressure P in (1) as a Lagrange multiplier of this constraint. Note that this kinetic model was also used in [7] to derive a k − ε type turbulence model or Modified k−ε turbulence model (M-k−ε). Note that the collisionless version of (1) with the constraint (3) has been studied in [3] where it is proved that the collisionless model is not always well-posed, but there exist initial conditions f (t = 0) for which the problem is well-posed Defining the mean velocityū = R d vf (v)dv, we remark that in general, we have u f =ū. However when the relaxation time τ is a (nonzero) constant function, then these two velocities coincide (see (2)). Now we emphasize that the collision operator Q u f (f ) satisfies the important properties of conservation and leads to an entropy function (see [7] for details). In particular it preserves the mass, momentum and energy (see [7]) and it is worthwhile to see that the constraint (2) is necessary to ensure the momentum conservation property.
The goal of this work is to provide numerical simulations to efficiently solve the model (1)-(2)-(3) satisfied by f, u f and P. There are two main difficulties in constructing such efficient numerical schemes for this model. First, the determination of the velocity u f from formula (2) requires the use of an iterative method like Newton-Raphson technique at each time step. This turns out to be numerically expensive and avoiding it would significantly accelerate the simulation, especially if one deals with multidimensional models. The second difficulty is with the stiffness of the collision operator. The usual way of solving stiffness problem is to use an implicit discretization for the collision operator, but this is not possible here since u f should be determined before. We note that advancing the mean velocityū is possible but advancing in time the velocity u f is not possible since it is determined by (2). This means in particular that classical approaches based on Asymptotic Preserving schemes for a class of other stiff kinetic equations (see [16,18,12,13,17,5,1,19,15,21]) do not work in this context unless some modifications are performed upon the equation as we shall see below.
The appellation Asymptotic Preserving has been introduced in [16] for numerical schemes that are stable with respect to a small parameter (η in this work) and degenerate into consistent numerical schemes for the limit model when η goes to zero. This class of numerical schemes is particularly well adapted to our framework since we want to solve (1)-(2)-(3) for arbitrary small values of η, which may lead to a severe constraint on the time step (it has to be of order η for stability) when a standard explicit numerical scheme is employed.
The above mentioned difficulties in designing an AP scheme are overcome in this framework with an efficient strategy as described in the following steps. First we perform a suitable change of frame to make possible the use of a fully implicit scheme for the collision term. Usually, this implicit treatment is at the heart of designing AP schemes. Unfortunately this is not sufficient to get such AP schemes in our context, since an additional stiff transport term arises after the change of frame (as we shall see in the next section) and therefore this requires a specific treatment. In fact, this new stiff term originates with the constraint (2) in the new frame. The numerical scheme we propose is then based on two main steps. Firstly, we focus on the stiffness of the collision operator and make it fully implicit by following for example the strategies in [18,13]. As mentioned above, the obtained scheme still contains a stiff transport term. We then propose a suitable way to consider making this new stiff term implicit: a part of it is actually made implicit. We emphasize that the whole obtained implicit scheme is computationally explicit, which means that no additional computational step is needed to solve the implicit schemes compared to the resolution of an explicit one. In particular, this results in a first order numerical scheme which enjoys the Asymptotic Preserving property, and has the same computational cost as a standard explicit solver (it does not require any linear system to invert).
The rest of the paper is organized as follows. In section 2, we give a reformulation of the continuous kinetic model under study through a suitable change of frame. The asymptotic limit η → 0 is also studied. In section 3, a suitable time discretization is performed on the new formulation, and the spatial discretization is discussed. Finally, in section 6, some numerical results are presented in a one-dimensional setting to illustrate the behavior of the new scheme in various configurations.

1.
A reformulation of the continuous kinetic model. In this section, the model is reformulated by using a suitable change of frame. In particular, the constraints (3) and (2) are reformulated in a more tractable way in view of the construction of efficient numerical schemes. The asymptotic limit η → 0 is then studied for this new formulation.
1.1. A change of frame. As underlined in the introduction above, the collision time τ allows to incorporate some important turbulence phenomenology (viscosity regime and turbulence regime following the Kolmogorov law). For instance it is shown that in the molecular viscosity regime (small turbulent energies ξ, with ξ = |v − u f | 2 /2), the relaxation time τ must satisfy τ (ξ) ∼ C/ξ whereas in the turbulent regime (large energies ξ), τ should satisfy τ (ξ) ∼ C/ √ ξ. Then a general power-law is considered for τ in [7], τ (ξ) = ξ α with ξ = |v − u f | 2 /2, − 1 2 ≥ α ≥ −1, and extended k − ε models are derived from the corresponding kinetic equations. We will not detail this study here and refer to [7] for a thorough description of these aspects.
Here we will then consider τ as a general function of the particle kinetic energy . One of the main difficulties at the numerical level is to determine u f which is the solution of the following nonlinear equation In the following we will use the notation u := u f . Hence, for a known function f , the computation of u needs typically a Newton type algorithm, which can be very expensive. Since we are working on a model where most of the parameters depend on (v − u), it makes sense to rewrite the equations after using the transformation v → (v − u) = v (see [8]). The greatest advantage of using this change of frame is that both τ and the projector Π can now be considered independent of time and space variables. In particular, the projector in the new frame will commute with the transport operator and the time derivative, a fact of great importance in the construction of efficient numerical schemes. Using this transformation, the constraints (3) and (2) can be rewritten as: where Let us now write the equation satisfied by F . The derivatives of F in terms of f can be written as Using the above relations we can rewrite (1) as follows where E = ∂ t u + ∂ x P + (∂ x u) T u and Hence, the unknowns (F, u, P, E) satisfy the model given by the equations This model is equivalent to the original model (1)-(3)-(2) satisfied by (f, u, P).
In the following, we rewrite the constraints R d F dv = 1 and R d v τ F dv = 0 in a more tractable way to construct numerical schemes.
Let us first consider the constraint R d F dv = 1. Integrating the first equation of (8) with respect to v leads to Then, taking the divergence of the first moment of the first equation of (8) and summing the obtained equation to the divergence of the second equation of (8) leads to (using (9)) (10) The left hand side gives ∂ 2 which finally gives an equation for the pressure P − ∆P = ∂ 2 x : where " : " denotes the contracted product of two tensors. We refer the reader to appendix 6.1 for the details of the computations.
Let us now consider the constraint Multiplying the first equation of (8) by v/τ and integrating with respect to v gives (using the constraint . Finally, the system to solve is the following, satisfied by (F, u) Remark 1. The case τ constant (equal to 1) leads to specific expression of E since the constraint (5) becomes R vF dv = 0 and as a consequence the stiff term in the last equation of (13) vanishes. Then, in this configuration, the equation for E reduces to

and the equation for F becomes
This equation has the same structure as the Vlasov-Poisson-BGK model in the hydrodynamic regime. Applying Π 0 to this equation enables to derive a non stiff macro equation. At the numerical level, IMEX methods proposed in [1,6,12,13,15,21] can be applied in a straightforward way. (13) can be formulated in a conservative form

Remark 2. The equation for F in
from which the conservation of mass d dt f (t, x, v)dxdv = 0 is deduced easily.
1.2. Asymptotic models from the reformulation. The aim of this section is to derive the asymptotic model, i.e., the model obtained formally from (1) when η goes to zero, or equivalently, from (13). It is of importance to understand the asymptotic limit at the continuous level when one wants to design Asymptotic Preserving schemes. To alleviate the notation, in the following we use brackets for the velocity integrals · := R d · dv. First of all, we rewrite (13) using the following notation so that the first equation of (13) reads where E is given by This term makes an additional stiff term appear so that the operator τ 2 F dv needs to be studied. This is the purpose of the following proposition.
Then, (15) becomes The function τ being positive, we conclude that if L(F )F = 0, then F = Π 0 F , which means that F is isotropic, i.e., F depends on |v|.
We can now derive formally the asymptotic model from (13) satisfied by (G 0 ,ū) which are the limits when η goes to zero of (F, u). This is presented in the following proposition.
Proposition 2. The asymptotic model obtained from (13) in the limit η goes to zero, satisfied by the limit (G 0 ,ū) of (F, u), writes Proof. We first derive the asymptotic model from the first equation of (13). To do that, we consider the Chapman-Enskog expansion of F : where G 0 = Π 0 F . Inserting this expansion in the first equation of (13) leads to Applying Π 0 to this latter equation leads to the asymptotic equation for G 0 The Rewritting the isotropic function G 0 (v) := g(ξ) with ξ = |v| 2 /2, we have ∂ v G 0 = v∂ ξ g so that the last term becomes We then conclude (∂ When η goes to zero, we have u converging towardsū since vF dv =ū − u with F which tends towards G 0 and vG 0 dv = 0; hence, the incompressibility condition ∂ x ·ū = 0 makes this term null when η goes to zero. The asymptotic model of the first equation of (13) then writes which is the first equation of the asymptotic model (16).
To get the limit of the second equation of (13), we need to derive the limit of E and ∂ x P when η goes to zero. The so-obtained model will be satisfied by the limit u of u.
Let us first focus on the derivation of the limit of E. To do that, we apply (I −Π 0 ) to (17) and get (neglecting O(η) terms) We multiply by v and integrate with respect to v to get The expression of the limit E 0 of E when η goes to zero then becomes We then focus on the limit equation satisfied by the pressure P 0 , limit of P when η goes to zero. Considering the third equation of (13) when F is replaced by its limit G 0 and u by its limitū, we have The asymptotic of the second equation of (13) satisfied byū which is the limit of u when η goes to zero, finally writes where E 0 and P 0 satisfy (19) and (20). This corresponds to the second equation of (16), and with (19) and (20) correspond to the third and fourth equation of (16). This concludes the proof.
2. Numerical schemes. In this section, we propose some numerical schemes for the kinetic models described previously. We divide our approach into the cases τ constant and τ non constant, the latter case being more complicated and needing different techniques. Our main goal is to provide semi-discretized (in time) numerical schemes which enjoy the Asymptotic Preserving (AP) property, that is: (i) they are stable with respect to the parameter η, (ii) they degenerate when η goes to zero towards numerical schemes which are consistent with the asymptotic model (16). Such schemes have been initiated in [16] and have become, during the last decade, very popular to approximate stiff kinetic equations.
One important problem here to derive an AP scheme is to make the relaxation operator implicit since the determination of u requires the knowledge of f . This involves a nonlinear system to solve which would be computationally very costly and difficult to handle. As presented in the previous section, one issue is to consider the change of variable v → v −u and the equation on F (7). However, an additional stiff term arises after this change of frame and the derivation of a numerical scheme where only the relaxation term 1 τ (Π 0 F − F ) is considered implicit does not provide the AP property of this scheme, as we shall see. We will also look at an exponential integrator following [18] but in this case also, this is not sufficient to get the AP property. Finally, we propose a new strategy which is a modification of the latter approaches, that ensures the numerical scheme to degenerate as η goes to zero towards a consistent numerical scheme of the asymptotic model (16). We emphasize that this numerical scheme has the same computational cost as a standard explicit solver.
In the following, we present two strategies from the literature (see [12,13,19,21]) and explain why they cannot be applied in our framework. We then propose our new scheme in two different versions. As we deal with a semi-discretization in time in this section, we introduce a time discretization t n = n∆t, n ∈ N, ∆t > 0 and we denote by F n an approximation of F at time t n .
Let us recall (13) where A is given by (14) and We introduce the following notations for the time discretization of the term E and where E n 1 is given by State of the art. In this part, we look at two standard strategies to derive an AP numerical scheme for (21) where A is given by (14).
As we shall see, the presence of the stiff term in the left hand side is the major obstacle in a direct application of standard strategies.
The first idea would be to consider implicit the relaxation term To get an expression of Π 0 F n+1 , we apply Π 0 to this last equation This macro equation still involves a stiff term in the term E n,n given by (22) which requires the time step ∆t to be lower than η. This is due to the fact that the other stiff term has not been considered implicit. A second idea to design an AP scheme is to use an exponential integrator ( [18,13]). From (25), one has Integrating between t n and t n+1 leads to The key point of this strategy is the approximation of the integrals (see [18]). For the two first integrals, we choose the left rectangle quadrature whereas the last integral is computed exactly after considering Π 0 F implicit. This leads to We can remark that when η goes to zero, F n+1 goes to Π 0 F n+1 , which was not the case in the previous approach. One needs now to compute Π 0 F n+1 . To do that, we apply Π 0 to this last equation to get (recalling Π 0 commutes with isotropic functions) where E n,n is given by (22). Here again, the macro equation (26) still contains a stiff term which cannot be considered implicit without requiring the inversion of a nonlocal operator. Let us remark that other choices of time integrals approximation do not lead to a semi-implicit Asymptotic Preserving scheme as well.
A slightly different strategy is then required since standard ones do not work. In the following, two schemes are proposed which enjoy the AP property. We point out that even if these new schemes involve an additional implicit term, namely E n,n+1 instead of E n,n , they have the same computational cost as an explicit numerical scheme.
New strategy: 1 (EXPO). In the sequel, we present a modification of the previous exponential scheme. The numerical scheme for F writes The main problem now is to determine (v/τ 2 )F n+1 in E n,n+1 given by (23) without requiring the inversion of a linear system. Integrating this latter equation with respect to v, after multiplying by (v/τ 2 ), leads to since (v/τ 2 )(1 − e −∆t/(ητ ) )Π 0 F n+1 = 0. Let us focus on the last term to deduce (v/τ 2 )F n+1 in an explicit way where E n 1 is given by (24). We can inject this term in (28) to get an expression of (v/τ 2 )F n+1 Then, we have determined in an explicit way (v/τ 2 )F n+1 and we can remark that it is of order η. Then, the term E n,n+1 can be computed with an explicit complexity thanks to where E n 1 is given by (24). The last thing to do is to express Π 0 F n+1 . To do that, we apply Π 0 to (27) (recalling that Π 0 commutes with isotropic functions) This macro equation is the same as (26) except that the term (v/τ 2 )F n+1 in E n,n+1 is now implicit. Thanks to (29), this term is of order η making the last term of (30) non-stiff.
To advance the equation for u, the stiff term contained in E also has to be considered implicit. This can be done easily since v/τ 2 F n+1 has been computed previously and E n,n+1 is known. Hence, the numerical scheme for the equation on u writes u n+1 = u n − ∆t(u n · ∂ x )u n + ∆tE n,n+1 − ∂ x P n , where E n,n+1 is computed from (23) and P n is computed from From (F n , u n ), the algorithm to compute (F n+1 , u n+1 can be summarized as follows Let us emphasize that the computational complexity of this algorithm is the same as the complexity of an explicit algorithm. The second step (computation of (v/τ 2 )F n+1 with (29)  • is consistent and stable for (7); • degenerates as η goes to zero towards a numerical scheme which is consistent with (16).
Proof. We observe that formally, (27) degenerates when η goes to zero to F n+1 = Π 0 F n+1 . Thanks to (29), the term E is of order one since the stiffness has been stabilized by considering E n,n+1 for its approximation. Hence, applying Π 0 to (27) and considering the limit η goes to zero, the term E n,n+1 · Π 0 (∂ v F n ) goes to zero since F n converges towards Π 0 F n and Π 0 (∂ v Π 0 F n ) = 0. Finally (30) degenerates as η goes to zero to with AF n given by (14). The same computations as in the continuous case leads to which is an explicit discretization of the asymptotic model for G 0 := Π 0 F (first equation of (16)). Applying now (I − Π 0 ) to (27), and since F n+1 = Π 0 F n+1 + O(η), one obtains Multiplying by ve −∆t/(ητ ) /∆t and integrating with respect to v (following the computations of the continuous case) leads to

CROUSEILLES, LEMOU, RAO, RUHI AND SEKHAR
or, after an integration by parts (using Π 0 F n = 1) which is the same expression as the limit equation satisfied by E (last equation of (16)). The equation for the pressure P does not present any difficulty since (32) degenerates into (u n converges towardsū n as η goes to zero) Finally, the equation (31) for u becomes, with the previous notations which degenerates when η goes to zero intō which is an explicit time discretization of the asymptotic equation onū (second equation of (16)). Finally, we check that this numerical scheme enjoys the Asymptotic Preserving property.
New strategy: 2 (DIMP). Here, a slight modification of the previous scheme is presented. Instead of using exponential integrators, we use here a simple implicit scheme, including the same treatment as before of the stiff transport term. The numerical scheme for F writes where A is given by (14) and E n,n+1 is given by (23), so that, with δ = 1/(1 + ∆t/(τ η)), one obtains As previously, one needs to determine (v/τ 2 )F n+1 ; this is done as in the previous case by integrating this latter equation against (v/τ 2 ) (34) Then, the term E n,n+1 can be computed thanks to where E n 1 is given by (24). The last thing to do is to express Π 0 F n+1 . To to that, we apply Π 0 to (33) and get an expression of Π 0 F n+1 The rest of the numerical scheme (computation of u n+1 , E n 1 and P n ) is the same as in the previous version. Hence, from (F n , u n ), the algorithm to compute (F n+1 , u n+1 can be written as follows • compute P n with (32), • compute (v/τ 2 )F n+1 with (34), • compute E n,n+1 with (35), • compute Π 0 F n+1 with (36), • compute F n+1 with (33), • compute u n+1 with (31). As a conclusion, we derive the Asymptotic Preserving scheme (36)-(34)-(31)-(33)-(32)-(35) for (7), the properties of which are summarized in the following proposition. • is consistent and stable for (7); • degenerates as η goes to zero towards a numerical scheme which is consistent with (16).
Proof. The proof is similar to the previous one and is left to the reader.
3. Spatial discretization. In this section, we focus on the one-dimensional case and we briefly detail the phase space discretization used to approximate the differential operators in (27) or in (33). Let us write the "EXPO" numerical scheme in the one-dimensional configuration (the DIMP version is similar) where AF n is given by and E n,n+1 is given by with E n 1 given by If moreover periodic boundary conditions are used, one can show thatū is a constant. Indeed, considering the first moment of (1) enables to write an equation for and using the change of frame, we get The incompressibility condition ∂ xū = 0 leads to the equation satisfied by the pressure P (applying ∂ x to (40)) −∂ 2 x P = ∂ 2 x (v + u) 2 F . Periodic boundary conditions enable us to write −∂ x P = ∂ x (v + u) 2 F and we deduce from (40) ∂ tū = 0. In addition to the incompressibility condition ∂ xū = 0, we conclude thatū is constant. As a consequence, since vF =ū − u, u is then determined directly from this latter relation and the equation (31) is not needed in this simplified one-dimensional framework.
In the rest of this section, we describe the phase space discretization. We introduce a phase space uniform grid of size ∆x > 0 in the x-direction and ∆v > 0 in the v-direction: x i = x min +i∆x, i = 0, . . . , N x and v j = v min +j∆v, j = 0, . . . , N v , with N x , N v ∈ N and x min , v min ∈ R. In the sequel, we denote by F n i,j an approximation of F (t n , x i , v j ).
We consider upwind (forward/backward difference depending on wave-speed) schemes. Typically, to approximate AF = (v + u)∂ x F − (∂ x u)v∂ v F (see (14)) and E n 1 given by (39), we use the following approximations with the standard notation a ± = a±|a| 2 , where D x,v ± denotes a first order one-sided difference operator (for example (D x . All the integrals with respect to v are approximated by a rectangle quadrature.
More sophisticated spatial discretizations can be used, such as WENO schemes to achieve high orders, but our goal here is to validate the whole strategy.
An important point is to ensure the conservations of R F dv and R (v/τ )F dv at the discrete level. The previous spatial discretizations do not ensure these conservations. To overcome this problem, the scheme should be extended using some correction so that the above conservations are ensured irrespective of the discretization used. Let us mention that the above conservations can be ensured without an a posteriori correction by computing the discrete version of E and P by following the same strategy as in the continuous level. However, this is more restrictive and we choose the first approach which needs an a posteriori correction which we detail below.
To do that, we propose to a posteriori modify the obtained F n+1 i,j so that conservations are ensured. Say we use discretization B to the scheme to getF n+1 = BF n following the time and spatial discretizations described above. In general, even if F n satisfies the conservations j F n i,j ∆v = 1 and j (v j /τ j )F n i,j ∆v = 0, the obtainedF n+1 does not preserve the properties which the corresponding equation in continuous level has. Then, we introduce an orthogonal projector P onto span{(1, (v/τ )) F n } in L 2 (dv/F n ), such that We hence modifyF n+1 as follows An example of such a projector can be given as from which we easily deduce P(g) = g and v τ P(g) = 0. At the discrete level, the integrals are replaced by the sums on the index j.
Introducing the defined projector P along with the scheme gives the benefit of using any higher order upwind discretization for space and velocity in the scheme. This way of correction does not depend on the scheme used, i.e., the scheme may be implicit or explicit.
4. Numerical results. This section is devoted to the validation of the numerical scheme presented above. Different schemes will be compared here in the onedimensional configuration (d = 1) noting that the strategy presented in section 2 can be applied to arbitrary dimensions. First, the scheme (30)-(29)-(27) will be called "EXPO" and compared with explicit schemes such as Runge-Kutta schemes for (7). Let us remark that the "DIMP" scheme gives results which are very similar to those obtained by "EXPO" so that only results obtained by EXPO will be shown here. For these latter schemes, the time step ∆t will be chosen to satisfy ∆t ≤ η which is not the case for "EXPO" scheme.
The numerical tests will be of increasing difficulty: firstly, the homogeneous case with non constant τ will be tested and second, the non homogeneous case with a constant and non constant τ will be also validated.
The first constraint F = 1 is automatically preserved by this model. The second constraint (v/τ )F = 0 should be reformulated. Taking (v/τ ) moment leads to The model to solve in this context is then originating from (1) in the homogeneous case. We detail the numerical scheme in this simplified configuration in Appendix 6.3.
The following initial condition is chosen F (t = 0, v)) = F 0 (v) such that In Figure 1, the time history of E(t) = F (t) − Π 0 F (t) L 2 is plotted for different numerical schemes (explicit RK schemes and EXPO) for η = 1. The same conclusions as before follow since EXPO and RK's are nearly superimposed. Note that the rate of relaxation is quite different from the previous case. Figure 2 displays the distribution function F (t, v) for different times t = 0, 1, 2, 5, 100, 1000 using EXPO scheme for η = 1. This illustrates the relaxation towards the isotropic state. The rate of isotropisation is not the same for all the velocities since we are dealing with a non constant relaxation time τ (see also [20]). In Figure 3, E is plotted for different values of η using the EXPO scheme and RK3. For this latter scheme, the time step is fixed as ∆t = 0.01 for stability, whereas ∆t = 0.1 for EXPO. We observe for η = 0.5 that EXPO is very close to RK3, even if the time step is 10 times smaller for RK3. For η = 0.1, the two curves are not so close (Figure 4, right) and we also display the result obtained by EXPO with ∆t = 0.01; we then observe that the two curves are nearly the same, which confirms the convergence in time of EXPO. In Figure 4, results obtained with EXPO (∆t = 0.1) are shown, where the value of η is chosen equal to 10 −2 and 10 −4 ; we observe that at the first iteration, the solution is very close to its projection (since E is nearly zero), even if the initial condition is not at equilibrium. This emphasizes the AP character of our numerical scheme.
The numerical parameters are the following: ∆t = 0.2∆x/v max , N x = 64, N v = 256 and v max = 16.
In Figure 5, the time history of E(t) = (F − Π 0 F )(t) L 2 (the L 2 norm is considered in x and v) is shown, for different values of η. The new scheme "EXPO" is compared to RK3 and RK1. The new scheme turns out to be competitive in the kinetic regime since the results are very close to those obtained with RK3 or RK1. When η is smaller (η < 10 −3 ), the considered ∆t = 0.2∆x/v max is larger than η so that the explicit RK methods become unstable and provide wrong solutions. In this     case, a resolved time step is needed for these explicit methods whereas the method "EXPO" still provides the right solution with the time step ∆t = 0.2∆x/v max (see Figure 5 where η = 10 −5 ). Only results obtained with EXPO scheme is shown. A very strong relaxation can be observed, confirming the fact that the solution is projected onto the set of isotropic functions at the first iteration.
In Figure 6, the time history of E(t) = (F − Π 0 F )(t) L 2 (the L 2 norm is considered in x an v), for different values of η and different numerical schemes (Runge-Kutta 1 and 3, EXPO) is shown. For η = 10 −2 , comparison with explicit RK method needs ∆t = 0.05∆x/v max whereas the "EXPO" method still uses ∆t = 0.2∆x/v max . We observe that the three methods have a similar behavior which validates the EXPO method for η = 1, 10 −1 , 10 −2 . When η is smaller, explicit methods are too costly so that only results obtained by the EXPO method are presented. In particular, in Figure 7 when η = 10 −4 , one can see that the initial condition is immediately projected onto isotropic functions.

5.
Conclusion. In this work, a new scheme is proposed for the numerical simulation of a kinetic theory based model of turbulence. The derivation of this scheme needed two important tools: (i) the change of frame enables to avoid the delicate issue of the computation of u; (ii) the improvement of standard AP schemes since the change of frame makes an additional stiff transport term appear. Even if this scheme is implicit to avoid the severe constraint coming from the stiff relaxation source term and the additional stiff transport term, its numerical cost is the same as an explicit scheme; moreover, it enjoys the Asymptotic Preserving property, i.e., it degenerates when the stiff parameter goes to zero, to a numerical scheme consistent with the asymptotic model. The current work presents first results obtained with a prototype one-dimensional model. More general configurations can be tackled in future, in order to study turbulence from kinetic theory in a deeper way. In addition, the spatial discretization can be improved to capture physical turbulence phenomena; as mentioned above, high order WENO schemes can be quite easily adapted to our time discretization strategy.
6.1. Computation of pressure. The equation for the computation of pressure P is presented here in detail. The starting equation is (10). Forgetting for a while the ∂ x ·, the left hand side of (10) becomes Indeed, using the relations Hence, from (10), we deduce by taking the divergence of (42) the equation for P −∆P = ∂ 2 x : (v + u) ⊗ (v + u)F dv .