Vincent Labarre
email: vincent.labarre@lsce.ipsl.fr

Bérengère Dubrulle
email: berengere.dubrulle@cea.fr

Didier Paillard

Sub-grid modelling for a diffusive lattice gas

Keywords: Lattice Gas, Statistical Physics, Sub-Grid Modelling

come

Introduction

Turbulent flows are ubiquitous in natural and industrial flows. Despite their large variety, they all share common characteristics, such as chaotic dynamics (i.e. sensitivity to initial condition) and huge range of scales. The combination of the two properties forbids the prediction of their evolution in a deterministic way, and claims for a statistical description [START_REF] Stephen B Pope | Turbulent flows[END_REF]. On the other hand, the smallest scales of the flow are irrelevant for practical purposes, while their resolution represents the most important part of the computational cost. This justifies procedures where only the evolution of the coarsegrained fields is considered, while sub-grid scales are modelled. There is still no universal approach to represent the effect of unresolved scales on the resolved ones. The problem couples the dynamics and the coarse-graining procedure. All methods need, at some point, to introduce empirical parameters (eddy viscosity models [START_REF] Lesieur | New trends in large-eddy simulations of turbulence[END_REF]) and/or idealized physical hypotheses (see for example [START_REF] Chavanis | A parametrization of two-dimensional turbulence based on a maximum entropy production principle with a local conservation of energy[END_REF]). This may lead to unwanted drawbacks, like non-conservation of energy or uncontrolled dependence of the empirical parameters with the coarse graining (in the Smagorinski model for example [START_REF] Smagorinsky | General circulation experiments with the primitive equations[END_REF], the size of the mesh in the numerical grid is introduced for dimensional consistency).

