L Bois

E Franck

L Navoret

V Vigon

A neural network closure for the Euler-Poisson system based on kinetic simulations

This work deals with the modeling of plasmas, which are chargedparticle fluids. Thanks to machine leaning, we construct a closure for the one-dimensional Euler-Poisson system valid for a wide range of collision regimes. This closure, based on a fully convolutional neural network called V-net, takes as input the whole spatial density, mean velocity and temperature and predicts as output the whole heat flux. It is learned from data coming from kinetic simulations of the Vlasov-Poisson equations. Data generation and preprocessings are designed to ensure an almost uniform accuracy over the chosen range of Knudsen numbers (which parametrize collision regimes). Finally, several numerical tests are carried out to assess validity and flexibility of the whole pipeline.

Introduction

Plasmas are gases composed of charged particles, i.e. ions and electrons, which are the subject of many studies because they are a very common state of matter in the universe (e.g. in stars, ionosphere). They are also present in industrial devices, like in fusion reactors. Plasma dynamics is complex since particles have both long range interactions, through the self-consistent electromagnetic fields, and short-range collisions. The most complete model to describe plasmas is given by the Vlasov equation, which is a kinetic equation satisfied by the distribution function of particles in the position-velocity (x, v)-phase space and which is coupled with the electromagnetic equations. To reduce the dimention, it is possible to use fluid models which are derived from the kinetic Vlasov equation with a collision operator by taking the velocity moments corresponding to the monomials 1, v, v 2 . This raises up three equations linking the four following quantities: the density (the zeroth order moment), the mean velocity (involving the two first moments), the temperature (involving the three first order moments) and the heat flow (involving the four first moments). Fluid models are very attractive as they are much computationally cheaper as they involve only spatial quantities. However, they require an additional equation called a closure to link the heat flow to the other three moments.

The first two fluid models obtained in the collisional regime are the Euler and the Navier-Stokes system. The closure is obtained from an asymptotic Chapman-Enskog procedure, assuming that the distribution function is close to be at thermodynamical equilibrium due to collisions. These equations are thus valid for small Knudsen number ε > 0, which is the mean free path between two collisions divided by the characteristic length studied. The Navier-Stokes system provides O(ε) corrections of the pressure tensor and heat flux. Higher order corrections have been considered but lead to ill-posed systems (Burnett equations). We refer to [START_REF] Braginskii | Transport phenomena in plasma[END_REF] for a specific derivation in the case of plasmas and to [START_REF] Degond | Macroscopic limits of the boltzmann equation: a review[END_REF] for a review.

To extend the validity of fluid models for larger Knudsen numbers, larger moment systems have been proposed for the Boltzmann kinetic equation, with neutral particles: the Grad 13 order moment model using perturbative theory [START_REF] Grad | On the kinetic theory of rarefied gases[END_REF], the Levermore 14 order moment model using entropy maximisation closure [START_REF] Levermore | Moment closure hierarchies for kinetic theories[END_REF][START_REF] Levermore | The gaussian moment closure for gas dynamics[END_REF]. They both suffer from intrinsic mathematical imperfections: a lack of hyperbolicity for the former [START_REF] Cai | Globally hyperbolic regularization of grad's moment system[END_REF], so-called realizibility issues for the latter [START_REF] Junk | Domain of definition of levermore's five-moment system[END_REF][START_REF] Schneider | Entropic approximation in kinetic theory[END_REF]. However, recent developments propose corrections for the Grad model [START_REF] Cai | Globally hyperbolic regularization of grad's moment system[END_REF] and the entropy methodology and has been very successful in the context of radiative transfer with compact velocity domain [START_REF] Dubroca | Etude théorique et numérique d'une hiérarchie de modèles aux moments pour le transfert radiatif[END_REF][START_REF] Garrett | A comparison of moment closures for linear kinetic transport equations: the line source benchmark[END_REF].

Specific closures have been developed for plasma fluid models. In electrostatic regime, the Hammett-Perkins closure [START_REF] Hammett | Fluid moment models for landau damping with application to the ion-temperature-gradient instability[END_REF][START_REF] Hammett | Fluid models of phase mixing, landau damping, and nonlinear gyrokinetic dynamics[END_REF] is designed to recover the Landau damping effect, which results from the transfer from spatial modes to velocity modes and which leads to the damping of the electrostatic energy. To this end, the heat flux depends on the temperature through a non-local integral dissipative operator. Dissipation is also introduced in the momentum equation. Let us mention the recent work [START_REF] Manfredi | Density functional theory for collisionless plasmasequivalence of fluid and kinetic approaches[END_REF] for a discussion on such closure. When considering magnetized plasma, several closures have also been proposed for fluid models, like for Magnetohydrodynamics (MHD) or gyrofluid equations. See for instance [START_REF] Chew | Application of dispersion relations to low-energy meson-nucleon scattering[END_REF][START_REF] Snyder | Landau fluid models of collisionless magnetohydrodynamics[END_REF][START_REF] Brizard | Nonlinear gyrofluid description of turbulent magnetized plasmas[END_REF][START_REF] Negulescu | Closure of the strongly magnetized electron fluid equations in the adiabatic regime[END_REF]. Recently, fluid closures that preserve the Hamiltonian structure of the Vlasov equation have also been developed [START_REF] Perin | Hamiltonian closures for fluid models with four moments by dimensional analysis[END_REF][START_REF] Tassi | Hamiltonian closures in fluid models for plasmas[END_REF].