Arbitrariness of the modelling could be decreased using systematic and well controlled procedure, that includes coarse-graining. One major difficulty relies on the nonequilibrium nature of the turbulence, in which all physical quantities (mass, momentum, energy, tracers) flow into a system due to an external forcing. While there is presently no systematic out-of-equilibrium theory for turbulence, there has been huge progress in characterization of out-of-equilibrium properties of simplified statistical models. Those models fall into two categories. First, Lattice Gas Cellular Automata (LGCA) are a discrete space-time models version of kinetic theory of gas [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Grosfils | Spontaneous fluctuation correlations in thermal lattice-gas automata[END_REF][START_REF] Dieter | Lattice-gas cellular automata and lattice boltzmann models -an introduction[END_REF][START_REF] Rivet | Lattice Gas Hydrodynamics[END_REF]. It was shown, by using classical tools of kinetic theory, that they are able to represent real flow under reasonable approximations. Second, Stochastic Lattice Gases (SLG) have been introduced in several contexts (DNA modelling, traffic flows, paradigmatic non-equilibrium systems) [START_REF] Mallick | The exclusion process: A paradigm for non-equilibrium behaviour[END_REF]. They have the particularity to be conceptually simple enough to allow detailed computations and complicated enough to represent some interesting features like phase transitions. For example, the detailed statistics of density and current has been derived for some one-dimensional models using the Matrix Anzatz [START_REF] Derrida | An exact solution of a one-dimensional asymmetric exclusion model with open boundaries[END_REF], or the Bethe Anzatz [START_REF] Lazarescu | The physicist's companion to current fluctuations: one-dimensional bulkdriven lattice gases[END_REF]. It was also shown that we can compute the probability of observing atypical density and current profiles directly at the macroscopic level (i.e. computing the large deviation functional of the density and current) in the stationary state using the Macroscopic Fluctuation Theory [START_REF] Bertini | Macroscopic fluctuation theory for stationary non-equilibrium states[END_REF][START_REF] Derrida | Non-equilibrium steady states: fluctuations and large deviations of the density and of the current[END_REF][START_REF] Derrida | Microscopic versus macroscopic approaches to non-equilibrium systems[END_REF][START_REF] Bertini | Macroscopic fluctuation theory[END_REF]. Both LGCA and SLG represent ideal bench test for conceptual ideas on sub-grid modelling because they are conceptually far more simple than a real fluid and are easily simulated. The main difference between LGCA and SLG is that LGCA are not solvable analytically (except in the lattice Boltzmann approximation in some case [START_REF] Zou | Analytical solutions of the lattice boltzmann bgk model[END_REF]).

The goal of the present study is to look for an improved understanding of basic subgrid modelling, by using a simplified model of diffusion based upon a two-dimensional diffusive LGCA. In classical turbulence theory, such systems would be naturally modelled using an eddy viscosity. Here, we show that the model is actually more involved, while still conceptually simple. Our procedure relies upon numerical analysis of the typical behaviour of the coarse-grained current in a two-dimensional diffusive LGCA as a function of the coarse-grained density, its gradient, and the size of the cubic spatiotemporal coarse-graining window. We show that the temporal variation of the coarsegrained current, at fixed coarse-grained density and its gradient, is a Gaussian random variable whose average is given by a relaxation to the hydrodynamic current at a given rate and root mean square that are strongly influenced by the coarse-graining. Using these observations, we propose a simple sub-grid model to describe the evolution of the coarse-grained current.

The outline of the paper is as follows. In a first part, we present our model (part 2). In a second part, we present our numerical results for the conditional PDF of the temporal variation of the coarse-grained current (part 3) and we propose a simple stochastic sub-grid model for the coarse-grained current. We discuss our results and perspectives in part 4. The derivation of the evolution equation of the average density is given in the Appendix A.

Model

Microscopic dynamics

We study a variant of the model introduced by Hardy, de Pazzis and Pomeau [START_REF] Hardy | Time evolution of a two-dimensional classical lattice system[END_REF][START_REF] Hardy | Molecular dynamics of a classical lattice gas: Transport properties and time correlation functions[END_REF] (HPP model). It consists on a square lattice L of L × L nodes (see figure 1). Each node consists in four channels which correspond to the velocities of the model (c 1 = (1, 0), c 2 = (0, 1), c 3 = (-1, 0), c 4 = (0, -1)). Occupation of the channel i at discrete time t * and discrete position r * = (x * , y * ) is noted n i (r * , t * ). A channel can be occupied by one particle of unit mass or can be empty (exclusion principle), such that n(t * ) = (n i (r * , t * ), i = 1, ..., 4, r * ∈ L) is the Boolean field that gives the configuration of the system at time t * . At each discrete time step, particles at a same node are collided and propagated to nearest nodes according to their velocities. Collision changes a precollision state s = (s i , i = 1, ..., 4) to a post-collision state s = (s i , i = 1, ..., 4) with probability A s s where s i is the occupation of channel i (Boolean variable). Here, we consider two types of collisions. Frontal collision occurs with a probability p when there is two particles with opposite velocities at the node. This collision conserves mass and momentum. Three particles collision occurs with a probability q. It consists to flip the velocity of the particle which does not face another one. This collision conserves mass but not momentum. It induces a tendency of particles to turn round more frequently when going in high density regions (where three particles collisions occurs more frequently). This collision destroys spurious invariants like the particles-holes symmetry (duality) and obviously the conservation of horizontal (vertical) current along horizontal (vertical) lines. The three particle collision is the main difference between our model and the HPP model. These collision rules respect the detailed balance (A s s = A ss ∀s, s ), and hence semi-detailed balance:

s A ss = s A s s = 1, ∀s, (1)
Since only the number of particles is conserved by collisions, we observe a diffusive behaviour at large scales. Propagation consists to move particles to neighbouring nodes according to their velocity n i (r * ) → n i (r * + c i ). If p = q = 1 then the evolution of the system inside the bulk is deterministic. At the boundaries, particles are injected in vacant channels (i.e. those that are not filled by particles propagated from the bulk) with a probability equal to the density of the reservoir divided by four (the number of channels per node). We introduce the microscopic observables

ρ * (r * , t * ) = 4 i=1 n i (r * , t * ) = 0, 1, 2, 3, 4, (2)
j * α (r * , t * ) = 4 i=1 n i (r * , t * ) c iα = -1, 0, 1, (3)
which are the microscopic density and the microscopic current of particles in the α's direction.

To obtain the macroscopic evolution of Lattice gas [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Dieter | Lattice-gas cellular automata and lattice boltzmann models -an introduction[END_REF][START_REF] Rivet | Lattice Gas Hydrodynamics[END_REF], we introduce an ensemble average procedure • that connects the discrete, microscopic state with the macrostate which is assumed to vary slowly (at least in the limit of large lattice size and large time). We introduce

N i = n i , ρ = ρ * , j α = j * α , (4)
that are the average population of particle with speed c i , the average density, and the α's component of the average current of particles. We introduce the following macroscopic space and time variables

x = x * , y = y * , t = t * (5)
where = 1/L. This convective scaling is more relevant to study the coarse-grained dynamics than the diffusive scaling used in the appendix. It can be shown (see Appendix A), that we have the following evolution equation for the average density:

∂ t ρ + ∂ α j α = 0, (6)
with average current (changed by a factor compared to the appendix due to the convective scaling).

j α (ρ, ∂ α ρ) = 1 4 - 4 qρ 2 ∂ α ρ. (7)
We note the presence of an anti-diffusive term 4 ∂ α ρ in the expression of the current. This is a spurious term that always appears in the transports coefficients of lattice gas models due to the discrete nature of the dynamics (see, for example, equation 7.4 of [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF].). This term is not problematic here since the lattice gas is viewed as a dynamical system that allows illustration of the sub-grid modeling problem while keeping numerical and theoretical analysis simple. The important relevant characteristic of our system is its non-equilibrium nature. We show in figure 2 the microscopic and average density field in stationary state for our simulation: L = 300, p = q = 1, ρ L = ρ B = 2.4 and ρ R = ρ T = 1.6 and 10 6 time steps in the statistically stationary state. It is also possible to show that the variance of the microscopic current depends on the average density:

(j * α -j * α ) 2 = ρ 2 1 - ρ 4 ≡ σ 2 * (ρ). (8)

Coarse-grained dynamics

We coarse grain space and time in cubic mesocells of τ time steps and τ × τ nodes where τ is called the coarse graining factor. It is therefore possible to define the coarse-grained density ρ τ and the coarse-grained current j τ as the average values of the microscopic density and the microscopic current over the mesocell

M τ (r, t) = {(x * , y * , t * ) | 0 ≤ x * -Lx < τ, 0 ≤ y * -Ly < τ, 0 ≤ t * -Lt < τ } at position r = (x, y)
and time t:

τ 3 ρ τ (r, t) = (r * ,t * )∈Mτ (r,t) 4 i=1 n i (r * , t * ) ( 9
)
τ 3 j τ α (r, t) = (r * ,t * )∈Mτ (r,t) 4 i=1 n i (r * , t * ) c iα . (10)
An instantaneous coarse-grained density field for our simulation with τ = 5 is shown in figure 3. Since there is no net separation of scales (i.e. τ cannot be considered infinite), the fluctuations of the coarse-grained density field are still important. To simulate the evolution of the coarse-grained field directly, we need a relation that gives the coarse-grained current as a function of the coarse-grained density, its gradient and the coarse graining factor to close the continuity equation

∂ t ρ τ + ∂ α j τ α = 0. ( 11
)
Since the collisions are local (i.e. involve only particles at one node at each time step), the information does not travel faster than the particles which propagate at a fixed speed 1. It follows that what happen in a mesocell at position (x, y) and a given time t only depends on the configurations of the neighbouring mesocells (x ± 0, τ, y ± 0, τ ) at times t -τ . It means that the law governing the evolution of the coarse-grained density is local. Because there is no scale separation, we expect to observe a lag between the coarse-grained current response and the applied coarse-grained density (this rheological property is shared by turbulence and viscoelastic fluids). We therefore have to model the temporal variations of the coarse-grained current. In general, these variations are expected to vary with the coarse-grained density, its gradient, the coarse graining factor, but also the dynamics (i.e. p and q). Assuming that the temporal variation of one component of the current does not depend on the coarse-grained gradient in the transverse direction, we search a constitutive relation for the current as:

∂ t j τ α = F τ α (ρ τ , ∂ α ρ τ , p, q, η) ( 12
)
where η is the noise representing the sub-mesocell variability.

To this end, we save for each mesocell of the simulation: the coarse-grained current at two successive times and the coarse-grained density measured in the mesocell, and the coarse-grained density of the eight neighbouring mesocells. It allows to compute the quantities

v τ α (r, t) ≡ j τ α (r, t + τ ) -j τ α (r, t) τ α = x, y, (13)
g τ x (r, t) ≡ ρ τ (x + τ, y, t) -ρ τ (x -τ, y, t) 2τ , (14)
g τ y (r, t) ≡ ρ τ (x, y + τ, t) -ρ τ (x, y -τ, t) 2τ , (15)
which are respectively the first-order approximation of the temporal derivative of the coarse-grained current and the second-order approximation of the spatial derivatives of the coarse-grained density in the microscopic unit. It follows that v τ α / and g τ α / are the finite differences approximations of the spatial derivative in the macroscopic scale estimated with the coarse-grained fields. These are the quantities that need to be modeled in order to perform simulation of the coarse-grained field directly.

Results

Conditional probability distribution of the current

We define p τ α (v| ĵ, ρ, ĝ) the conditional PDFs of observing

v τ α ∈ [v -dv 2 , v + dv 2 ] knowing that j τ α ∈ [ ĵ -d ĵ 2 , ĵ + d ĵ 2 ], ρ τ ∈ [ρ -dρ 2 , ρ + dρ 2 ]
, and

g τ α ∈ [ĝ -dĝ 2 , ĝ + dĝ 2 ]
. In order to identify the structure of these PDFs, it is quite useful to work with normalized random variables. Therefore, we introduce the following random variables:

δ τ α ( ĵ, ρ, ĝ) = v τ α -v τ α ĵ,ρ,ĝ (v τ α -v τ α ĵ,ρ,ĝ ) 2 ĵ,ρ,ĝ α = x, y, τ > 1 ( 16
)
where

v n τ α ĵ,ρ,ĝ = vn p τ α (v| ĵ, ρ, ĝ) dv (17)
is the conditional n th moment of the variation of the coarse-grained current knowing the coarse-grained density and its gradient.

In figure 4, we have represented the PDFs of δ τ α ( ĵ, ρ, ĝ) for various values of ρ, ĝ, τ = 10, 15, 20, 20, 30 and α = x, y. We emphasize that the analysis relies on the numerical simulation in the statistically stationary steady-state for ρ L = ρ B = 2.4, ρ R = ρ T = 1.6, and p = q = 1. The histograms are computed dividing ĵ into 19 windows between -0.04 and 0.04, ρ divided into 16 windows between 1.2 and 2.8, and ĝ divided into 19 windows between -0.03 and 0.03. The distributions are colored differently depending on ĵ. We have only retained the samples of size > 2000 to compute the empirical PDFs in order to limit noise. For small values of the current | ĵ| 1 (in blue in figure 4), one observes a universal Normal behavior meaning that v τ α is a Gaussian random variable fully determined by its mean value and its root mean square (that a priori both depend on the coarse-grained density, its gradient, and the coarse-graining factor) in the limit of low current (or equivalently the low mean speed of particles). We also observe that higher values of the coarse-grained current (in red in figure 4) are associated with a non-Gaussian behavior. This non-Gaussian behavior for high-values of the current is of secondary importance for two reasons. The first is that lattice gas is usually studied in the low-speed limit, in particular, to ensure the isotropy of the macroscopic dynamics. The second reason is that the high values of the current are atypical. When we compute the empirical conditional probability distributions of observing v τ α only knowing ρ τ and g τ α (but not j τ α ), we recover a Gaussian behavior. Then, it is of prime importance to model the coarse-grained dynamics in the low current limit.

The Gaussian behavior leads us to assume that v τ α can be modelled by the following equation

v τ α = - j τ α -µ τ α r τ + σ τ η α , (18)
where µ τ α is the relaxed current (i.e. in stationary state in the absence of fluctuations), r τ is the relaxation time, σ τ is the root mean square of v τ α , and η α α = x, y are decorrelated normal random variables. Both µ τ α , r τ and σ τ have to be modelled to construct the sub-grid model. If they are considered as constant, the Gaussian behavior will follow. Note however that if we allow σ τ and r τ to be statistical variables, we could get a non-Gaussian solution of equation ( 18) as a superstatistics [19]. In the sequel, however, we concentrate on low-current modeling and assume that σ τ and r τ are deterministic functions that need to be modeled. For a given value of ρ τ and g τ α , the linear regression for the points v τ α vs j τ α has a slope -1/r τ , an intercept µ τ α /r τ and an error term σ τ . This is how we estimate µ τ α , r τ and σ τ numerically. In the following, we discuss the behaviour of these functions with ρ τ , g τ α , q and τ .

Relaxed current µ τ α

The figure 5 shows µ τ α for τ = 10 in our simulation and the average current predicted in the hydrodynamic limit with a convective scaling (given by equation ( 7)). We do not show the results for other values of τ because our statistics are not sufficient to provide a converged estimate of µ τ α for large values of g τ α . Quite remarkably, we observe that the coarse-grained current tends to relax to the average current (obtained in the hydrodynamic limit) estimated by the local coarse-grained density and its gradient. It means that the relaxed current is invariant under the effect of coarse graining. This suggests that the relaxed coarse-grained current obeys the same functional relation as a function of the (coarse-grained) gradient density than the microscopic gradient, a property that can be viewed as a kind of Germano identity [START_REF] Germano | A dynamic subgridscale eddy viscosity model[END_REF]. In turbulence, such identity is used via a hypothesis of scale invariance of the turbulence in the inertial range, resulting in a Reynolds stress that keeps the same shape above and below the cut-off. In our case, the Germano identity amounts to take

µ τ α (ρ τ , g τ α , q) = 1 4 - 4 qρ 2 τ g τ α (19)
for our sub-grid model.

Relaxation time r τ and typical fluctuations σ τ

We show on figure 6 the relaxation time (divided by the coarse graining factor) and the root mean square of v τ α as a function of the coarse-grained density and its gradient for τ = 10. We observe that both r τ and σ τ depend little on the gradient of the coarse-grained density (at least for the range of parameter we have investigated), suggesting that the relaxation and the fluctuations of v τ α can be modelled by the equilibrium fluctuations at the coarse-grained size at a first approximation. We propose a simple model where we simplify the dynamics by considering that at each discrete time step, the mesocell Table 1. Our model for the gain/loss X = -1, 0, 1 of particle with speed c 1 at a given time at each node in the mesocell for the population N s1 . Source/Sink nodes probability distribution entry (left boundary) τ

   P (X = -1) = 0, P (X = 0) = 1 -ρ τ /4, P (X = 1) = ρ τ /4. exit (right boundary) τ    P (X = -1) = N s1 , P (X = 0) = 1 -N s1 , P (X = 1) = 0. 2 particules collision τ 2    P (X = -1) = N s1 (1 -N s2 )N s3 (1 -N s4 )p, P (X = 0) = 1 -P (X = -1) -P (X = 1), P (X = 1) = (1 -N s1 )N s2 (1 -N s3 )N s4 p. 3 particules collision τ 2    P (X = -1) = N s1 N s2 (1 -N s3 )N s4 q, P (X = 0) = 1 -P (X = -1) -P (X = 1), P (X = 1) = (1 -N s1 )N s2 N s3 N s4 q.
is spatially statistically homogeneous (i.e. we neglect the spatial inhomogeneities at all time step of the mesocell). This crude approximation allows to describe the state of the box at each time step by the mean occupations of channels in the mesocell, noted N si (t * ), i = 1, ..., 4. It implicitly assumes an instantaneous mixing of particles at each discrete time step inside the bulk of the mesocell. To write the evolution equation for the N si 's, one must distinguish the surface effects (rate at which particles enter/leave the box), and the volume effect due to collisions. Naturally, the scalings of these processes are different since they do not involve the same number of nodes. We model the variations of the populations at each time step by considering the different sources and sinks as stochastic processes. For a given population i, particles enter at a boundary with a rate ρ τ /4 at each node. They exit the mesocells at the opposite boundary with a rate N si at each node. By doing so, we implicitly neglect the collision step at boundaries. Particles populations are mixed by collisions that occur in the all bulk. It is taken into account by considering the effect of collisions of the mean population at one node. For example, the sources/sinks for N s1 are given in table 1.

Considering τ to be sufficiently large, we can approximate the random change of populations with Gaussian random variables using the central limit theorem. For example, the resulting stochastic equation describing the evolution of N s1 is

τ 2 (N s1 (t * + 1) -N s1 (t * )) = τ ρ τ 4 + τ ρ τ 4 1 - ρ τ 4 η in 1 -τ N s1 -τ N s1 (1 -N s1 ) η out 1 +τ 2 pN p + τ 2 σ 2 p η p + τ 2 qN qx + τ 2 σ 2 qx η qx ( 20
)
where η in 1 , η out 1 , η p and η qx are all independents, identically distributed normal random variables, and

N p = (1 -N s1 )N s2 (1 -N s3 )N s4 -N s1 (1 -N s2 )N s3 (1 -N s4 ), (21)
σ 2 p = (1 -N s1 )N s2 (1 -N s3 )N s4 p + N s1 (1 -N s2 )N s3 (1 -N s4 )p -(pN p ) 2 , (22)
N qx = (1 -N s1 )N s2 N s3 N s4 -N s1 N s2 (1 -N s3 )N s4 , (23)
σ 2 qx = (1 -N s1 )N s2 N s3 N s4 q + N s1 N s2 (1 -N s3 )N s4 q -(qN qx ) 2 . ( 24
)
Similarly, one can write down the equation for N s3

τ 2 (N s3 (t * + 1) -N s3 (t * )) = τ ρ τ 4 + τ ρ τ 4 1 - ρ τ 4 η in 3 -τ N s3 -τ N s3 (1 -N s3 ) η out 3 +τ 2 pN p + τ 2 σ 2 p η p -τ 2 qN qx -τ 2 σ 2 qx η qx . ( 25
)
We point out that η in 3 and η out 3 are also identically distributed normal random variables, and that η p and η qx are the same random variables involved in [START_REF] Germano | A dynamic subgridscale eddy viscosity model[END_REF] due to the conservation of the number of particles in the mesocell during the collision step. Now, we introduce the spatially coarse-grained current (over a square mesocell) in the x direction j sx ≡ N s1 -N s3 . One obtains the evolution equation for j sx by taking the difference of ( 20) and ( 25), divided by τ 2 :

j sx (t * + 1) -j sx (t * ) = - 1 τ j sx + 2qN qx + 1 τ 3 ρ τ 4 1 - ρ τ 4 η in 1 -η in 3 - 1 τ 3 N s1 (1 -N s1 ) η out 1 + 1 τ 3 N s3 (1 -N s3 ) η out 3 + 2 1 τ 2 σ 2 qx η qx . ( 26
)
To simplify the last equation, we linearise N qx around the mean occupation of channels ρ τ /4: N si (t * ) ρ τ /4 + δN si (t * ), i = 1, ..., 4 with δN si (t * ) << 1. Moreover, we consider only the additive part of the noises. We therefore evaluate

N qx ρ 2 τ 16 (δN s3 -δN s1 ) = - ρ 2 τ 16 j sx , (27)
N si (1 -N si ) ρ τ 4 1 - ρ τ 4 , i = 1, ..., 4 (28)
σ 2 qx 2q ρ τ 4 3 1 - ρ τ 4 . (29)
It follows that

j sx (t * + 1) -j sx (t * ) = - 1 τ + q 8 ρ 2 τ j sx + 1 τ 3 ρ τ 4 1 - ρ τ 4 η in 1 -η out 1 -η in 3 + η out 3 + 1 τ 2 q 8 ρ 3 τ 1 - ρ τ 4 η qx . (30)
Now, we use basic properties of independent Gaussian random variables to write this equation in a more compact form:

j sx (t * + 1) -j sx (t * ) = - 1 τ + q 8 ρ 2 τ j sx + 1 τ 2 ρ τ 1 - ρ τ 4 1 τ + q 8 ρ 2 τ η x (31)
where η x is a normal random variable. The equation ( 31) is a Langevin equation that arise in many contexts in Physics (see for example [START_REF] Gardiner | Stochastic Methods: A Handbook for the Natural and Social Sciences[END_REF] for an introduction to stochastic methods). The large deviation function of the equivalent of the time integrated current in our model, given by the limit for τ → ∞ of

j τ x (t) = 1 τ t≤ t * <t+ τ j sx (t * ), (32)
has already been computed for this process (see, for example, the last section of [START_REF] Touchette | Introduction to dynamical large deviations of markov processes[END_REF]).

One can deduce from these results that the distribution of j τ x satisfies a large deviation principle

lim τ →∞ - 1 τ ln P (j τ x = j) = I(j) (33)
where the rate function is given by

I(j) = 1 2 ( 1 τ + q 8 ρ 2 τ ) 2 1 τ 2 ρ τ (1 -ρτ 4 )( 1 τ + q 8 ρ 2 τ ) j 2 = τ 2 2 ( 1 τ + q 8 ρ 2 τ ) ρ τ (1 -ρτ 4 ) j 2 . (34)
Therefore, our model states that j τ α is a Gaussian random variable with mean zero and root mean square

ρτ (1-ρτ 4 ) τ 3 ( 1 τ + q 8 ρ 2 τ ) . We suggest that v τ α fluctuates like j τ α /τ . It is equivalent to choose σ τ (ρ τ , q) = ρ τ (1 -ρτ 4 ) τ 5 ( 1 τ + q 8 ρ 2 τ ) (35)
for our sub-grid model. In figure 7, we compare this model to the numerical results. The numerical results show that as far as we coarse grain, we change the shape and the scaling of the fluctuations of the coarse-grained current as a function of the coarse-grained density. Strikingly, σ τ is no more symmetric around the density ρ τ = 2. It goes from ρ 2 (1 -ρ 4 ) for the microscopic fluctuations (corresponding to τ = 1) to a limit function that we have not attained with our limited simulation. Our model provides the correct variation of σ τ (ρ) with respect to (ρ 2). To quantify the relative error, we compute the ratio between the numerical value of σ τ and our theoretical prediction noted R τ . It is provided in figure 7. We observe that the variations of R τ with ρ τ are rather mild, and that the deviations between the simulation result and the model remain within about 10% for sufficiently large τ (in the regime of 1.5 < ρ r < 2.0). This suggests that observed deviations from our theoretical modeling are finite-size effects.

Averaging the equation (31), the evolution of j sx in the mesocell can be written

j sx (t * + 1) = a j sx (t * ), a = 1 - 1 τ + q 8 ρ 2 τ . ( 36
)
It yields to the relation

j sx (t * ) = j sx (t * )a t * -t * ∀t * , t * ∈ N. ( 37
)
Using the last relation in the average of equation ( 32), we show that

j τ x (t + τ ) = a τ j τ x (t). ( 38
)
It follows that

j τ x (t + τ ) -j τ x (t) τ = - 1 -a τ τ j τ x (t). (39)
This equation suggests that the relaxation time for the coarse-grained current is

r τ (ρ τ , q) = τ 1 -a τ = τ 1 -1 -1 τ + q 8 ρ 2 τ τ . (40)
Since -1 -1 τ = a(ρ τ = 4) ≤ a ≤ (ρ τ = 0) = 1 -1 τ , we expect that in normal conditions (ρ τ not too close from 4) |a| < 1 so a τ → 0 when τ → ∞. We observe on figure 8 that this model correctly predicts that the relaxation time converges to τ as τ increases. The agreement between the prediction and the model is better than 10 percent.

Discussion and conclusion

In this study, we have addressed the question of sub-grid modelling in a simple diffusive Lattice Gas. In the low current limit, we observed that the first order approximation of the temporal derivative of the coarse-grained current, knowing the coarse-grained density and its gradient, is a Gaussian random variable whose mean and root mean square depend only on the size of the mesocell, the coarse-grained density, its gradient, and eventually on the dynamical parameters of the model. We have proposed a simple model where the temporal variation of the coarse-grained current is modelled as the sum of an average relaxation to a relaxed current at a given relaxation time, and a Gaussian noise. Observing the fact that the root mean square and the relaxation time depend little on the gradient of density, we have constructed a simple equilibrium model where the spatial inhomogeneities inside a mesocell are neglected to determine the variations of the root mean square of the grained current and the relaxation time. Our model agrees quantitatively with the numerical results, providing a agreement better that 20 percent for the fluctuations, and better than 10 percent for the relaxation time. In particular, it explains the variations of the relaxation time and the root mean square as a function of the coarse-grained density. The relaxed current is simply given by taking the hydrodynamic limit.

LGCA and SLG have been studied extensively as simplified model of kinetic theory [START_REF] Rivet | Lattice Gas Hydrodynamics[END_REF] or paradigmatic models of non-equilibrium physics [START_REF] Mallick | The exclusion process: A paradigm for non-equilibrium behaviour[END_REF]. However, up to our knowledge, these models are not used as conceptual object to learn about sub-grid modelling (we can note one attempt to compute a turbulent viscosity with such model [START_REF] Pierre Boon | A lattice gas automaton approach to 'turbulent diffusion[END_REF]). Yet, deriving sub-grid models for Lattice Gases seems to be an interesting problem to deal with before addressing more difficult problems like the turbulence modelling. In particular to understand the non-trivial variations of the coefficients that are introduced empirically in turbulence models. As an example, our simple model for the fluctuations of the current shows that the shape of the root mean square of the coarse-grained current (as a function of the density) is changed when we coarse grain.

Our model for the evolution of the coarse-grained density is finally

∂ t ρ τ + ∂ α j τ α = 0, ( 41
) ∂ t j τ α = - j τ α -µ τ α r τ + σ τ η α , (42)
where = 1/L is the ratio between the microscopic scale and the size of the system, η α is a Gaussian white noise and µ τ α , r τ and σ τ are specified in equations ( 19), ( 40) and (35).

For this diffusive lattice gas, there remains the question of the non-Gaussian fluctuations observed for the atypical large values of the coarse-grained current. We expect that a refined sub-grid model (that takes into account the spatial inhomogeneities inside the mesocell) can be constructed by extending recent works on SLG (like for the Zero Range Process [START_REF] Evans | Nonequilibrium statistical mechanics of the zero-range process and related models[END_REF]) to compute the statistics of coarse-grained observables, taking into account fluctuations at boundaries of the mesocell. Alternatively, one could use the superstatistics formalism to derive more elaborate sub-grid models [19,[START_REF] Beck | Superstatistics: theory and applications[END_REF][START_REF] Beck | Generalized statistical mechanics for superstatistical systems[END_REF]. It would also be interesting to perform simulations with different boundary conditions, for different values of the dynamical parameters p and q, and for transient simulation to test the robustness of the sub-grid model to various conditions. An important question is whether one can employ such a method to other systems of practical interest, like turbulent flows. For example, does a relaxation equation similar to the one used here could be used to model the Reynolds stress tensor using the analysis of the coarse-grained velocity field? To our knowledge, such Eulerian model has not been tested for other systems.

(ii) The Hydrodynamic limit consists to introduce the (dimensionless) variables

t = t * /T, x = x * /L, y = y * /L, (A.9)
where T → ∞ and L → ∞ are the macroscopic time and space (typically the time of the simulation and the size of the Lattice), such that t, r = (x, y) are considered as continuous variables. We then introduce a macroscopic population that varies continuously in space and time

f i (r, t) = n i (r * , t * ) = N i (r * , t * ). (A.10)
We keep the same notation for the continuous density and current. In the following, we will consider the diffusive scaling T = L 2 and introduce the Knudsen number = 1/L → 0.

(iii) We can show that the equilibrium solution (i.e. the solution in stationary state without perturbation) of the system is f (0) i = ρ/4 (we can invoke the Universality Theorem for lattice gases, valid in case of semi-detailed balance and no spurious invariant, in the simple case where only mass is conserved by collisions [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Rivet | Lattice Gas Hydrodynamics[END_REF]). We assume small perturbations around this equilibrium such that

f i = ρ 4 + f i , 4 i=1 f i = 0. (A.11)
Averaging (A.1-A.4), using (A.8), substituting (A.11) in (A.5-A.7) and replacing the finite difference by derivatives up to second order in with the diffusive scaling, we obtain

∂ x + 2 ∂ t + 2 1 2 ∂ 2 x ρ 4 + f 1 = p ρ 4 1 - ρ 4 (f 2 + f 4 -f 1 -f 3 ) (A.12) + 2 p 1 - ρ 2 (f 2 f 4 -f 1 f 3 ) + q (f 3 -f 1 ) ρ 2 16 + 2 q ρ 4 (f 3 -f 1 ) (f 2 + f 4 ) + O( 3 ), -∂ x + 2 ∂ t + 2 1 2 ∂ 2 x ρ 4 + f 3 = p ρ 4 1 - ρ 4 (f 2 + f 4 -f 1 -f 3 ) (A.13) + 2 p 1 - ρ 2 (f 2 f 4 -f 1 f 3 ) -q (f 3 -f 1 ) ρ 2 16 -2 q ρ 4 (f 3 -f 1 ) (f 2 + f 4 ) + O( 3 ).
Here and in the following, we only treat f 1 and f 3 , similar equations hold f 2 and f 4 .

Recasting the terms of the same order in , one has

1 : ∂ x ρ 4 = p ρ 4 1 - ρ 4 (f 2 + f 4 -f 1 -f 3 ) + q (f 3 -f 1 ) ρ 2 16

Figure 1 .

 1 Figure 1. (a) Collision rules of the model. The two particles collision occurs with a probability p, and three particles collision occurs with a probability q. (b) Bottom boundary condition: at y * = 0, channels 1, 2 and 3 are filled randomly with probability ρ B /4 and particle in cell 4 come from the bulk so ρ B represents the density of the bottom reservoir. Similar conditions are applied on the other boundaries with different densities for the reservoirs. (c) View of the lattice of L × L nodes between reservoirs at densities ρ L (x * = 0), ρ B (y * = 0), ρ R (x * = L) and ρ T (y * = L).

Figure 2 .

 2 Figure 2. Simulation for L = 300, p = q = 1 ρ L = ρ B = 2.4 and ρ R = ρ T = 1.6. Left: instantaneous microscopic density field in stationary state. Right: average density field measured in stationary state over 10 6 time steps.

Figure 3 .

 3 Figure 3. Instantaneous coarse-grained density field for our simulation with a coarse graining factor τ = 5.

Figure 4 .

 4 Figure 4. Empirical PDFs of the variables δ τ α computed in our simulation for various ĵ, ρ, ĝ, τ = 10, 15, 20, 25, 30 and α = x, y. The Normal law behaviour is indicated by a blue thick line.

Figure 5 .

 5 Figure 5. Left: relaxed coarse-grained current µ τ α as a function of the coarse-grained density and its gradient for τ = 10 in our simulation. Right: theoretical relation obtained by the Chapman-Enskog expansion for the averaged observables.

Figure 6 .

 6 Figure 6. Relaxation time r τ divided by the coarse graining factor (left) and the root mean square of the variation of the coarse-grained current σ τ (right) as functions of the coarse-grained density and its gradient for τ = 10 in our simulation.

Figure 7 .

 7 Figure 7. Left: τ multiplied by τ 5/2 as a function of the coarse-grained density. The model (35) (in full lines) and the root mean square of the microscopic current (in light blue full line) are compared to the numerical results (in stars) for various values of τ . Right: Ratio between the numerical value of σ τ and our model R τ .

Figure 8 .

 8 Figure 8. Relaxation time r τ divided by τ for our simulation with various values of τ (represented with stars) compared to our model (40) (represented with full lines).

Figure A1 .

 A1 Figure A1. Comparison between the stationary density field solution obtained by a Laplace solver ρ theo and the average density field from simulation ρ. Left: the average density computed over 10 6 time steps. Right: difference with the solution obtained by solving the Laplace equation ρ theo .

Figure A2 .

 A2 Figure A2. Left: Relation between the average current j α and its theoretical prediction (A.25): -D(ρ)∂ α ρ. The line j α = -D(ρ)∂ α ρ is shown in yellow. Right: Relation between the variance of the microscopic current and local average density. The theoretical relation (A.30): σ 2 * (ρ) = ρ 2 1 -ρ 4 is shown in yellow. Blue circles represent α = x and orange circles represent α = y.

Figure A3 .

 A3 Figure A3. Relative error coefficients defined in (A.35-A.38) that quantify the departure from the Lattice Boltzmann Approximation for our simulation. The fields presented here are smoothed by a Gaussian filter.

Acknowledgments

This work has been supported by the ANR EXPLOIT, grant agreement no. ANR-16-CE06-0006-01. We thank the anonymous referees for their useful comments that have helped to improve significantly the article.

Appendix A. Hydrodynamic limit and microscopic fluctuations of the current Appendix A.1. Average current For a LGCA (i.e. for a given lattice and collision rules), one can derive the macroscopic behaviour of observables at large scales via the equation verified by the average quantities (here ρ and j) in the limit of infinite lattice size. The most popular way to it is to apply the multi-scale analysis used in the continuous kinetic theory (the so-called Chapman-Enskog expansion) to the discrete lattice gas. Here, a simple diffusive scaling is sufficient to derive the macroscopic equation for the evolution of the average density. To do this, we start with the local micro-dynamic equation for the Boolean field which is given by the Chapman-Kolmogorov equation [START_REF] Frisch | Lattice gas hydrodynamics in two and three dimensions[END_REF][START_REF] Rivet | Lattice Gas Hydrodynamics[END_REF]. For our simple model with only four velocities, we can explicitly write it: with

)

)

where ζ p (ζ qx or ζ qy ) represents a Bernoulli random variables with mean p (q), that are decorrelated in space and time, and with the n i 's. We have omitted the spatio-temporal dependence (ζ(x * , y * , t * ) → ζ and n i (x * , y * , t * ) → n i ) in the three last equations for brevity.

We now use three assumptions:

(i) The lattice Boltzmann (Mean Field) approximation: the correlations between particles before collisions are neglected. This allows to write

where I is any set of non repeated indices of the channels.

Using (A.14) and (A.15) with (A.11), we show that

With the same procedure,

Substituting (A.18) in (A.16), we obtain

and, similarly,

Summing the two last equation and replacing f 1 and f 2 by their expressions in term of density, one obtains the evolution equation for the average density

The stationary solution of the density field can be obtained by solving the Laplace equation for ρ/4 + 4/qρ. We compare this solution to the results of the simulation with ρ L = ρ B = 2.4, ρ R = ρ T = 1.6 and deterministic collisions p = q = 1 (see figure A1).

We naturally observe that the errors are more important near the corners where density of boundary reservoir is discontinuous. The equation (A.22) corresponds to a non-linear diffusion equation

with diffusivity

The diffusivity has two contributions. The factor -1/4 is a spurious (anti-diffusive) term called propagation diffusivity [START_REF] Hénon | Viscosity of a lattice gas[END_REF] which emerges from the discrete nature of the velocity set. The term 4 qρ 2 will be dominant for q ρ 4 2

1 and is due to the three particles collisions. The hydrodynamic current of particles is then

The last equation has been checked numerically (see figure A2). The points that don't fit well with the theoretical relation correspond to the nodes near the boundaries.

Appendix A.2. Microscopic fluctuations of the current

Using the Boltzmann approximation (A.8), it is easy to write the probability to observe a local state s = (s i , i = 1, ..., 4) at node r * and time t * using the definition of the mean occupation of channels N i (r * , t * ):

We can then write the fluctuations of the two components of the current

Using the small perturbation hypothesis (A.11) and keeping only the zero'th order in , we obtain lim

We point out that the fluctuation function σ is the same as the one computed in the Simple Symmetric Exclusion Process (SSEP) (see [START_REF] Derrida | Microscopic versus macroscopic approaches to non-equilibrium systems[END_REF]) and is typical for models with exclusion. Numerically, the current verifies the relation (A.30) well (see figure A2).

As for the mean current, the points that don't fit well with the theoretical relation correspond to the boundary nodes.

Appendix A.3. Violation of the Boltzmann hypothesis

In practice, the particles can be correlated before collision. To quantify these correlations, we look at the local observables which are the relative errors (compared to the Lattice Boltzmann Approximation) of n p (n qx , n qy and j * x j * y ). These coefficients are evaluated at each node (see figure A3). For re p and re qx the fluctuations are important. The statistics of re p are homogeneous inside the system, except near the boundaries where there is a discontinuity of density of the reservoirs. It is globally positive. However, the statistics of re qx are not homogeneous and we observe important fluctuations near the top and bottom boundaries. These correlations between particles are not significantly modified when we change the size of the lattice. A similar behaviour is observed for re qy .


Jérôme Pocachard

Nathalie Klement

Christophe Jouve

Anthony Quenehen

Methodology to select the best part presentation in cobotics

Keywords: Part presentation, collaborative robot, methodology, Industry 4.0

come   L'archive ouverte pluridisciplinaire

Introduction

A new industrial revolution towards Industry 4.0 has been ongoing for ten years with the introduction of new technologies. Automation, integration, collaboration, flexibility, safety and security are the main fields of research attention [START_REF] Liao | Past, present and future of Industry 4.0 -a systematic literature review and research agenda proposal[END_REF]. New technologies are also going to impact the management of shopfloor operators in terms of requirements and demand [START_REF] Holm | The future shop-floor operators, demands, requirements and interpretations[END_REF]. In a new paradigm of reconfigurability requested by the markets [START_REF] Koren | The rapid responsiveness of RMS[END_REF], system has to evolve. Collaborative robots, named cobots, have been clearly identified as promising tools for enabling human operators to bring flexibility to the robotized workplace [START_REF] Jocelyn | Collaborative Robotics : Assessment of Safety Functions and Feedback from Workers , Users and Integrators in Quebec (Report 1030)[END_REF] [5] [START_REF] Schonberger | Modeling Workflows for Industrial Robots Considering Human-Robot-Collaboration[END_REF] and to increase efficiency and ergonomics [START_REF] Realyvásquez-Vargas | Introduction and configuration of a collaborative robot in an assembly task as a means to decrease occupational risks and increase efficiency in a manufacturing company[END_REF]. The existing literature focuses on technical solutions to increase safety, productivity and decrease costs [START_REF] Faccio | Collaborative and traditional robotic assembly: a comparison model[END_REF]. In this article, a mathematical model was made comparing collaborative and traditional robotic assembly by measuring the objectives of throughput and unit direct production costs. However, Faccio et al. [START_REF] Faccio | Collaborative and traditional robotic assembly: a comparison model[END_REF] concluded that development costs should be addressed by further research. The shop-floor operators in a human-robot collaboration have to be able to change the full manufacturing configuration [START_REF] Schou | A plug and produce framework for industrial collaborative robots[END_REF] quickly. As though by Gambao et al. [START_REF] Gambao | A new generation of collaborative robots for material handling[END_REF] through their new generation of human-robot solutions, new paradigms and approaches have to be developed for handling technology.

Collaborative robots allow a human-robot collaboration (HRC). This HRC takes place inside workstations manufacturing semi-finished or finished products, which are fed from warehouses. Bortolini et al. [START_REF] Bortolini | Assembly system design in the Industry 4.0 era: a general framework[END_REF] divided the material feeding, when dealing with the design of the assembly operation in the frame of Industry 4.0, into storage (packaging type and dimensions) and feeding policy (line side stocking, kitting feeding, or kanban feeding). Once the part has been fed to the workstation, the part must be presented to the manufacturing resources. In the paper context, the resources are specifically cobots, performing some operations of picking and placing. To perform these operations, a cobot needs to grasp parts reliably and efficiently. The parts are presented in the workstation thanks to a part presenter, waiting for the picking by the cobot. The question comes regarding the best selection of part presenter.

Boothroyd [START_REF] Boothroyd | Use of Robots in Assembly Automation[END_REF] defined the purpose of part presentation as to present the parts in the same orientation at the same location. Indeed, the cobot itself with a simple gripper cannot grasp the part if the part is not in a given position and orientation. However, with the recent development of new technologies, some hightech devices can give some senses to cobots, through vision systems or torque sensors. Part presentation is consequently important for the efficiency of the cobot system. Best practice would be to include the design of the part presenter during its development at the design phase to avoid rework, delays and cost. Indeed the cost of part presentation can significantly influence the total manufacturing costs by 10%-30% [START_REF] Bi | Flexible fixture design and automation: Review, issues and future directions[END_REF] [START_REF] Lee | Flexible part-feeding system for machine loading and assembly. Part I. A state-of-the-art survey[END_REF]. Part dedicated-fixtures do not bring the flexibility requested to deal with a variety of work-pieces. So improving fixtures can help reduce unit costs. Moreover, the optimum fixture can evolve during the lifespan of the part. The fixture at the start of a new product may differ from maturity step.

Consequently, part presentation is a link between part feeding and the manufacturing workstation. The introduction of new technologies through Industry 4.0 influences part presentation. Moreover, reconfigurability is required by the markets demanding frequent product redesigns, so shorter part life cycle. Collaborative robots are a response to it. The objective of the paper is to determine if a shorter part life cycles impacts the choice of the best part presenter, which feeds a cobot. Reconfigurability being an objective, the performance criterion to select the best part presenter is the time spent due to a new part presentation device (presenter redesign, cobot programming, change of part loading). The methodology used for the research is to develop a complete method (flowchart, inputs, processes, outputs) selecting the best part presenter, and then to test it with a concrete case. The paper is organized as follows. Section 2 presents a literature review. Section 3 organizes the characteristics of part presenters, parts, and manufacturing operations. Based on these characteristics, a method is suggested in Section 4 to select the best part presenter. In Section 5, this method is tested through an industrial use case. Section 6 concludes this paper and points out the stages to follow.

Literature review

The literature review starts with part feeding, which is the process before part presentation as explained above. Then, literature review is extended to other relative topics, which are part disorder and fixture & jig. Motivations for these topics, definitions and literature reviews are detailed hereunder.

Part feeding is transporting a part from the warehouse into the manufacturing workstation. The existing literature concerning part feeding is really intensive [START_REF] Kilic | Advances in assembly line parts feeding policies: A literature review[END_REF]. Linked with this paper objective, Battini et al. [START_REF] Battini | Design of the optimal feeding policy in an assembly system[END_REF] and Caputo et al. [START_REF] Caputo | Selection of assembly lines feeding policies based on parts features[END_REF] evaluated the best feeding policy (kitting, line storage, and just in time delivery) according to the part features. Part presenter can be supplied by one of the three above part feeding strategies. Hanson et al. [START_REF] Hanson | Assembly station design: A quantitative comparison of the effects of kitting and continuous supply[END_REF] experimented with the time spent for fetching parts and for assembling parts between the feeding principles of continuous supply and kitting. Time spent searching parts and distance are two factors influencing the fetching time.

Literature concerning part presenter and its design choice is cost-oriented. Boothroyd [START_REF] Boothroyd | Use of Robots in Assembly Automation[END_REF] described a decisional economical graph to select the right part-presentation system. The costs taken into account were the investment cost (general-purpose equipment plus part-dedicated tooling) and the labour cost. Four different automatic feeder systems were studied: vibratory bowl feeder, Hitachi multipart feeder (vibratory system and vision system), Salford multipart feeder (interchangeable mechanical tooling) and double-belt feeder. Based on cycle time and quantity to be manufactured during payback period, the best feeder system among the four can be selected on a graph calculated by the authors. Below 6500 pieces, none of the automatic feeders was economically interesting. Ho and El-Gizawy [START_REF] Ho | A Programmable, Multi-Part Presentation System for Robot Assembly[END_REF] also performed an analysis on economic interests based on the costs of basic equipment, tooling and changeover.

When a part is dropped off on a surface as in a part presenter, the part can be located in various positions. Positions can be from well-defined ones to fully random ones if the degrees of freedom are not constrained. Sanderson [START_REF] Sanderson | Parts entropy methods for robotic assembly system design[END_REF] studied mathematically part disorder by defining the concept of part entropy. Chirikjian [START_REF] Chirikjian | Parts entropy and the principal kinematic formula[END_REF] extended this concept to the multiple interacting parts. The entropy H was calculated from the probability distribution of part position and orientation. If H is equal to zero, there is no ambiguity of the part position.

Fixtures and jigs are sometimes mentioned in literature when searching publications on part presenters. Nevertheless, fixtures are really different to part presenters as they clamp one or few parts with technical devices. The functions of a fixture are locating, supporting and clamping, before manufacturing or measurement operations. Bi and Zhang [START_REF] Bi | Flexible fixture design and automation: Review, issues and future directions[END_REF] detailed the design of flexible fixtures. The fixture has not to be mixed up with jigs, which are devices to guide another tool for repeatability purposes (e.g. drill). The flexible strategies for fixture systems are either a modular structure (modules can be adapted depending on the part), or a single structure but with material phasechange (e.g. temperature-induced) or with adaptive clamps (internal variables fitting with various parts). As these fixtures are somewhat complex, dynamic simulation can help to design the part feeders and its configurations [START_REF] Berkowitz | Designing parts feeders using dynamic simulation[END_REF]. The fixture design and its computer-aided tools were extensively addressed by the researchers [START_REF] Bi | Flexible fixture design and automation: Review, issues and future directions[END_REF].

As a conclusion of the review, as opposed to part feeding literature, part presentation has been less attractive for the researchers for a long time. This can be easily explained: traditional robots are used for high volumes of production, so the choice for part presentation is obviously some sophisticated devices to achieve very high throughput. But it is not the case with cobots dealing with low personalized throughput. On top of that, cost impacts are often addressed, meanwhile the impact of development has not been addressed. The novelty of the study is to explore the impact of this development time on part presenter choice. Indeed, development is a key point for reconfigurability objective as referenced in the introduction section. Consequently, this paper aims to continue previous research to select the proper part presenter, for cobots in the scope of In-2 dustry 4.0, including the development of the part presenter as a new criterion of selection.

Manufacturing characteristics

The characteristics of the part presenter

Boothroyd [START_REF] Boothroyd | Use of Robots in Assembly Automation[END_REF] breakdown the part presentation into two categories: single location (robot grips part always at the same location) or multiple location (robot grips part based on a predetermined pattern or randomly thanks to a vision system). Applying the part entropy methodology [START_REF] Sanderson | Parts entropy methods for robotic assembly system design[END_REF], these two categories (single and multiple location) have different impacts on part entropy and reconfigurability. This leads to switching from Boothroyd's two categories to three (single, multiple, random). The study keeps this taxon 'location' with these three possibilities.

In order to select the best part presenter, the first step is obviously to know the potential part presenters already used by robots and cobots worldwide. Fig. 1 classified them using a Function Analysis System Technique (FAST) graph. A FAST graph organizes the functions that need to be performed by the process into "how?/why?" relationships. The right side of the graph consequently lists the existing technical solutions to reach the part presentation objective. The loading of the parts into the presenter can be either on-line (no stopping of the robot) or off-line (stopping of the robot required). Additionally to the taxon 'location' previously explained linked to part entropy, the taxon 'energy' is first introduced on the left side. Indeed, the graph is divided into two sets depending if energy is required. The following paragraphs list and detail this new classification of the part presenters using the taxons 'energy' and 'location'.

NO ENERGY REQUIRED:

The following part presenters are systems without electronic or moving devices. The design is only mechanical, quite quick (a few hours), and cheap (hundreds of euros). 3D printing can produce them. Their contents (number of parts) are limited. Possible part presenters are: a) One-piece structure, to present one part at once. It is the simplest way; b) Gravity track, a structure so that parts can easily slide once part at the bottom is picked up. It is one dimension structure, and may be two-dimension one if the part is cylindrical; c) Tray, a two-dimension structure; d) Box, a three-dimension structure. It is mandatory that parts can be stacked one above the other.

REQUIRING ENERGY:

The following part presenters require electrical or pneumatic energy. Indeed, they are complex systems with different moving devices. Their design is consequently quite long (a few weeks) and expensive (some thousands of euros or even more). Qualified multidisciplinary f) Vibratory, rotating and mechanical feeders are specialtooled systems. The bulky parts in a container are then positioned 1-by-1 through a vibration or well-engineered mechanism and then moved to the workstation through an escapement. This system is dedicated to a given part (driven by the geometry and weight barycentre of each part). A vision system may be integrated; g) Tray magazine loader can present a new full tray when the previous one has been emptied by the process; h) Conveyor carries linearly the parts (randomly oriented and positioned) on a belt until a specific place triggered by sensors. An appropriate 2D vision system detects the parts and transmits the positions precisely to the cobot so that the cobot can pick up each part. The parts at the entry of the conveyor or directly below the vision system may come from a vibratory plane; i) Bin picking: a bin contains the bulky parts. The gripper of a cobot picks up appropriately each part through a 3D vision system. It leads to complexity to avoid some picking errors or to crushing parts.

A part presenter with energy is technologically more complex but is more efficient. It will consequently be dedicated to parts with a high volume of production and low variety, to avoid reworking the design of the presenter. This energy-required part presenter is an antagonist to the flexibility targeted with the usage of cobots in HRC. Nevertheless, a vision system may reduce the reworking impact as only reprogramming would be requested. Traditional industrial robots usually perform operations in a context of fixed special purpose automation dedicated to high production volumes and a low variety of products.