Other strategies consist in adding kinetic effects. Let us mention model reduction technic to obtain fluid models like the water-bag method [START_REF] Besse | The multi-water-bag equations for collisionless kinetic modeling[END_REF][START_REF] Besse | Gyro-water-bag approach in nonlinear gyrokinetic turbulence[END_REF][START_REF] Desjardins | A quadrature-based moment method for dilute fluid-particle flows[END_REF], where the kinetic distribution function is approximated with piece-wise constant functions in velocity. Kinetic effects can also be numerically introduced through micro-macro decomposition methods [START_REF] Crouseilles | A hybrid kinetic/fluid model for solving the gas dynamics boltzmann-bgk equation[END_REF][START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF][START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for vlasov-poisson-bgk equation using particles[END_REF].

A more recent approach consists in using supervised machine learning to find a data-driven closure. Artificial neural networks (ANN) have proven very efficient to interpolate data and detect underlying structures. In particular, convolutional neural networks are very efficient to analyse images and thus the outputs of numerical simulations. There are numerous works for designing physics-based models using neural networks (see for instance [START_REF] Beck | Deep neural networks for data-driven les closure models[END_REF][START_REF] Wang | Physics-informed machine learning approach for reconstructing reynolds stress modeling discrepancies based on dns data[END_REF][START_REF] Duraisamy | Turbulence modeling in the age of data[END_REF]). Such tools have already been used in two works concerning moment closure. In [START_REF] Han | Uniformly accurate machine learningbased hydrodynamic models for kinetic equations[END_REF], in the context of neutral gases, the authors propose to learn the appropriate higher moments required in the model. In [START_REF] Maulik | Neural network representability of fully ionized plasma fluid model closures[END_REF], which concerns plasmas models, the closure is directly learned from the analytic closure like the Hamett-Perkins one.

Description of our work. Here, we introduce a data-driven closure, based on a fully-convolutionnal neural network, which is valid for a large spectrum of collisional regimes (i.e. a large range of ε). Data are obtained from numerical simulations of the Vlasov-Poisson system satisfied by the distribution function and the electric potential. This method has already been mentioned in [START_REF] Han | Uniformly accurate machine learningbased hydrodynamic models for kinetic equations[END_REF] for neutral gases, but the authors preferred to develop another approach. Here we further investigate the design of the closure to obtain accurate predictions.

Collisions between charged particles would normally be modeled with the Fokker-Planck-Landau operator. However, as it is usually done, we replace it by a BGK (Bhatnagar-Gross-Krook) operator, that models the relaxation of the distribution function towards the equilibrium distributions, the Maxwellians. The inverse of the Knudsen number ε is in prefactor of this operator.

Knudsen numbers ε range is chosen equal to [0.01, 1]. Indeed, for ε < 0.01, the Navier-Stokes closure already gives good results. An upper bound of the Knudsen interval has to be prescribed. Here we choose ε = 1 for embrassing mildly collisional regimes, but larger values could be considered.

The closure takes as input three one-dimensional vectors corresponding to the spatial discretized density, mean velocity and temperature, and one vector corresponding to the parameter ε. The closure returns an estimation of the spatial discretized heat flux.

The core of the closure is a fully convolutional neural network which has a V-Net architecture. It was first described and developed by Milletari et al. [START_REF] Milletari | V-net: Fully convolutional neural networks for volumetric medical image segmentation[END_REF] to treat medical images (3D signals). It is an evolution of the U-Net, developed for biomedical 2D images [START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF]. It consists in a succession of several convolution kernels, down-samplings and up-samplings that perform multi-scale analysis of the signal. No fully connected layer is used, so the whole network has relatively few parameters. The architecture of the network can mostly be described by three hyperparameters: the number of levels of the "V", the depth (which rules the number of channels at each level), and the size of the kernels of convolution.

The influence of these parameters on the performance is investigated. As usual in neural networks construction, some processing of the data is required to predict the heat flux in the widest possible range. Therefore multiple steps of processing are included in the closure, on top of the central neural network. One of these is the normalization of the heat flux with the Navier-Stokes approximation, that helps lowering the otherwise large relative error on predictions of low heat fluxes. Another weakens the dependency of the closure on the underlying mesh size, by performing a data slicing before applying the neural network part. We also use a resampling strategy for discretizations with a resolution different from the one used to train the neural network.

The insertion of this closure in a fluid solver raises mathematical issues. Indeed, it is not guaranteed that the neural network closure is dissipative and thus it could lead to instabilities. To prevent such scenario, a smoothing of the prediction is added into the closure.

The paper is organised as follows. In Section 2, we introduce the Vlasov-Poisson system and the moment closure issue. In particular, we give the Euler and Navier-Stokes closures. Then the neural network closure strategy is explained in Section 3: prediction strategy, architecture of the network and training methodology are presented. Then Section 4 details the data generation and Section 5 their processing. Finally, in Section 6, we carry out several numerical tests to quantify and analyze the accuracy of the closure.

Fluid closure for Vlasov-Poisson

In this section we present the kinetic description of the plasma dynamics in one space dimension and its fluid approximations.

Kinetic model

A kinetic model is a model interested in the function f : (x, v, t) → f (x, v, t)
that describes the evolution of the distribution of the particles in the (x, v)phase space, where x ∈ [0, L] denotes the space variable, v ∈ R the velocity variable and t ∈ R the time. We will consider periodic boundary conditions in space. From this distribution function can be computed some macroscopic physical quantities. The first three moments give the particle density ρ(x, t), the mean velocity u(x, t) and the total energy w(x, t) defined as:

ρ(x, t) = R f (x, v, t)dv, ρ(x, t)u(x, t) = R f (x, v, t)vdv, ( 1
)
w(x, t) = 1 2 R f (x, v, t)v 2 dv. ( 2
)
We can also define the pressure p(x, t), the temperature T (x, t) and the heat flux q(x, t):

p(x, t) = R f (x, v, t)(v -u(x, t)) 2 dv, ρ(x, t)T (x, t) = p(x, t), (3) q(x, t) = R 1 2 f (x, v, t)(v -u(x, t)) 3 dv. ( 4
)
Note that we have the following relation:

w = ρu 2 /2 + ρT /2 = ρu 2 /2 + p/2.
The evolution of the distribution is described by the Vlasov equation:

∂ t f + v∂ x f -E∂ v f = Q(f ), (5)
where E(x, t) is the self-induced electric field, which satisfies the Poisson equation:

E = -∂ x φ , ∂ xx φ = ρ - [0,L] ρ dx. ( 6
)
Here φ(x, t) denotes the electric potential.

The source term Q is called a collision operator and allows the model to take into account the collisions between particles. Different collision operators can be considered to deal with different situations. In this work we use the BGK operator (Bhatnagar, Gross and Krook), built to conserve the mass, momentum and kinetic energy of the system, and model the relaxation induced by the collisions toward a local equilibrium distribution M (f ). This operator simply reads

Q(f ) = 1 ε (M (f ) -f ),
where M (f ) is called the Maxwellian of f and is given by

M (f )(x, v, t) = ρ(x, t) 2πT (x, t) e -(v-u(x,t)) 2 2T (x,t) , ( 7
)
with ρ, u and T the density, mean velocity and temperature associated to f and defined in ( 1)-( 3). The parameter ε > 0 is called the Knudsen number and represents the mean free path between two collisions divided by the characteristic length studied. In the limit ε → 0, the distribution function is expected to be closed to local equilibria. In this regime, the phase-space dynamics could be infered from the spatial dynamics of its macroscopic moments ρ, u, T . The kinetic equation ( 5) can be easily solved numerically in this one dimensional case. Several numerical methods have been proposed in the literature [START_REF] Sonnendrücker | Numerical Methods for the Vlasov-Maxwell equations[END_REF][START_REF] Dimarco | Numerical methods for kinetic equations[END_REF]. In this work, we use a Finite Difference/Finite Volume method similar to the one introduced in [START_REF] Pham | Space-only hyperbolic approximation of the vlasov equation[END_REF][START_REF] Helluy | Reduced vlasov-maxwell simulations[END_REF] and described in appendix A. Although it is possible to extend such computations in dimension 2 or 3, the computational cost becomes very prohibitive. We are therefore led to consider fluid models.

Fluid model

A fluid model of a plasma describes the evolution of its density ρ(x, t), mean velocity u(x, t) and kinetic energy w(x, t), instead of the distribution of its particles f (x, v, t) in a kinetic model. Such models are far less expensive to simulate numerically than the original kinetic model. Indeed, this requires the computation of only three spatial quantities (3 × N x unknowns), instead of the full kinetic distribution (N x × N v unknowns), where N x and N v stands for the number of discretization points in space and velocity.

Fluid models can be obtained by using the moment method [START_REF] Degond | Macroscopic limits of the boltzmann equation: a review[END_REF]. Formally, it consists in computing the first three moments in velocity of the Vlasov equation, i.e. multiplying Equation ( 5) by 1, v and v 2 /2 and then integrating in velocity:

for p = 0, 1, 2, R v p (∂ t f + v∂ x f -E∂ v f ) dv = R v p Q(f ) dv.
Since the collision operator conserves mass, momentum and energy, the right-hand sides vanish. Therefore, it results in the following system:

     ∂ t ρ + ∂ x (ρu) = 0, ∂ t (ρu) + ∂ x (ρu 2 + p) = -Eρ, ∂ t w + ∂ x (wu + pu + q) = -Eρu, ( 8
)
where p is the pressure and q the heat flux defined in (3)- [START_REF] Braginskii | Transport phenomena in plasma[END_REF]. Note that the pressure p is actually a function of ρ, u, w since we have the following relation:

p = 2w -ρu 2 .
The electric field E is still given by the Poisson equation [START_REF] Cai | Globally hyperbolic regularization of grad's moment system[END_REF].

We thus get a system of three equations on the four variables ρ, u, w and q. For it to be closed, we need a fourth equation connecting these four unknowns, called a closure. Usually this closure consists in replacing the true heat flux q with a simplified one q given as a function of the other quantities , and that can be written q = C(ε, ρ, u, T ).

Note that this closure depend on the physical regime we consider through the parameter ε. Let us also mention here that in higher dimension, the pressure stress tensor is not completely determined by ρ, u and w and an additionnal closure relation is necessary.

In fluid regimes, where most of the kinetic effects can be neglected, two closures are classically considered: the Euler or the Navier-Stokes closures. The Euler closure is valid in regimes where the distribution function is closed to be Maxwellian: f = M (f ) + O(ε). Using this ansatz, we obtain the Euler closure: q = 0.

When we are interested in regime with O(ε) deviation from the Maxwellians, f = M (f ) + εg + O(ε 2 ), we get the Navier-Stokes closure:

q = - 3 2 ε p ∂ x T. ( 9
)
Note that it is a O(ε) correction of the Euler closure. We refer to [START_REF] Degond | Macroscopic limits of the boltzmann equation: a review[END_REF] for more details. More complex closures have been developped to capture more kinetic effects, in regimes where the disbribution function is more distant from the Maxwellians. Such closures are chosen non local, meaning that the heat flux q(x) at location x does not rely on the values of ρ, u, T and their derivatives at location x only, but on their values on all the domain. Indeed, the true heat flux, given by a third order moment of the kinetic distribution (Eq. ( 4)), depends on the full distribution in velocity at location x and thus retains non-local information about the dynamics.

The goal of our work is to provide a method to provide a new closure, where the function C is implemented using a neural network. The next section describes this approach in more details.

Closure with a neural network

In this section we introduce the overall principle of our method to build the fluid closure presented in Section 2.2 with a neural network. First we describe the basic functioning of the neural network, before introducing the other operations included in the closure, and that surround the neural network. We end this section with the detailed description of our neural network and of its training.

Interpolation of the heat flux with a neural network

In this work we are interested in the ability of machine learning to find a function C that maps the Knudsen number, the density, the mean velocity and the temperature of the plasma to its heat flux. Such a mapping cannot be exact as the heat flux is not a function of these four quantities but a function of the particles' distribution in the phase space. Nonetheless, the ability of neural networks to find non obvious patterns and correlations can be used to provide a good approximation of the heat flux, that can then be used as closure for the fluid model described in Section 2.2.

Neural networks enable to build functions relying on many parameters. The obtained closure can be formally written as:

q = C θ (ε, ρ, u, T ),
where θ denotes the set of parameters of the obtained network. Here the closure is chosen to be non-local: ρ, u, T as well as q are all spatial discretized quantities. The closure takes as input four vectors or, in the neural network terminology, a 1D signal with four channels corresponding respectively to the Knudsen number (turned into a constant vector), the density, the mean velocity and the temperature at spatial discretization points:

X = (ε, ρ, u, T ) ∈ (R Nx ) 4 ,
where N x is the number of spatial discretization points. The output of the closure is one vector of size

N x Y = C θ (ε, ρ, u, T ) ∈ R Nx ,
that will correspond to the estimated heat flux q ∈ R Nx .

To find such a function, there are two main steps:

1. First, by setting the architecture of the neural network, we define a family of functions C θ parametrized by θ ∈ Θ, where Θ represents the set of all the possible parameters of the neural network.

2. Then, by training the neural network, we find a set of parameters θ so as to minimize the error of the neural network predictions on a dataset, for which the true heat flux is known.

The neural network is the core component of the closure C θ , but it is not its only one. Indeed, for the data to be usable by the neural network and to provide satisfying results, it needs to go through a certain amount of processing. Section 3.2 introduces these other components and their articulation in the closure. Sections 3.3 and 3.4 then describe respectively the architecture and the training of the neural network.

Detailed composition of the closure

In practice, the closure is used at each iteration in time to compute the heat flux q on all the mesh from the input (ε, ρ, u, T ) on all the mesh. Such a use of a neural network based closure into a numerical scheme requires some work to transform the data handled by the solver into data usable by the network, and vice-versa. This section briefly introduces the different steps involved in the closure and their role in relation to the network.

There are three main processing operations designed to integrate the closure into the numerical scheme. The first one deals with the difference in resolution between the numerical scheme and the data used to train the network, by resampling the input of the network to its training resolution, and its output back to the original resolution. The second one deals with the difference in length between the resampled data of the scheme and the inputs handled by the network, by slicing the data into several pieces to give to the network, and then aggregating the outputs to reconstruct the full heat flux. The third one deals with the stability of the resulting numerical scheme, by smoothing the data from the network to prevent any oscillation to propagate into the numerical scheme. These three steps are part of a whole process, that can be broken downs as follows:

C θ : X (Re)+(P) -→ X (Sl) -→ ( Xj ) j (NN θ ) -→ ( Ỹj ) j (R) -→ Ỹ (P')+(Sm)+(Re) -→ Y.
This process is illustrated by figure 1, and the different steps are described below.

Resampling (Re) and pre-processing (P)

The input X of size N x is resampled with a Fourier method (that relies on the periodicity of the signal) and pre-processed into a signal X of size N x . As explained in section 6.3.1, the resampling changes the resolution of the input to match the one used for training and get better results. The pre-processing step that follows is a common step in the use of neural networks that allows to improve the accuracy. It is detailed in section 5.1.

Slicing (Sl)

The inputs are sliced into several pieces ( Xj ) j of size N before being handed over to the network, where N is the size of the input the neural network works with.

Neural network (NN θ )

For each piece of input Xj , the network predicts a piece of heat flux Ỹj , of size N .

NN θ NN θ NN θ ε, ρ, u, T q (Re)+(P) (Sl) (R) (P')+(Sm)+(Re')

Reconstructing (R)

The signals Ỹj are aggregated in order to reconstruct the whole heat flux. An added benefit of this slicing and reconstructing process is that it allows to prevent some edge effect introduced by the network: by using overlapping pieces, the ends of each piece can be ignored when reconstructing the whole signal. As described in section 5.3, the amount of overlapping is a parameter that can be tweaked.

Post-processing (P'), smoothing (Sm) and resampling (Re')

The post-processing is an operation of normalization designed to improve the accuracy of the network and detailed in section 5.2. The smoothing of the signal is here to avoid oscillations coming from the predictions to propagate and amplify, thus making the numerical scheme unstable. The intensity of this smoothing can be tweaked, and is discussed in section 6.2.4. Finally it is followed by a resampling to recover the original resolution, for the output to be used in the rest of the computations.

These different operations will be fully exposed in Section 5. We now turn to the heart of the closure: the neural network.

Architecture of the neural network

The closure C θ depends on the architecture of the network. In our work we use a 1D version of the V-Net [START_REF] Milletari | V-net: Fully convolutional neural networks for volumetric medical image segmentation[END_REF][START_REF] Ronneberger | U-net: Convolutional networks for biomedical image segmentation[END_REF] implemented with the Tensorflow 2 library through its python interface. The V-net is a fully convolutional neural network meaning that it is based only on a successive sequence of convolutions (no fully connected layers). It consists of two parts. First a descending part that acts as an "encoder", and decomposes the input X into multiple features, performing a multiscale analysis. Then an ascending part that acts as a "decoder", and synthesizes the created features to predict the output Y . As we will see below, it can be characterized by three hyperparameters: the number of levels , the depth d of the first convolution and the kernel size p of all convolutions. The composition of the V-net illustrated Figure 2 Input As described previously, it is a 1D signal X made of 4 channels (ε, ρ, u, T ), and with a given length called the window size that we set to N .

Initialization

We first perform successively two 1D convolutions 1 with a kernel size of p (we tried p = 5, 7, 9, 11), both followed by the softplus activation function s(x) = ln(1 + e x ), and change the number of channels from 4 to a depth d (we tried d = 4, 5, 6, 8). All convolutions preserve the size of the signal with the help of a constant padding that extends the signal by continuity. The original V-Net uses the ReLu activation function, but the regularity of the softplus activation function seems to give us better results.

Descending part

The V-net is made of levels (we tried = 3, 4, 5), so -1 descents and as many ascents. Let us describe one descent.

• The input is down-sampled by a convolution of kernel size 2, stride 2, and that doubles the depth of its input. The stride allows to halve the length of the input, incidentally blurring the precise locations of the highlighted features.

• Then are applied two successive convolutions of kernel size p and that conserve the depth of their input, both followed by the softplus function.

• The output of this double convolution is added to its input. This "shortcut" is represented by the dotted arrows in Figure 2. This way the convolutions produce additive (or residual) modifications. This is the principle of the famous Resnet. It has been shown that it accelerates the training process and limits the problems of vanishing or exploding gradient [START_REF] He | Deep residual learning for image recognition[END_REF].

The output of this descent is a signal with half the length and double the depth of the input signal.

Ascending part It works as a mirror of the descending part. Each ascent simultaneously increases the length of the signal and decreases its depth by a factor of 2. Technically, we up-sample the data using a transposedconvolution [START_REF] Dumoulin | A guide to convolution arithmetic for deep learning[END_REF] of size 2 followed by two convolutions of kernel size p with softplus activations.

Long shortcuts At each level, the features obtained in the ascents are combined with the features created by the descents, using a summation. This

1 For the sake of completeness, we recall the formula for a one dimensional convolution: applied on an input X of shape (N, d) to get an output Y of shape (N, d ), a 1D convolution with an odd kernel size p uses a kernel K of shape (p, d, d ) and we have

Y i,k = d j=1 p di=1 Xi+di,j K di,j,k ,
where X is the input X padded on both ends so that Y i, k is well defined for i up to N . allows to retrieve some fine resolution details, which are lost by the blurring produced by down and up samplings.

Weighted mean

The resulting signal of depth d is turned into a one dimensional output Y by simply taking a weighted mean of all d channels. This operation is implemented with a convolution of kernel size 1. No activation function is added, as we deal with a regression problem.

All the coefficients of the convolutions constitute the set of parameters θ of the neural network. Table 1 summarizes the hyper-parameters we choose in practice as a reference. This choice was mostly motivated by the results given in Section 6.1.2. In order to avoid the oscillations that can be produced by the upsampling oper-

Hyper

Training of the neural network

The training consists in finding an optimal set of parameters θ for the neural network. It is defined as the minimizer of a loss function that measures the error of the network: θ = argmin θ∈Θ Loss(θ).

For this loss function we choose the mean absolute error (MAE) between the output of the network and the expected heat flux over a training dataset:

Loss(θ) = j∈training dataset |NN θ ( Xj ) -Ỹj |,
where ( Xj ) (resp. Ỹj ) are data of size N × 4 (resp. N ), obtained from kinetic simulations after pre-processing (P) and slicing (Sl) of vectors (ε, ρ, u, T ) (resp. q). In practice, the minimization is carried out with a gradient descent algorithm. This whole process falls under the category of supervised learning, as it requires a labelled dataset, including both the inputs and the corresponding expected outputs. The first step to train the neural network is thus to generate such a labelled dataset. To do so, we consider the kinetic model ( 5)- [START_REF] Cai | Globally hyperbolic regularization of grad's moment system[END_REF], that describes the particles' distribution from which can be derived their density, mean velocity, temperature, and heat flux. We run many simulations with a given time step and given number of discretization points in space N x and velocity N v , with the numerical method described in appendix A. We save the results in our dataset at selected times.

The making of this dataset is described in more details in Section 4. In particular, it requires crucial choices in the initial distributions, the recording times and the range of the parameter ε. It has a direct impact on the range of validity of the closure obtained and it will be numerically discussed in Section 6.

In the end, the network is trained on a dataset of 10 000 entries (

X k ; Y k ) = (ε k , ρ k , u k , T k ; q k )
, that are processed as described in Section 5. In particular, each vector (of size N x = 1 024) is sliced into overlapping windows of size N = 512, resulting in 80 000 inputs and labels. 4% of these are isolated as a small test dataset and 10% of what remains are used for validation. The rest is used for the actual training.

This training consists in 5 series of 120 epochs with a decaying learning rate, reset to its initial value of 0.005 between each series. We use the mini-batch gradient descent method, with the Adam optimizer and mini-batches of size 1 024, to minimize the mean absolute error (MAE).

Data generation

In order to train a neural network to predict the heat flux q knowing the Knudsen number ε, the density ρ, the mean velocity u and the temperature T , we need to generate a training dataset. This dataset must meet two criteria. First, it must be a labelled dataset, i.e. a dataset with both the input X = (ε, ρ, u, T ) and the expected output Y = q. The expected output is the "correct" output of the given input, from which we want the network to interpolate, or generalize, for new inputs. Then, this dataset must contain as much diversity as possible, for the interpolation to be usable in many different situations.

Consequently, the generation of the training dataset relies on two key ingredients. First, a solver for the kinetic model that can reliably estimate the distribution of particles in the phase space, from which can be derived all the fluid quantities ρ, u, T and q. Second, a mechanism to produce a variety of distributions that can be given to the solver as initial solutions. This whole process of data generation is illustrated in Figure 3 and described in more details in this section ; except for the solver for the kinetic model that is not specific to our work and is described in appendix A.

Global description

Let us first get a glance at the global process illustrated by Figure 3. To begin with, a randomized process allows us to produce a Knudsen number ε, an initial solution f 0 , and a recording time t. These three elements are given to the kinetic model that computes the evolution of the distribution up to t. The resulting distribution f is then used to compute the density ρ, the mean velocity u, the temperature T and the heat flux q of the system. On top of these four quantities, we add the Knudsen number ε as a constant vector of the same size as the other four, so that it can be used as input by the network. This could also allow us to generalize our method to a problem with a Knudsen number that varies in space. In the end, these five vectors (ε, ρ, u, T ; q) form an entry in the dataset. Actually, in order to capitalize on the simulations that can be computationally expensive -especially in higher dimensions-, we decide to produce several entries with each simulation. For a given Knudsen number and a given initial solution, we generate not one but 20 different recording times t 1 < • • • < t 20 , and compute the distribution at each one of these times, respectively f 1 , ..., f 20 . From each one of these distributions can be derived ρ, u, T , q, to which we prepend ε to form 20 different entries to be stored in the dataset. Thus, for the cost of one simulation up to t 20 , we produce 20 different entries corresponding to different intermediate steps in the simulation.

To build a whole dataset, we repeat this process with 100 different values of ε, and with 5 different randomly selected initial solutions for each one of these ε. Thus, we get a dataset with 10 000 entries, based on 500 different initial solutions. We produce two such datasets: a training dataset used to train the neural networks, and a test dataset used to measure and compare the performance of the trained networks on new data.

Random initial conditions

The initial solution f 0 is a Maxwellian as in [START_REF] Chew | Application of dispersion relations to low-energy meson-nucleon scattering[END_REF], with density ρ, mean velocity u and temperature T randomly generated as partial Fourier series under the form

α × a 0 2 + 0.5 N n=1 (a n cos(nx) + b n sin(nx)) , x ∈ [0, 2π].
We set N to 20 and we randomly generate the a n and b n coefficients for n 1 according to a uniform distribution on -1 n , 1 n . This choice for the a n and b n coefficients set the high frequencies to low amplitudes, resulting in more regular functions. Then, the choice for α and a 0 depends on the function to be generated. For the density, a0 2 = 1 and α = 1. For the temperature, a0 2 = 1 and α is chosen uniformly in [0.1, 1]. For the mean velocity, a0

2 is chosen uniformly in [-1, 1], and α is chosen so that the maximum Mach number max

x∈[0,2π] |u(x)| 2T (x) ,
falls uniformly in [10 -4 , 5 • 10 -1 ] in logarithmic scale. We also make sure that the density and temperature generated are positive functions. Possible functions generated with this process are given Figure 4. For their discrete form, these functions are sampled with N x = 1 024 points to get vectors, that are turned into a discrete Maxwellian our solver can work with, with N v = 141 points is velocity.

Random Knudsen numbers and recording times

In this work we focus on Knudsen numbers in the range [ε min , ε max ] = [0.01, 1]. First because below 0.01 the fluid model with the Navier-Stokes estimation is already very accurate, and also because we do not want a range too wide. Regarding the distribution of the values in this interval, we want a decent amount of entries with ε really close to 0.01, but we do not want too few entries for ε close to 1. As a result, a uniform or logarithmic distribution over the interval [0.01, 1] is not satisfying. Instead we opt for a uniform distribution for

√ ε in the interval [ √ ε min , √ ε max ]
. This allows us to have about one fourth of the Knudsen numbers between 0.01 and 0.1, and three fourths between 0.1 and 1. Also we use a deterministic distribution for the train dataset in order to have a nice coverage of the interval, while we use a random distribution for the test dataset in order to test the ability of the trained networks to generalize to new Knudsen numbers.

For the recording times t 1 < • • • < t 20 , we do not want t 1 to be too close to zero as the initial maxwellian state always has a heat flux of zero and does not carry much physical information, and we do not want t 20 to be too big so that the simulations do not last too long. For these reasons we choose 20 times uniformly in the interval [0.1, 2], that are then sorted and given successively to the solver for the kinetic model.

Computing of the fluid quantities

Once the Knudsen number, the initial solution and the times are generated, a simulation is computed with the numerical method described in appendix A. It relies on a discretization of the phase space

(x i , v j ) , i ∈ {1, ..., N x }, j ∈ {1, ..., N v },
and results in an approximation (f i,j ) of the real distribution f :

f i,j f (x i , v j ).
From this distribution can be derived the density ρ, mean velocity u, temperature T and heat flux q needed for the dataset, as mentioned in Section 2.1. In their discrete form, these quantities are vectors computed as follows: for any

1 i N x , ρ i = Nv j=1 f i,j , ρ i u i = Nv j=1 v j f i,j , ρ i T i = Nv j=1 (v j -u i ) 2 f i,j
and

q i = Nv j=1 (v j -u i ) 3 f i,j .
For our datasets, we use a discretization with N x = 1 024 points in space, and

N v = 141 points in velocity. Training Network Network Network Loss ε, ρ, u, T q (S) (Sl) (Sl) (N) Prediction Network Network Network ε, ρ, u, T q (S) (Sl) (R) (I) (Sm)

Data processing

In this section we describe how the data is processed on both sides of the neural network, both for the training process and the prediction process. At this point we have generated some data with the kinetic model, and want to process it for a neural network that takes as input a 1D signal with four channels -the Knudsen number, the density, the mean velocity and the temperature-, and outputs a 1D signal with one channel -the heat flux. This data processing consists in the transformations described below and summarized in Figure 5.

In the remainder of this section, we use the word standardized to describe a quantity with mean zero and unit variance, and normalized a quantity with norm one or that has been applied the transformation described in Section 5.2.

Input standardization

Neural networks are known to be easier to train with standardized inputs, where all channels have the same order of magnitude. Consequently, we standardize the inputs channel by channel, by removing the mean and scaling to unit variance the whole training dataset. For example for the density ρ, if (ρ k ) is the family of densities of the whole training dataset, this standardization applied on a given entry ρ k0 can be written

ρ k0 standard = ρ k0 -mean k,i ρ k i std k,i ρ k i ,
where mean is the empirical mean and std the empirical standard deviation.

Since the neural network is trained with this standardized data, it is designed to work only with inputs that received that specific standardization: after the training, any new input given to the neural network must be applied that very same standardization. In particular, the means and standard deviation computed in the training dataset have to be stored with the neural network, in order to be available later to make predictions for the fluid model.

Output normalization

To prevent the outputs in our datasets from being smaller than the typical error of our neural network, we choose to train the network to predict normalized output. Otherwise, the network would produce pure noise when trying to predict them, which turns to be problematic when used by the fluid model, especially regarding its stability. Another solution would be to use a cost fonction measuring a relative error instead of an absolute one but it turns out to be less efficient. Let us describe how we proceed.

The main issue when training the network with normalized output is the reversibility of the normalization. Since in the end we want to predict a heat flux, and not a normalized heat flux, we need to be able to apply an inverse normalization at the output of the neural network. As a consequence, this transformation cannot use any knowledge on the heat flux to be predicted, that would be available in the labelled training dataset, but would not in a real case scenario. For this reason, we choose to normalize the heat flux with its estimation given by the Navier-Stokes approximation [START_REF] Crouseilles | A hybrid kinetic/fluid model for solving the gas dynamics boltzmann-bgk equation[END_REF], and that can be computed from the Knudsen number, the density and the temperature, all available in a real case scenario.

But using an estimation of the heat flux brings another concern: when this estimation is way off, the normalization (or inverse normalization) would distort the information a lot. And this is expected to happen, as the Navier-Stokes approximation tends to greatly overestimate big heat fluxes with Knudsen numbers in the range we consider. To prevent this from happening, we choose to normalize heat fluxes only when the corresponding Navier-Stokes estimation is below a given threshold θ, that we set to 0.1. This threshold is also a way of only addressing the heat fluxes that were problematic in the first place.

Thus, for a given input (ε k0 , ρ k0 , u k0 , T k0 ) in the training dataset with expected output q k0 , the normalization we use can be written

q k0 norm =          q k0 × θ q k0 N S , if 0 < q k0 N S θ, q k0 , otherwise,
with

q k0 N S = max i=1,...,N 3 2 ε k0 ρ k0 i (∂ x T ) k0 i ,
where i refers to the index of the vectors and the spatial derivative is approximated with a centered finite difference formula. Note that we choose not to normalize at all heat fluxes such that q N S = 0. Once trained to predict these normalized heat fluxes, the neural network must be followed by an inverse normalization when used to make predictions of non normalized heat fluxes. This inverse normalization simply reads

q k0 =        q k0 norm × q k0 N S θ , if 0 < q k0 N S
θ, q k0 norm , otherwise.

.

Slicing and reconstructing

In order to be usable with meshes of different sizes, the neural network is designed to work with only one portion of the signal at a time as input, and to return the corresponding portion as output. As a consequence, each signal given to the network has first to be sliced into several windows of the right size. For the training process, both the input and the expected output have to be applied the same slicing. For the prediction process, the input has to be sliced into windows, and the resulting predictions have then to be aggregated to reconstruct a complete signal.

In one case or the other, we slice the signal into overlapping windows. The overlapping serves two purposes. For the training process, it allows us to artificially increase the amount of data we can give to the network from the training dataset. Though it introduces redundancy, such data augmentation is known to help prevent overfitting and improve the accuracy of the network [START_REF] Shorten | A survey on image data augmentation for deep learning[END_REF]. For the prediction process, we use this redundancy introduced by the overlapping when reconstructing the output signal, to provide a first smoothing process and to dismiss the defects on the borders of the predictions.

These defects come from the necessary padding in our V-Net architecture, that can cause oscillations of big amplitudes at the ends of the predictions. As a consequence, each window must have a margin that will be ignored when reconstructing the signal. We set this margin to 10% of the output on each side, leaving 80% actually useful for the reconstruction, later referred to as the useful part. This observation on its own is enough to require some overlapping between the windows to reconstruct the whole signal, but we introduce even more overlapping. We slice the original signal such that each one of its points is found in a fixed number r of windows, margins aside. We call this number the redundancy parameter. The bigger it is, the bigger the overlapping, and the more windows are produced. This slicing is illustrated Figure 6.

For the reconstruction of an entire signal from the predictions of the network, we get rid of the margins and multiply the useful parts by the kernel

1 n cos 2π L u x -π + 1 , x ∈ [0, L u ]
where r is the redundancy parameter and L u the length of the useful part in the physical space. We then sum up the resulting windows. This kernel has the property to sum up to one when duplicated every Lu r (see Figure 7). As a consequence, each value in the reconstructed signal is a weighted mean of the r values predicted by the neural network, which results in a somewhat smooth estimation of the heat flux. 0

L u 0 1 0 Lu 3 2Lu 3 L u 0 1 sum
Figure 7: Kernel used for the reconstruction of an entire signal from the windows predicted by the neural network, with the redundancy parameter r = 3.

Smoothing of the output

At this point we have a process that allows us to estimate the heat flux q with the neural network, and that can be used in a fluid model. But as it is, even with the smoothing provided by the aggregation of different predictions as described in the previous section, this estimation can show some small irregularities. In the fluid model, these irregularities on q are amplified by the computation of ∂ x q, and can cause the instability of the hyperbolic system. This instability cannot be prevented by refining the time discretization, we have to control the irregularities in the estimated heat flux. To do so, we smooth the reconstruted heat flux thanks to a convolution with a gaussian kernel. This kernel with standard deviation σ reads

w : t ∈ [-3σ, 3σ] → e -1 2 t 2 σ 2 3σ -3σ e -1 2 t 2 σ 2 dt
, and the smoothed heat flux is

q(x) = 3σ -3σ q(x + t)w(t) dt.
We show in Section 6.2.4 that σ can be chosen relatively large (σ 0.05) without impacting the accuracy of the estimation, and has to be for the method to be stable.

Results

This section is divided in three parts: first we introduce the results regarding the neural network only, then when it is used in the fluid model, and finally how the whole model performs in configurations that differ from the one used for the training. Unless otherwise stated, the hyperparameters of the neural network and the training are those presented in section 3.3.

Accuracy of the network

In this section we show numerical results of the neural network, independently of its use with the fluid model. To this end we use the test dataset with entries (ε, ρ, u, T ; q) computed with the kinetic model, or directly compare the predictions with a simulation of the kinetic model. To measure the accuracy of the network on a given entry, we use the relative error in norm

L 2 ||q -q|| 2 ||q|| 2 ,
where q is the prediction of the neural network with inputs (ε, ρ, u, T ). We also observed the relative error in norm L ∞ , but it was essentially the same, that is why we focus on the L 2 norm in what follows. Unless stated otherwise, the results given use a redundancy parameter of 2 and a smoothing of σ 0.06. This last choice is motivated by the results of Section 6.2.4.

Comparison with the Navier-Stokes estimation

Figure 8 shows some examples of predictions, compared to the real heat flux and the estimation of Navier-Stokes [START_REF] Crouseilles | A hybrid kinetic/fluid model for solving the gas dynamics boltzmann-bgk equation[END_REF]. The first example illustrates how the Navier-Stokes estimation tends to perform better on data with small Knudsen numbers. The second and third example show how on the contrary it very often overestimate the heat flux when the Knudsen number increases.

The relative errors of both the network and the Navier-Stokes estimation over the whole test dataset are summarized in the histogram Figure 9. On this dataset, the overall error of the neural network is about an order of magnitude below that of the Navier-Stokes estimation.

Figure 10 gives a closer look at these errors and highlights two important factors in the performance of the network over the Navier-Stokes estimation: the Knudsen number and the norm of the real flux. The first scatter plot (Fig. 10, left) shows that the error of the Navier-Stokes estimation heavily depends on the Knudsen number ε, which is not the case for the neural network predictions, except for small values. This is even better illustrated by Figure 11, that uses the same data but where the error is not set in logarithmic scale. The second scatter plot (Fig. 10, right) shows that there is a strong correlation between the norm of the real heat flux and the relative error of the network: the smaller the heat flux, the bigger the error of the network. This could be explained by a lack of small heat fluxes in the training dataset or by the normalization we chose for small heat fluxes. In any case, it should not be that big of a problem when used with the fluid model, as smaller heat fluxes have a smaller impact on the result, at least as long as they are quite smooth.

It can also be interesting to look at the evolution of the error during a simulation, to get a better idea of what will happen when used with the fluid model. To do so, we choose 50 random initial solutions and 50 random Knudsen numbers in the same way we did to build the test dataset, and run the kinetic model 0 with these data. At regular times, the heat flux of the kinetic model is compared to the heat flux the network would have predicted and to the Navier-Stokes estimation. Figure 12 shows the statistical results of this experiment. We observe that the network has a significant advantage over the Navier-Stokes estimation during the first 2.5 time units, before decreasing to a smaller advantage. This result can be explained by the fact that the heat flux is bigger during the first three time units (the dissipation seems decrease the size of q in time), and by the strong correlation between the norm of the heat flux and the error of the network, also visible in this figure.

π 2 π 3π 2 2π -0.02 0 0.02 ε = 0.019 q(x) Real Network Navier-Stokes 0 π 2 π 3π 2 2π -0.2 0 0.2 ε = 0.098 q(x) 0 π 2 π 3π 2 2π 0 0.5 ε = 0.555 x q(x)

Hyper-parameters of the neural network

In this section we show the results that motivated the choice for some hyperparameters of the V-Net. We focus on the number of levels, the depth d, and the size p of the kernels of the convolutions. Figure 13 compares the median relative error of V-Nets with different sets of hyper-parameters. For these results the window size was set to 256, but we later decided to set it to 512 as it gave slightly better results.

We observe that all hyper-parameters do not affect the performance of the V-Net in the same way. For instance, decreasing the size of the kernels significantly decreases the accuracy of the trained network, while not decreasing the number of parameters as much. On the other hand, decreasing the depth allows to decrease the number of parameters with a very small effect on the accuracy. In the remainder of this paper we work with the "l5,d4,p11" architecture, as it seems to be a good compromise between the number of parameters and the accuracy of the network. 0.115 0.120 0.125 0.130 0.135 0.140 0.145

1 • 2 • 3 • 4 • 5 • 6 • 7 • 5,d8,p11 4,d8,p11 3,d8,p11 5,d6,p11 5,d5,p11 5,d4,p11 5,d8,p9 5,d8,p7 5,d8,
• 2 • 3 • 4 • 5 • 6 • 7 • 5,d8,p11 4,d8,p11 3,d8,p11 5,d6,p11 5,d5,p11 5,d4,p11 5,d8,p9 5,d8,p7 5,d8,p5
L ∞ relative error Number of parameters Figure 13: Median relative error and number of learnable parameters for V-Nets with different hyper-parameters : ( ) number of levels, (d) depth, and (p) size of the kernels. Points sharing the same color represent networks that differ by only one hyper-parameter.

Fluid model with the neural network

In this section we look at how the neural network performs when used in the fluid model. We denote by "Fluid+Network" this method. The fluid model is solved using an explicit finite volume method as presented in appendix B.1. We compare it to three others: Kinetic : the kinetic model. It is the most accurate model and serves as a reference for the real target. The numerical method is described in appendix A.

Fluid+Kinetic : the fluid model with the heat flux from the kinetic model. It is the result we would get with a perfect neural network that makes no error, and helps to distinguish the error of the fluid model from the error of the neural network on the heat flux. The numerical method is identical to the one of the Fluid+Network method (see appendix B.1). Theoretically this model should give the same result as the kinetic model, but since the numerical schemes and viscosities are different it is not always true in practice.

Navier-Stokes : the fluid model with the Navier-Stokes estimation. It does not use the same numerical method as our model, since the formula for the heat flux requires an implicit scheme to avoid too stringent stability condition (see appendix B.2).

Our criteria to compare the fluid models is the L 2 relative error on the logarithm of the electric energy E

E(t) = [0,2π] E(x, t) 2 dx,
compared to the kinetic model. For all the following tests we use a discretization of N x = 512 points on [0, 2π] in physical space, and N v = 101 points on [-7, 7] in velocity space (for the kinetic model). In particular, the data is resampled on 1 024 points at each iteration for the neural network as explained in Section 3.2. The simulations are computed up to t = 8.

Examples

Figure 14 shows examples of electric energies obtained with the four models described above. The first one uses a very small Knudsen number, and with no surprise the Navier-Stokes model performs better than the fluid model with the network. On the second example on the other hand, the Knudsen number is bigger and we observe that the oscillations are slightly shifted: Navier-Stokes shows some dispersion that the fluid model with the network does not. This dispersion often increases as the Knudsen number gets bigger as it shows on the third example. Finally the fourth one illustrates a case where the Navier-Stokes model really struggle while the fluid model with the network does pretty good.

Global performances

The L 2 relative errors of the three fluid models over 200 simulations with different initial solutions and different Knudsen numbers are summarized in Figure 15. We observe that our closure with the neural network does not reach relative errors of 0.01 or below, contrary to the closure with the real heat flux or even with the Navier-Stokes estimation. But it successfully maintains the relative error below 0.2, which is close to what the closure with the real heat flux achieves, and much better than the Navier-Stokes estimation that often exceeds 0.2. As already noted, the Fluid+Kinetic relative error should be zero but is not, due to numerical diffusion errors.

Influence of the Knudsen number

In Figure 16 close to the fluid model with the kinetic heat flux. This would seem to indicate that the network appropriately plays its role and that its error on the heat flux does not impact the whole model too much.

Smoothing of the prediction and stability

In this section we show the impact of the smoothing on the accuracy of the predictions, as well as its impact on the stability of the fluid model using these predictions. Figure 17 gives a visual example of how the smoothing modifies a prediction of the neural network.

To plot Figure 18, the network was used to predict the heat flux of each entry in the test dataset, and different quantities of smoothing were applied before computing the relative error. We observe that the accuracy of the resulting heat flux does not deteriorate for standard deviations lower than σ 0.04, and then slightly decreases as σ increases.

This smoothing was introduced to cope with the instability of the fluid model relying on the predictions of the network. To measure how this stability is improved by smoothing the output of the network, we choose 30 random initial solutions as we did to build the test dataset and 30 Knudsen numbers uniformly distributed between 0.01 and 1. Then we choose a set of smoothing parameters σ, and for each one of these we run the 30 simulations with the fluid model, using the network with this quantity of smoothing. Finally, for each σ we look at the proportion of simulations that reached t = 3, from which we assume the model will remain stable. The results are shown Figure 19. It appears that all 30 simulations reach t = 3 for σ above approximately 0.05. We use σ 0.06 for our tests with the fluid model shown in the next section, and it appears to be enough since no simulation failed so far with this quantity of smoothing.

Using the network with different configurations

In this section we are interested in the flexibility of our approach to different configurations. First we take a look at different discretizations of space to measure the ability of the network to generalize to other resolutions than the one it has been trained with. Then we introduce discontinuities in the initial solutions to see how the fluid model with the network reacts in terms of stability and accuracy.

Different resolutions

The network was trained with data using 1 024 points on [0, 2π] in space. The fact that it uses windows of size 512 would allow it to be used with any discretizations with at least 512 points, or even less if we use the periodicity to make windows of size 512. However, with initial solutions of the same form than those used to build the datasets, the change in resolution modifies the way the data is interpreted by the neural network. For instance if we sample an initial condition on 512 points instead of 1 024, the signal frequencies would be multiplied by two in the eyes of the neural network. This can result in inputs differing from the training inputs, and a decrease in accuracy is to be expected.

One way to measure this effect is to evaluate the network on the test dataset but resampled at different resolutions, to mimic the data that could be generated by the fluid model with these resolutions. We can then compute the relative errors as in Section 6.1. The results Figure 20 show that the accuracy of the network does indeed decrease as the resolution of the data moves away from its To avoid this issue, if the fluid model uses a different resolution it is best to resample the data for the network to match the resolution of the training data, as this resampling operation only introduces a negligible error. However, this phenomenon shows some of the limits of the network in terms of generalization. With initial solutions that differ from the ones used to train the network -in frequency for instance-, we could expect to see a significant loss of accuracy.

Discontinuities

The datasets were built from simulations starting with a continuous initial solution, but using the network in the fluid model with discontinuous initial solutions could be an option. In practice, the initial discontinuities quickly fade away in this physical system, but it might be enough for the predictions of the network to cause some oscillations and the instability that goes with it.

To see if that is the case, we run multiple simulations with different Knudsen numbers and different initial solutions, similar to those used to build the dataset, except that each function (density, velocity and temperature) can be multiplied by a function

x ∈ [0, 2π] → c π (x -x d ) + (1 + c), if x < x d , c π (x -x d ) + (1 -c), if x x d , adding a discontinuity at a random location x d ∈ [0, 2π] with a random ampli- tude c ∈ [-1, 1]
. With these initial solutions, no simulation failed to reach t = 8, and as can be seen Figure 21, the results are comparable to those presented in Section 6.2.

Conclusion

We construct a fluid closure for the Vlasov-Poisson dynamics based on V-net neural network and supervised learning from kinetic simulations. Slicing process of the data is introduced in order to manage meshes of different sizes. Several data processing have been also designed to improve the quality and regularity of the heat flux estimation. The numerical results show that the closure predicts the heat flux with a uniform relative error on the Knudsen interval [0.01, 1], while the Navier-Stokes closure does not as expected. We also observe that the prediction is better at the beginning of the numerical simulations where the distribution function farthest from the equilibria set and so the real heat flux is larger. Surprisingly, we numerically observe that the neural network closure does not introduce instabilities when inserted in the fluid simulations, provided, however, that the outputs are regularized. Finally, the closure is quasi-optimal as the relative error between a full kinetic and a Fluid+Network closure behaves like the one between a full kinetic and a Fluid+Kinetic closure. This work raises several issues from the numerical point of view. The first question is its efficiency in terms of computing time or its energy cost. Presently, the Fluid+Network model requires about the same computing time as the kinetic model. However, we can hope that it will be more efficient in higher dimension, as the number of computations for a V-Net with levels, depth d and kernel size p in dimension 1 is about O(2 d 2 pN ), while it is about O( d 2 p 2 N 2 ) in dimension 2 and O(d 2 p 3 N 3 ) in dimension 3. Compared with the O(N m N m v ) computations required for the kinetic model where m denotes the space dimension, it increases much more slowly since p << N v .Moreover, the neural network approach can greatly benefit from GPU parallelism.

The stability of the Fluid+Network closure is also an important issue. Despite the good numerical results obtained, a mathematical guaranty is lacking. Constructing neural networks closure ensuring such stability remains to be done.

Finally, we plan to apply this method to both the Vlasov-Poisson system in higher dimension and also to other closure problems arising in the MHD or gyrofluid design.

A Numerical scheme for the kinetic model

In this appendix we describe the numerical method used to solve the Vlasov-Poisson equations ( 5)-( 6) resulting from the kinetic model. This numerical method is used to produce data that can in turn be used by the neural network to interpolate the heat flux. It also serves as a reference to compute the error of the other methods, allowing us to compare them. For better readability, let us remind the one dimensional Vlasov-Poisson equations:

∂ t f + v ∂ x f -E ∂ v f = 1 ε (M (f ) -f ), E = -∂ x φ, ∂ xx φ = ρ - L 0 ρ dx.
The spatial domain is given by [0, L], where L > 0 is the spatial length. For numerical purpose, the velocity domain is restricted to the bounded interval [-v max , v max ]. Thus we complement the equation with the following boundary conditions: (±E) -f = 0, at v = ±v max .

We consider a time discretization (t n ) n with variable time step ∆t and a discretized phase space (x i , v j ) i,j with constant steps ∆x and ∆v respectively. The number of discretization points in space (resp. in velocity) is denoted N x (resp.

N v ). We denote by f n i,j the approximation

f n i,j f (x i , v j , t n ).
We also use the notations f n for the matrix (f n i,j ) i,j , f n i for the vector (f n i,j ) j and f n j for the vector (f n i,j ) i .

A.1 Time discretization

To solve the Vlasov-Poisson equations over the time interval [t n , t n + 1] we use a splitting between three stages:

1. Compute the electric field at time t n by solving the Poisson system:

E = -∂ x φ , -∂ xx φ = 1 L L 0 ρdx, 2.
Transport the distribution function by solving the Vlasov equation over the time inverval [t n , t n + 1]:

∂ t f + v∂ x f -E∂ v f = 0,
3. Update the distribution function by taking into account the collision operator:

∂ t f = 1 ε (M (f ) -f ).
The first two stages rely on the space discretization and are discussed in the next section. For the third operator with a stiff source term, we use an implicit scheme:

f n+1 -f n ∆t = 1 ε (M (f n+1 ) -f n+1 ).
Knowing that the fluid quantities ρ, u, T are preserved by this operator, we have M (f n+1 ) = M (f n ), so the scheme can be rewritten as:

f n+1 = f n + ω (M (f n ) -f n ), with ω = ∆t ∆t + ε . ( 10
)

A.2 Spatial discretization

For the spatial discretization, we propose a method introduced in [START_REF] Pham | Space-only hyperbolic approximation of the vlasov equation[END_REF][START_REF] Helluy | Reduced vlasov-maxwell simulations[END_REF]. First we discretize in velocity with a centered finite difference scheme. After discretization, we get the following hyperbolic system:

f n+1 -f n ∆t + Λ∂ x f n + EB(f n ) = 0 (11)
with Λ the diagonal matrix of velocities and B(f n ) a vector given by

B(f n ) j = - f n j+1 -f n j-1 2∆v
, j ∈ {2, . . . , N v -1},

B 1 = - 1 2 max(E, 0)Ef n 1 , B Nv = 1 2 min(E, 0)Ef n Nv
The choice of the boundary terms allows to obtain a dissipative hyperbolic system and to ensure L 2 stability. Then we discretize in space the hyperbolic system [START_REF] Degond | A multiscale kinetic-fluid solver with dynamic localization of kinetic effects[END_REF]. We use a finite volume scheme with an upwind flux:

f n+1 i -f n i ∆t + f n i+ 1 2 -f n i-1 2 ∆x -EB(f n ) = 0, (12)
with

f n i+ 1 2 = 1 2 Λ(f n i+1 + f n i ) -1 2 | Λ | (f n i+1 -f n i ).
The Poisson equation is solved using a classical finite difference scheme:

E n i = - φ n i+1 -φ n i-1 2∆x , - φ n i+1 -2φ n i -φ n i-1 ∆x 2 = ρ j - 1 N x N i ρ i . ( 13
)
where the density ρ j is computed as follows: ρ i = Nv j=1 f i,j . The total scheme is given by ( 13)-( 12)- [START_REF] Degond | Macroscopic limits of the boltzmann equation: a review[END_REF]. The time step is chosen such as to satisfy the following stability condition: ∆t min ∆x

v max , ∆v max i |E i | .
Note that it is only first order accurate in order to avoid dispersive oscillations in the numerical solutions.

B Numerical scheme for the fluid models

In this section we introduce the numerical methods for solving fluid models [START_REF] Crestetto | Kinetic/fluid micro-macro numerical schemes for vlasov-poisson-bgk equation using particles[END_REF]. The time discretization depends on the considered closure. For the neural network based (Fluid+Network) or the kinetic one (Fluid+Network), we use an explicit scheme. For the Navier-Stokes closure (Navier-Stokes), an implicit scheme is required to avoid a too stringent stability condition. In one case or the other, the Poisson equation is solved at the beginning of each iteration in time to compute the electric field, exactly as in Eq. ( 13). This section therefore focuses on the fluid equations. The following schemes are classical. We briefly present them for the sake of completeness.

B.1 Explicit scheme for the neural network or the kinetic closure

Fluid equations ( 8) can be written

∂ t U + ∂ x F(U) = -EH(U), (14)
with U = (ρ, ρu, w), F(U) = (ρu, ρu 2 + p, wu + pu + q), H(U) = (0, ρ, ρu) and where the heat flux q is given by the neural network based closure (Fluid+Network) or the one obtained from full kinetic simulations (Fluid+Kinetic). We solve it on the spatial domain [0, L].

When q = 0 (Euler closure), equation 14 is an hyperbolic system that can be solved using a finite volume method with a local Lax-Friedrichs numerical flux and an explicit scheme in time. The hyperbolic system has characteristic speeds equal to u, u + c, uc where c = 3p/ρ is the sound speed. When considering an additional heat flux, we propose to use the same scheme and just add a centered approximation of the heat flux in the numerical flux, as explained below.

Like for the kinetic model, we consider a time discretization (t n ) n with variable time step ∆t and a discretized phase space (x i ) i with constant step ∆x. We denote by U n i the approximation

U n i U(x i , t n ).
The finite volume scheme results in the following formula: are chosen to be larger than the maximum characteristic speed of the hyperbolic system, as u n i is the local speed of particles and c n i the local sound speed. Comparing with the classical scheme for Euler system (q = 0), the numerical flux for the momentum have an additional term equal to 1 2 (q n i+1 +q n i ). The time step is chosen such as to satisfy the CFL stability condition:

U n+1 i -U n i ∆t + F(U) n
max i S n i+ 1 2 ∆t = 1 2 ∆x. ( 15
)
We refer to [START_REF] Leveque | Finite volume methods for hyperbolic problems[END_REF] for more details on finite volume methods.

For the Fluid+Network method, the heat flux q involved in the flux term is computed from ε, ρ, u and T at each iteration. For the Fluid+Kinetic method, it is obtained from an underlying kinetic simulation, using the same initial condition at t = 0. Note that the stability condition [START_REF] Dumoulin | A guide to convolution arithmetic for deep learning[END_REF] does not take into account the non-zero heat flux. However, as observed in Section 6.2, this term does not result in an additional stability condition for the Fluid+Kinetic and the Fluid+Network method, provided for the latter that the heat flux is sufficiently smoothed out.

Figure 1 :

 1 Figure 1: Graphic of the composition of the closure. The different operations are (Re)-(Re') resampling, (P) pre-processing, (Sl) slicing, (NN θ ) neural network, (R) reconstruction, (P') post-processing and (Sm) smoothing.

Figure 2 :

 2 Figure 2: Graphic of a 1D V-Net with a window of size N = 512, a depth k = 2 and = 3 levels.

Figure 3 :

 3 Figure 3: Scheme of the data generation process for one simulation with the kinetic model. With one Knudsen number and one initial solution, we produce 20 entries in the dataset.

Figure 4 :

 4 Figure 4: Possible functions ρ, u and T generated with the process described in Section 4.2.

Figure 5 :

 5 Figure 5: Scheme of the whole data processing for the training and the prediction processes respectivily. The different operations are (S) standardization, (Sl) slicing, (N) normalization, (I) inverse normalization, (R) reconstruction and (Sm) smoothing.

3 Figure 6 :

 36 Figure 6: Slicing of a 1D signal into overlapping windows with different redundancy parameters r. The original signal is represented vertically on the left, and the windows are on the right. The hatched parts of each window are the margins that are ignored when reconstructing the signal. The windows that exceed the span of the original signal are completed by periodicity.

Figure 8 :Figure 9 :Figure 10 :Figure 11 :Figure 12 :

 89101112 Figure 8: Three examples of predictions.

Figure 14 :Figure 15 :

 1415 Figure 14: Examples of the evolution of the electric energy with different initial solutions and different Knudsen numbers.

  we look at the influence of the Knudsen number on the result of the fluid models. All three seem to decrease in accuracy as the Knudsen number increases, but it is much more significant for the Navier-Stokes model. For Knudsen numbers above 0.1, not only the fluid model with the network seems to work systematically better than the Navier-Stokes model, but it is also

Figure 16 :

 16 Figure 16: L 2 relative errors of the fluid models over the kinetic model depending on the Knudsen number. The first plot shows the raw data and the second shows the median and interquartile interval for 20 uniform classes of Knudsen numbers between 0.01 and 1.

Figure 17 :Figure 18 :Figure 19 :

 171819 Figure 17: Example of a heat flux predicted by the neural network, with and without smoothing. The smoothing here uses σ 0.05.

Figure 20 :

 20 Figure 20: Distribution of the relative errors of the network on the test dataset with its original resolution, downsampled to 768 points and downsampled to 512 points.

Figure 21 :

 21 Figure 21: Distribution of the relative errors of the three fluid models on the kinetic model over 100 simulations with discontinuous initial solutions.

i+ 1 2 -

 2

  is the following:

	4	2 (depth)				2	1
	input X 512 (window size)	512	4				512	output Y 512
	3 levels	D e s c e n d i g n	256			256	4	= s y n t h e s i s
		p a r t				p a r t
			= a n a l y s i s	128	8	A s c e n d i n g
	Convolution Softplus	Down-sampling Up-sampling			Summation Weighted mean	Small shortcut Big shortcut

Table 1 :

 1

	-parameter	Value
	size of the input window (N )	512
	number of levels ( )	5
	depth (d)	4
	size of the kernels (p)	11
	activation function	softplus

Hyper-parameters of the reference neural network ation of the V-Net architecture, we initialize the transposed convolutions with constant kernels.

B.2 Semi-implicit scheme for the Navier-Stokes closure

With the Navier-Stokes closure q = -3 2 εp∂ x T , the first two equations remain the same as above and are solved using the same explicit finite volume method, while the third one reads

To solve this equation we use the relation w = 1 2 ρu 2 + 1 2 ρT to turn it into an equation on T , as ρ and u are known after solving the first two equations. We use a finite difference approximation for the term ∂ x (p∂ x T ), and T n+1 can then be computed by solving the following linear system:

Here the only unknown is T n+1 , and F(U) [START_REF] Besse | The multi-water-bag equations for collisionless kinetic modeling[END_REF] is the third coordinate of F(U). Finally, (w) n+1 can be computed using again the relation w = 1 2 ρu 2 + 1 2 ρT . Here the time step is chosen such as to satisfy the same CFL condition [START_REF] Dumoulin | A guide to convolution arithmetic for deep learning[END_REF] as the implicit Navier-Stokes term to not introduce additional stability constraint.


François Cretté

Antoine de Falguerolles. Looking for Growth Curves in the Situation Designed by

come

Looking for Growth Curves in the Situation Designed by

François Cretté de Palluel


Jacques Lacaze
email: jacques.lacaze@ensiacet.fr

Charbel Moussa

Yannick Thebault

Wilson L Guesser

Extrema of micro-hardness in fully pearlitic compacted graphite cast iron

Keywords: Iron alloys, casting, carbon and graphite, hardness, pearlite

come   L'archive ouverte pluridisciplinaire

Introduction

Use of compacted graphite cast iron (CGI) is strongly increasing nowadays as it shows better mechanical properties than lamellar graphite cast iron (LGI) and less micro-porosity than spheroidal graphite cast iron (SGI). CGI thus appears very well suited for further development of thin wall castings as required, in particular, in the automotive industry [START_REF] Dawson | Castability, product design, and production of compacted graphite irons[END_REF].

In order to have the highest possible mechanical properties, fully pearlitic CGIs are most often preferred. Such a matrix structure is obtained by alloying cast irons with copper, manganese and tin which are strong pearlite promoters. However, too high a level of any of these elements is known to be detrimental in one way or another, see e.g. Neumeier et al. [START_REF] Neumeier | Tin and copper combinations in as-cast ductile iron[END_REF].

Accordingly, optimizing the addition in pearlite promoter elements has long been of concern in the cases of SGI (see for example Yu et al. [START_REF] Yu | The effect of molybdenum, copper, and nickel on the microstructure, hardness, and hardenability of ductile cast irons[END_REF]) and LGI (see for example the review by Janowak and Gundlach [START_REF] Janowak | A modern approach to alloying gray iron[END_REF]), but has received only little attention in the case of CGI [START_REF] Gumienny | Effect of copper on the crystallization process, microstructure and selected properties of CGI[END_REF]. In fact, the characteristics looked for pearlite are the same for all kinds of cast irons and were stated simply by Janowak and Gunlach: "fine pearlite with uniform strength and hardness". Such an aim calls for: 1) increasing pearlite fineness, i.e. decreasing the interlamellar spacing, which increases hardness and improves strength, and: 2) achieving a uniform interlamellar spacing of pearlite as it is thought to ensure uniform matrix hardness.