Whereas, the robots associated with human operators (semi automation) are efficient for mid production volumes and a mid variety of products [START_REF] Heilala | Modular reconfigurable flexible final assembly systems[END_REF]. It is specifically adapted to shorter life cycle products in the current economy [START_REF] Bannat | Artificial cognition in production systems[END_REF]. For this semi automation, which is clearly the objective of collaborative robots, the design of the parts and the finished products evolves quickly. Consequently, the paper focuses only on the part presenter not requiring energy. On top of that, doing part presentation without energy contributes to the energy challenges of the 21 st century.

The characteristics of the parts

The objective of the part presenter is to present the part easily for the gripper. Each part has its own graspability, i.e. the qual- [START_REF] Malik | Framework to Implement Collaborative Robots in Manual Assembly: A Lean Automation Approach[END_REF] did a systematic framework for cobot deployment in assembly cells. Their process requirements include the definition of the physical properties of the parts: form, shape, dimensions, tolerances, material, weight, size and surface finish, etc. These physical properties are added to two other ones more specific to part presenters, leading to the following three types of part characteristics:

• Physical properties: shape, dimensions with their tolerance intervals, material, weight, surface finish

• Properties making presenter design easier: narrow tolerances to locate accurately and isostatically, shape allowing error-proofing system, symmetric part (cylindrical etc.)