On the other hand, it is well known that cast irons with high hardness and fine interlamellar spacing of pearlite give rise to machinability difficulties. In a recent work dealing with drilling test of gray irons with different alloying elements, Silva [START_REF] Silva | Furação de ferros fundidos de alta resistência aplicados em cabeçotes de motores[END_REF] showed that tool life decreases as Brinell hardness and microhardness increase and proposed that this is a result of decreasing pearlite interlamellar spacing. This suggested further investigating the homogeneity of the material at the scale of the microstructure using micro-hardness measurements. In the present study, the effect of copper and tin in fully pearlitic CGI is considered at constant manganese content.

Experimental details

A series of nine melts were prepared in an induction furnace using a base melt having the composition listed in Table 1 to which various copper and tin additions were done in the treatment ladle. Chemical analysis was performed by Optical Emission Spectroscopy for all elements but sulfur and carbon which were measured by combustion. The control of the CGI liquid metal was carried out using the Sintercast Process with a special thermal analysis procedure after the first addition of magnesium aimed at compacting graphite [START_REF] Guesser | Production experience with compacted graphite iron automotive components[END_REF]. The final correction of magnesium content in the melts to reach 0.015±0.004 wt.% and the inoculation treatment followed this thermal analysis. The melts were cast in stair castings with 7 stair thicknesses varying from 5 mm to 60 mm. The shake out time was 12 hours, meaning that all phase transformations happened with the part inside the mold. The stair 30 mm in thickness was selected for the present study and it was found that four of the nine melts showed a fully pearlitic matrix in the middle of the stair where all subsequent measurements were performed. The carbon, silicon, copper and tin contents of these four alloys are listed in Table 2. Figure 1 illustrates the microstructure of alloy B before (a) and after (b) Nital etching of the same area. For each alloy, three macro-hardness HB measurements have been performed with a load of 30 kg. On each sample, 100 microhardness HV measurements were carried out along a square grid having a 200 µm step size.

Thanks to a slight etching of the metallographic sections, the exact location of each measurement could be adjusted by having each indent on one pearlite grain. Preliminary tests suggested using a charge of 100 g that leads to indent diagonals of 20-30 µm in length, corresponding to a maximum depth at the center of 3-4 µm. This depth is large enough to assume that etching did not affect the measurements.

Results and discussion

The average macro-hardness values listed in table 2 are in the expected range of values [START_REF] Yu | The effect of molybdenum, copper, and nickel on the microstructure, hardness, and hardenability of ductile cast irons[END_REF][START_REF] Janowak | A modern approach to alloying gray iron[END_REF].

It is observed to increase slightly with Sn addition (compare alloys B and C) and with copper addition (compare alloys A and C or B and D).

Figure 2 shows an example of a micro-hardness indent which was properly located within a pearlite grain. Contrary to most observations during the present study, the interlamellar spacing was large enough to be resolved in this micrograph, certainly because of the inclination of the lamellae with respect to the metallographic section. Micro-hardness distributions appeared quite extended, which suggested representing them with percentiles as in figure 3. For each alloy, the average value is shown as a dot within a grey rectangle that contains 50% of the data while the upper and lower percentiles contain each 25% of the data.

It is seen first that the average micro-hardness does very slightly increase with copper content as expected from literature [START_REF] Gumienny | Effect of copper on the crystallization process, microstructure and selected properties of CGI[END_REF] and in agreement with the slight increase in macro-hardness.