• Exclusive properties for some presenters: regular shape for pattern logic, stackable, can slide on an inclined surface, cannot overlap with other identical part, cannot tangle with other identical part These latter exclusive properties drive the choice for some part presenters. For example, parts must be stackable to be stored in box. Table 1 presents the analysis for each part presenter.

The characteristics of the manufacturing workstation

When designing the workstation and its part presenters, one important data is the expected cycle time of the workstation, as a step in a global workshop. It is linked with the takt time from Lean methodology. Moreover, depending on the product routing and the abilities of the agents (i.e. humans and robots), a level of collaboration between them is defined. The levels of collaboration between the human worker and the cobot can be multiple [START_REF] Aaltonen | Refining levels of collaboration to support the design and evaluation of human-robot interaction in the manufacturing industry[END_REF]. Bauer et al. [START_REF] Bauer | Lightweight robots in manual assembly -best to start simply, Examining companies' initial experiences with lightweight robots[END_REF] defined four typologies: coexistence, synchronized, cooperation, or collaboration. It impacts the requirement regarding workstation autonomy as analyzed in Table 2 in terms of workstation autonomy. The workstation autonomy (t workstation-autonomy ) is the time during which the workstation can work without human intervention. Based on cycle time and the level of collaboration, the optimum workstation autonomy can be computed.