However, the striking result is that the upper percentile enlarges dramatically when copper content increases and surprisingly when tin content decreases. The observed variability in the micro-hardness values is illustrated in Appendix A. To explain this result, it could first be argued that the variability of the micro-hardness values may be due to differences in the local internal stresses within the different pearlite grains. This possibility was checked by measuring the area of a few indents. Referring to Appendix B, it was concluded that micro-hardness variability is not caused by internal stress changes.

The mechanical response underneath a Vickers indenter is mainly defined by the material strength for a given representative strain which is strain hardening dependent [START_REF] Hernot | Study of the concept of representative strain and constraint factor introduced by Vickers indentation[END_REF].

Accordingly, the increase of the maximum micro-hardness values with copper content could be thought to be due to an effect of interlamellar spacing because a smaller interlamellar spacing relates to a decreased average distance for dislocation motion as quantified by Ray and Mondal [START_REF] Ray | The effect of interlamellar spacing on strength of pearlite in annealed eutectoid and hypoeutectoid plain carbon steels[END_REF] for pearlitic steels. Therefore the interface density within the indented pearlite grains plays a major role in increasing the mechanical behavior and more precisely strain hardening. Copper decreases the interlamellar spacing of pearlite by about 10% in cast irons as reported by Svensson and Sjögren [START_REF] Svensson | The ffect of cooling rate, section size and alloying on matrix structure formation in pearlitic grey cast iron[END_REF] and Nayyar et al. [START_REF] Nayyar | Investigation of microstructure and material properties of 18 different graphitic cast iron model materials with focus on compacted graphite iron (CGI)[END_REF] who measured values changing from 300 to 100 nm when the cooling rate was increased from 0.005 to 1.9 °C/s.

This effect of copper is small and could hardly explain the occurrence of high extrema as observed in the present work.

The most probable reason for explaining the effect of copper in cast irons appears to be the strong decrease in copper solubility with temperature in both austenite and ferrite and the associated potential precipitation of nearly pure copper (epsilon-Cu in the following). Figure 4 presents an isopleth Fe-Cu section that illustrates this decrease for a 3.5 wt.% carbon, 2.5 wt.% silicon and 0.24 wt.% manganese alloy. From this figure, such a precipitation is very much expected in ferrite of usual cast irons, even though they normally contain less than 1-1.5

wt.% copper. The associated strengthening mechanism has been characterized for Fe-Cu alloys [START_REF] Nakashima | Interaction between dislocation and copper particles in Fe-Cu alloys[END_REF] and has long found application to strengthen martensitic steels such as alloy 15-5 PH. A low temperature heat treatment of such an alloy leads to a nano-scale precipitation of epsilon-Cu in martensite between 673 K (400°C) and 773 K (500°C) as evidenced by differential thermal analysis [START_REF] Herny | Study by differential thermal analysis of reverse spinodal transformation in 15-5 PH alloy[END_REF].

Some work has been devoted to the study of this nano-precipitation of copper in steels and Fe-rich alloys by transmission electron microscopy (TEM). Chairuangsri and Edmonds [START_REF] Chairuangsri | The precipitation of copper in abnormal ferrite and pearlite in hyper-eutectoid steels[END_REF] have studied in detail the precipitation and ageing of epsilon-Cu in pearlite of alloys with about 1.5 wt.% carbon and 1.18 or 1.96 wt.% copper. Most of the precipitation occurs in the ferrite and at the interface between cementite and ferrite, but a small portion may also occur in the cementite constituent of pearlite. Murakami et al. [START_REF] Murakami | Microstructure and mechanical property of copper bearing eutectoid steel[END_REF] have studied the precipitation of epsilon-Cu during isothermal growth of pearlite in steels (0.8 wt.% carbon, 0.26 wt.% silicon and 0.35 wt.% manganese) alloyed with 0.06, 1.06 and 2.06 wt.% copper. In the case of the alloy with 2.06 wt.% Cu, they observed that epsilon-Cu precipitates both at the ferrite/austenite interface and in ferrite at 923 K (650°C), whereas it precipitates mainly in ferrite at 873 K (600°C). This precipitation results in an increase in the micro-hardness of the pearlite compared to that of the alloy with 0.06 wt.% Cu. In more recent studies, Razumakov et al. [START_REF] Razumakov | The structure and properties of cast iron alloyed with copper[END_REF] and Bataev et al. [START_REF] Bataev | Ferrite and perlite hardening in copper-alloyed steels and irons[END_REF] have confirmed the effect of epsilon-Cu on hardening of ferrite and pearlite in standard steels and other low-silicon irons. They report a large increase in micro-hardness when the copper content varies up to 3.15 wt.%, followed by little change when the copper content is higher. Recent results by Gumienny et al. [START_REF] Gumienny | Effect of copper on the crystallization process, microstructure and selected properties of CGI[END_REF] show exactly the same evolution for compacted cast irons containing 1.0 to 3.8 wt.% copper. Considering the evolution of micro-hardness during heat-treatment demonstrated by Murakami et al. [START_REF] Murakami | Microstructure and mechanical property of copper bearing eutectoid steel[END_REF], an attempt to use differential thermal analysis on the cast irons under study -as previously done for 15-5 PH alloy [START_REF] Herny | Study by differential thermal analysis of reverse spinodal transformation in 15-5 PH alloy[END_REF] -was carried out. This was not successful, probably because most of the copper precipitation was completed during cooling of the castings. It should be mentioned that the proof of copper precipitation in cast irons has already been made by Gilbert [START_REF] Gilbert | Properties of pearlitic and annealed ferritic nodular-irons alloyed with copper[END_REF] in the case of alloys containing 0.5 to 1.0 wt.% copper. This was achieved by heat treatment of the alloys to make the precipitates coarser, which however leads to loss of information on the early stages of precipitation. Fortunately, although ferrite in conventional cast irons is magnetic at room temperature because of silicon, it has recently been shown that imaging of copper precipitates by TEM is possible [START_REF] Garcia | Identification of Cu-rich precipitates in pearlitic spheroidal graphite cast irons[END_REF]. This opens up possibilities for investigating and optimizing early stages of copper precipitation.