Another requirement concerning the manufacturing workstation is the available physical area for part presenters. The area is key data as the usage of cobots improves the footprint, by avoiding cage fencing when applications are safe. A reduced area, leading to shorter and quicker flows and lower fixed costs, is a strategic advantage. However, a reduced area for part presenters with fewer parts would decrease workstation autonomy. The optimization of the area is measured with the indicator named Area Efficiency as given by Equation [START_REF] Liao | Past, present and future of Industry 4.0 -a systematic literature review and research agenda proposal[END_REF].

Area E f f iciency = Area one-presenter Area one-part × c (1)
where (c) is the capacity of parts inside one part presenter.

Method

Before explaining the methodology, the paragraphs will explain the development time and lifespan load time.

The development time

When a new presenter has to be (re)designed due to part evolution, there is a development time (DT ). In a previous study [START_REF] Quenehen | Process optimisation using collaborative robots -comparative case study[END_REF], we divided a project into process engineering, implementing and operating steps. Development time is the duration of process engineering and implementing operations. For a new part presenter, Equation (2) concerns the design of the part presenter (t design ) and the programming of a cobot (t programming ). These two times are directly linked to the complexity of a part presenter.

DT = t design + t programming (2)

The lifespan load time

During the manufacturing cycle, parts are grasped by the cobot. So the number of parts in the part presenter decreases until it is empty. At that time, as shown in Fig. 2, the part presenter is to be reloaded, either by replacing it by a new full presenter or by refilling it with parts depending on the design of the presenter. This reloading is called a rotation of part presenter. Equation ( 3) calculates the number of rotations of the part presenter (r).

r = q/c (3)
where lifespan quantity (q) is the number of the same parts, which are dropped-off on the part presenter during part's lifespan, i.e. until its redesign.

Each part presenter requests a time (t li f espan-load ) to load it fully during all rotations (r) as in Equation (4).

t li f espan-load = r × t one-rotation (4)
where t one-rotation is the time to change the part presenter for an offline load, and the time to load q parts in the presenter for an online load. For the purpose of the calculation, all part presenters are considered to have the same exposure for the picking (angle, height). The dichotomy from Boothroyd et al. [START_REF] Boothroyd | Part Presentation Costs in Robot Assembly[END_REF] is used for the study: picking with one hand vs. picking with two hands. The picking time is valued with an average of previous studies [START_REF] Boothroyd | Part Presentation Costs in Robot Assembly[END_REF] [30] [START_REF] Finnsgård | Factors impacting manual picking on assembly lines: An experiment in the automotive industry[END_REF], so 2 seconds to load one part with one hand and 4 seconds if hand loading requires two hands. The autonomy of the workstation, i.e. until the part presenter is empty, is defined in Equation [START_REF] Matthias | Safety of collaborative industrial robots: Certification possibilities for a collaborative assembly robot concept[END_REF].

t workstation-autonomy = c × CT ( 5
)
where CT is the cycle time of the manufacturing operation itself. Operating lifespan (OT ) can be calculated as in Equation ( 6) if one part is requested at each cycle. Considering DT and OT , the times impacting the choice of the presenter are t design , t programming , and t li f espan-load . Only t li f espan-load depends on the lifespan quantity of the part.