Finally, the most unexpected and still unexplained result is the effect of tin addition which decreases the maximum pearlite micro-hardness values. In fact, when the tin content is increased, the macro-hardness increases, the average micro-hardness almost does not change and the micro-hardness maximum value decreases. This latter decrease is thus probably due to a local phenomenon. It is hypothesized that it is related to tin decreasing the undercooling for pearlite nucleation thanks to precipitation of a transient Fe 3 SnC phase at the graphite/matrix boundary [START_REF] Lacaze | Effect of tin on the phase transformation of cast irons[END_REF]. This would lead to an increase in the interlamellar spacing for given casting conditions and cooling rate and thus to a decrease in micro-hardness close to graphite precipitates. Note that this explanation would apply to fully pearlitic matrix and not to cast irons containing some ferrite located around graphite precipitates as in the recent work by Lyu et al. [START_REF] Lyu | Effect of tin on microstructure and mechanical properties of compacted graphite iron[END_REF]. This suggestion would be worth dedicated experiments involving isothermal growth of pearlite after austenitization.

Conclusion

Fully pearlitic compacted graphite cast irons alloyed with high copper and low tin contents showed maximum micro-hardness values which are significantly higher than for pearlite in alloys with low copper and tin contents. It was further observed that high values of tin decrease these maximum values. The effect of high copper content in increasing microhardness values may be convincingly related to precipitation of epsilon-Cu while the effect of high tin content is suggested to relate to coarser pearlite. In practice, precipitation of epsilon-Cu leading to hard points in the matrix may be detrimental to machining performances and the present results suggest limiting copper additions to 1 % mass in fully pearlitic cast irons.