OT = t li f espan-load + CT × q (6)

Methodology to select the part presenter

The data defined in the above paragraphs are used through a methodology to select the part presenter as detailed in Fig. 3.

Step 1: a first selection is made based on the characteristics of the parts and the exclusive properties for some presenters. As explained previously, some presenters are not possible due to some technical characteristics.

Step 2: the required capacity of the presenter is a key data. It is determined based on the characteristics of the workstation, in terms of autonomy and cycle time.

Step 3: this last step allows us to select the appropriate part presenter. Development and programming times are calculated to evaluate their impacts during the reconfigurability period.

Results from use case

In order to check the methodology of part presenter selection, a concrete industrial use case has been used. The manufacturing process is the assembly of a pneumatic cylinder as detailed by Quenehen et al. [START_REF] Quenehen | Process optimisation using collaborative robots -comparative case study[END_REF]. Its bill of material has six components, including a cylinder barrel. A cylinder barrel is selected for the study of part presentation because, as explained in Step 1, there are no blocking characteristics to design the four part presenters (one-piece, gravity, tray and box). Also, its shape is simple enough to avoid any impact of a complex grip-per on the design of the part presenter. The workstation is composed of a collaborative robot sharing the assembly tasks with a human operator. The objective is to use the above methodology and analyze how the various part presenters differentiate themselves.

Step 1: The characteristics of a cylinder barrel are assessed in Table 3 along the criteria defined in the methodology, including the exclusive part properties for some presenters from Table 1. Based on this comparison, Step 1 concludes that all part presenters are technically possible for the concerned part, which is a cylinder barrel.

Step 2: This step is to calculate the required capacity of the part presenter. The sequence of assembly is already determined and shared between a collaborative robot and a human operator. The level of collaboration is classified as cooperation based on Table 2. Consequently, the capacity should be sufficient to match with the cycle time. The cycle time to assemble one cylinder is known to be 53 s with a waiting time of the operator at 5 s, which is compatible with the time to change a presenter with two hands (4 s based on the literature review previously quoted). Based on the upstream flow of parts inside the plant, the autonomy of the workstation is required to be 420 s. Following the method, the required capacity of the presenter is 420 s divided by 53 s, so 7.9 parts. This capacity is used later on at

Step 3 when the design of the part presenter is known and can be compared with this required capacity.

Step 3: The four potential presenters have been computeraided designed (CAD) as reported in Fig. 4: one-piece (a), gravity track (b), tray (c), and box (d). Supervised fifth-year students in mechanical engineering with a basic background in CAD software and robotics have done the designs and evaluated the programming times of the cobot and the time to design the part presenters. Table 4 summarizes all these data. Concerning the required capacity from Step 1 (7.9 parts), obviously a one-piece structure cannot satisfy it. The workstation, due to space restrictions on the shop floor, has a possible area of 300 mm x 200 mm for the presenter of cylinder barrel. The tray fits this area but its capacity is below the required one. The box is fully dimensional OK and has the advantage to optimize area (higher area efficiency) thanks to the two layers of parts. The gravity track has its length exceeding the required area length and is below the required capacity. The three impacting times following the method (design, programming, lifespan-load) have been calculated as per Fig. 5. In this figure, the best presenter is clearly linked to lifespan quantity. For a quantity of thousands, the onepiece structure or tray provides the less time-consuming solution. For a quantity of a few tens of thousands, boxes or trays are the best ones.

As a conclusion, each presenter has their own advantages and drawbacks. The tray looks to be the best compromise between all parameters, while the gravity track is the worst one. The box is very sensitive to lifespan quantity. If a high quantity is anticipated in the short term, it could be the best presenter.

Conclusions, limitations, further development

Industry 4.0 is ongoing worldwide, bringing new technologies such as cobots. At the same time, commercial markets request higher flexibility. Manufacturing has to manage in their workstations smaller product volumes, higher numbers of product variants, and more frequent product redesigns. In manufacturing workstations, cobots have to grasp parts from part presenters. Consequently, part presenters contribute to this flexibil-ity. Previous studies, as detailed in the literature review Section, have addressed part presentation through economical analysis. This paper continues it by including the time of development as a questionable parameter. Indeed, this parameter influences the presenter redesign when the product is evolving. A methodology is proposed in order to list and analyze all parameters influencing the choice of the best part presenter. After having structurally listed all various part presenters, the study focuses only on presenters not requiring energy for flexibility and sustainability reasons. Regarding the objective of research (impact of shorter part life cycle on the choice of the best part presenter), the study demonstrates that lifespan quantity affects the choice of the presenter. The presenter is to be adapted all along the life cycle of the part, to be always the optimal one. An appropriate compromise is the usage of a tray as a part presenter, especially when the part life cycle is difficult to predict.

When considering the three impacting times (design, programming, lifespan-load), the resources of these tasks have not been differentiated. In a classical organization, the design would be performed by the manufacturing engineering department meanwhile programming and loading by the manufacturing operators and technicians. Further development would consider either the impact of different resources on the proposed methodology, or even investigate how to execute these tasks solely by manufacturing operators. It is perfectly aligned with the objective of Industry 4.0. Indeed this would lead to empower humans and let them focus on the true added value from their competencies and versatility versus robot's ones. Also, in this paper, part presenters requiring energy have not been addressed, as explained for reconfigurability and sustainability reasons. However, recent relevant advances on grasping objects using image processing and artificial intelligence address reconfigurability objectives [START_REF] Bergamini | Deep learning-based method for vision-guided robotic grasping of unknown objects[END_REF]. Future research could include part presenters requiring energy in the developed method, while evaluating sustainability objectives too.

Fig. 1 .

 1 Fig. 1. Function analysis system technique. (1) Cannot be supplied online

Fig. 2 .

 2 Fig. 2. The operating flow.

Fig. 3 .

 3 Fig. 3. Methodology to select the appropriate presenter.

Fig. 4 .

 4 Fig. 4. Design of presenter with parts inside.

Table 1 .

 1 Impact of the exclusive part properties.

	Exclusive part characteristics	a) One-piece	b) Gravity track	c) Tray	d) Box
	Regular shape for pattern logic	Optional	Mandatory	Mandatory	Mandatory
	Stackable	Optional	Optional	Optional	Mandatory
	Slide on an inclined surface	Optional	Mandatory	Optional	Optional
	Cannot overlap with other identical part	Optional	Mandatory	Optional	Optional
	Cannot tangle with other identical part	Optional	Mandatory	Optional	Optional

Table 2 .

 2 Workstation autonomy based on the level of collaboration.

	Level of col-	Situation of human	Workstation auton-
	laboration	and cobot		omy
	Coexistence	Work	alongside	High. The work-
		but	not	shared	station works until
		workspace. Human	the part presenter is
		does not interact	empty
		with cobot during
		operating time
	Synchronized Work sequentially on	One-piece flow. The
		one given product	human can feed from
					bulk the workplace
	Cooperation	Shared workspace,	Medium. The cobot
		non-simultaneous	needs to be au-
		tasks,	separate	tonomous but if
		product		needed the human
					can feed from bulk
					the workplace
	Collaboration Work simultaneously	One-piece flow. The
		on the same product	human can feed from
					bulk the workplace

Table 3 .

 3 Characteristics of the part, a cylinder barrel.

	Shape	parallelepiped
	Dimensions	60 mm x 60 mm x 76 mm
	Tolerances	0.1 mm
	Material	Aluminum
	Weight	232 g
	Surface finish	exterior: extruded, interior: machining
	Narrow tolerances to locate accurately	Yes
	Shape allowing error-proofing system	Yes
	Symmetric part (cylindrical etc.)	Partly
	Regular shape for pattern logic	Yes
	Stackable	Yes
	Slide on an inclined surface	Yes
	Cannot overlap with other identical part Yes
	Cannot tangle with other identical part	Yes

Table 4 .

 4 Data for the presenters designed in the use case.

	Presenter	Number	Capacity Length	Width	Height	Footprint	Area effi-	Design	Program.
		of		(mm)	(mm)	(mm)	(cm 2 )	ciency	time	time (s)
		parts							(s)
	a) One-piece	1x1	1	70	70	83	49	73%	7200	360
	b) Gravity track	6x1	6	350	63	130	221	98%	14400	360
	c) Tray	3x2	6	260	155	83	403	54%	10800	1800
	d) Box	2x2x2	8	182	157	112	286	101%	18000	7200