Appendix A

It is seen in figure 3 that the minimum micro-hardness values of the pearlite are nearly the same for all alloys. A micrograph of most of the corresponding largest indents was taken as illustrated in figure A1 for alloys B and C. In some cases, it was seen that the indent may have been affected by a nearby graphite particle as in figure A1-b, thus explaining the lower hardness value. However, this does not seem to be verified in most cases, e.g. the indent in figure A1-a where the large interlamellar spacing seen on this section plane could instead be the explanation.

Figure A2 presents similarly indents for highest micro-hardness values for the same alloys as in figure A1. While most of these small indents were made in grains containing pearlite too fine to be resolved by optical microscopy (e.g. figure A2-b), this is not the case of the one in figure A2-a. Thus, it seems there is no one to one relation between the lamellar spacing seen on the section plane and the micro-hardness value, and this has been seen described quantitatively elsewhere [START_REF] Freulon | Solidification and room temperature microstructure of a fully pearlitic compacted graphite cast iron[END_REF].

Figure 1 -

 1 Figure 1 -Microstructure of alloy B before (a) and after (b) Nital etching. These micrographs represent the same area and the scale is the same for both.

Figure 2 -

 2 Figure 2 -Example of a micro-hardness indent within a pearlite grain where interlamellar spacing can be resolved (Alloy B).

Figure 3 -

 3 Figure 3 -Percentiles of pearlite micro-hardness distributions as a function of copper content.

Figure 4 -

 4 Figure 4 -Isopleth Fe-Cu section for an alloy containing 3.5 wt.% carbon, 2.4 wt.% silicon and 0.24 wt.% manganese. Calculation made with the THERMOCALC software and the TCFE-8 database. FCC_A1 and BCC_A2 stand respectively for austenite and ferrite.

Figure

  Figure A1: micro-hardness prints with low micro-hardness value for alloys B (a) and C (b).

Figure

  Figure A2: micro-hardness prints with high micro-hardness value for alloys B (a) and C (b).

Figure

  Figure A3: micro-hardness prints in two pearlite grains showing close lamellar spacing (alloy B).

Figure A4 :

 A4 Figure A4: SEM micrograph of each of the micro-hardness prints seen in figure A3.

  To verify if the variability observed in the micro-hardness values during the present study could be due to internal stresses variability, micro-hardness was calculated using the diagonals and the projected area. Two indents giving different hardness values were analyzed.The micro-hardness was evaluated using diagonals, on the one hand, and the projected area, on the other hand, see figureB2. Both calculation methods gave similar results which indicate that the micro-hardness variability is not due to local internal stresses.

Figure B2 .

 B2 Figure B2. Two micro-hardness indents with HV=352 (a) and HV=293 (b). Solid lines indicate the diagonals and delimit the projected area used for hardness calculations..

Table 1 -

 1 Composition of the base melt (mass%).

	C	Si	Mn	P	S	Cr	Ti	B	Sn	Cu
	3.52	2.4	0.24 0.026 0.009 0.035 0.01 0.0001 0.018 0.183
	Pb	Bi	Al	Mo	Sb	Ni	Nb	Zn	V	As
	0.0003 0.0004 0.003 0.002 0.001 0.001 0.003 0.044 0.002 0.001

Table 2 -

 2 Carbon, silicon, copper and tin contents (mass %) measured on the four selected

		samples, and macro-hardness values
	Alloy reference	C	Si	Cu	Sn	HB
	A	3.57	2.26	0.58	0.094	257
	B	3.53	2.29	0.96	0.057	246
	C	3.53	2.25	0.99	0.095	266
	D	3.48	2.30	1.28	0.056	278

Acknowledgment

Thanks are due to David de Jonghe for performing the micro-hardness measurements during a one month internship at CIRIMAT laboratory. This research did not receive any specific grant from funding agency in the public, commercial, or not-for-profit sectors.

Appendix B

The influence of internal stresses on the measurement of Berkovich hardness was studied by Tsui et al. [START_REF] Tsui | Influences of stress on the measurement of mechanical properties using nanoindentation: Part I. Experimental studies in an aluminum alloy[END_REF]. In fact, internal stresses influence displacement depth of the indenter within the material and also the shape of the indent with more or less pile-up or sink-in, whatever the indenter shape. Compressive internal stresses tend to increase the pile-up or decrease the sinkin of the imprint while tension internal stresses act in the opposite way. The influence of internal stresses on the indenter displacement is thus compensated by the pile-up or sink-in of the imprint while the contact surface (or the projected area) between the indenter and the material is not influenced. However, Vickers hardness is generally calculated using the diagonals length which relies on the assumption that the indent presents no pile-up or sink-in.

A schematic is presented in Figure B1 to illustrate the errors on hardness measurement that can be made due to internal stresses. When no internal stresses are present, figure B1-a, there is no pile-up or sink-in. If compressive or tension internal stresses are present the imprint presents pile-up (figure B1-b) or sink-in (figure B1-c), respectively. The 3 cases of figure B1 have different diagonals, therefore they will have different hardness values when evaluated using the diagonals. However, the projected areas defined with the dashed lines are equal for the 3 cases. Therefore, using the projected area to calculate the hardness will lead to the same value for the 3 cases, which is correct since internal stresses should not influence the hardness value of the material as this is an intrinsic property.