Riccardo Artoni

Giovanni Loro

Patrick Richard

Fabio Gabrieli

Andrea C Santomaso

Fabio Gabrieli Geotechlab

Drag in wet granular materials

Keywords: wet granular material, drag, mixing, shear testing

come

Introduction

The study of forces induced on bodies immersed in granular media has several applications e.g. mixing of granular materials and pastes, predicting damage induced by landslides on structures and designing retaining or dissipative systems for geohazards. In all the cases, there is a strong need for the development of (1) conceptual models and (2) new sensors for monitoring and control. In this perspective, the measurement of the drag force exerted by the material on a moving intruder (whose size is much larger than the particle diameter) can be seen as an opportunity in order to gather information about the rheology of granular materials and their interaction with structures. For these reasons, in the last decades there has been interest on characterizing the drag force on objects in dense granular media. [START_REF] Wieghardt | Forces in granular flow[END_REF][START_REF] Wieghardt | Experiments in granular flow[END_REF] was one of the first to study "granular" drag. His experiments consisted in moving partially immersed rods through sand beds. The drag force is found to be dependent on both the geometry of the rod and its velocity although the latter parameter has to be relatively large to have a significant effect. The behavior at low velocity can be understood by assuming that the interactions between the rod and the surrounding grains are mostly governed by sliding friction, a velocity independent force. In contrast, collisions between the rod and the surrounding grains lead to a velocity dependent dissipation. They are no more negligible at important velocities and can explain the observed weak influence of the velocity.

A large amount of work in this area has been performed by Albert and coworkers [START_REF] Albert | Slow drag in a granular medium[END_REF](Albert et al., , 2001a,b),b). In particular, they first investigated the effect of velocity, particle size, object size and depth on the resistance provided by dense granular materials on a partially immersed, slowly moving cylinder [START_REF] Albert | Slow drag in a granular medium[END_REF]. The symmetry axis of the cylinder was vertically oriented and the motion was horizontal. Subsequently, they focused on the effect of the shape of the dragged object on the drag force and the relevant scalings (Al-bert et al., 2001a). The main results which were obtained in these works are that the drag is independent of velocity (in the quasistatic and dense flow regimes), and particle diameter (at least as soon as the intruder size is large with respect to the particle diameter), but depends on the depth of the object by a scaling which depends on object shape. Moreover, the drag force appears to increase linearly with the object size. These results stem from the frictional nature of granular materials, in which the yield stress is proportional to pressure. In a further work (Albert et al., 2001b), they also characterized drag force fluctuations, which they linked to jamming events occurring in the bulk of the material.

Later, [START_REF] Chehata | Dense granular flow around an immersed cylinder[END_REF] studied the drag force acting on a cylinder immersed in a dense, quasi twodimensional vertical chute flow. They varied the flow rate, the object diameter and particle size, finding a rateindependent behavior, a linear dependence on cylinder size and a weak dependence on particle size. By means of 2D numerical simulations, [START_REF] Wassgren | Dilute granular flow around an immersed cylinder[END_REF] studied the flow of a dilute granular system around a disk. They found that the drag coefficient decreases with decreasing normal coefficient of restitution and is relatively insensitive to the friction coefficient. Interestingly they also reported shock and expansion fan structures, similar to those observed in gas dynamics. [START_REF] Seguin | Dense granular flow around a penetrating object: Experiment and hydrodynamic model[END_REF] coupled drag force measurements with Particle Image Velocimetry (PIV) to characterize the case of a horizontal cylinder moving into a granular medium. The drag force behavior is in agreement with previous works, apart from a nonlinear dependence of the force on particle size, which seems to be present when the latter is below 1 mm. They compared the experimental velocity profiles around the object with predictions from an extended kinetic theory. In a different configuration, a rod rotating at a given depth in a box filled with grains, [START_REF] Guillard | Depthindependent drag force induced by stirring in granular media[END_REF] obtained that the torque needed for stirring the material increased at the beginning, then dropped down, eventually reaching a stationary value which appeared to be independent on the depth. Clearly, the scaling on depth is completely different from the literature referenced above, but may be ascribed to the peculiar geometry used by the authors. More recently, [START_REF] Kumar | Scaling law of the drag force in dense granular media[END_REF] numerically studied the scaling of the drag force on velocity. They extended previous studies by including flows well beyond the slow quasistatic regime, and identified the existence of three regimes (rate-independent, linear and quadratic dependence of the drag force on velocity).

All the works referenced above have dealt with dry, cohesionless granular materials. Rheology of partially wet conditions have been discussed with reference to rheometer experiments on plastic beads [START_REF] Badetti | Rheology and microstructure of unsaturated wet granular materials: Experiments and simulations[END_REF], rapid shear cell tests with glass beads [START_REF] Liao | Experimental analysis of dynamic properties in wet sheared granular matter[END_REF], quasistatic tests in a Schulze shear cell with glass beads [START_REF] Louati | Apparent friction and cohesion of a partially wet granular material in steady-state shear[END_REF][START_REF] Louati | Effect of the particle size and the liquid content on the shear behaviour of wet granular material[END_REF] and column collapse tests with glass beads [START_REF] Artoni | Collapse of quasi-two-dimensional wet granular columns[END_REF][START_REF] Gabrieli | Discrete particle simulations and experiments on the collapse of wet granular columns[END_REF] and different types of sand [START_REF] Santomaso | Collapse and runout of granular columns in pendular state[END_REF]. [START_REF] Collet | Mixing of a wet granular medium: Effect of the particle size, the liquid and the granular compacity on the intensity consumption[END_REF][START_REF] Collet | Mixing of a wet granular medium: Influence of the liquid addition method[END_REF]) also studied the effect of liquid addition during mixing of fine powders and glass beads by analysing the electric current absorbed by the mixer.

An analogous approach in the pharmaceutical field was initially proposed by [START_REF] Leuenberger | Granulation, new techniques[END_REF] for high shear wet granulation experiments. Leuenberger linked the motor power consumption profiles to the level of liquid saturation of the wet agglomerate and proposed a criterion for determining the amount of liquid required for granulating. Power consumption is also used to monitor and control the evolution of wet mixing in other application fields, such as concrete production [START_REF] Cazacliu | Concrete mixing kinetics by means of power measurement[END_REF]. However, power-consumption has a number of drawbacks (current drift with motor temperature, lag between current and voltage, use of over-sized motors and gearmotors which hinder current variations) that makes direct force or mixer torque measurements a reliable alternative for monitoring wet mixing processes [START_REF] Jørgensen | Comparison of torque measurements and near-infrared spectroscopy in characterization of a wet granulation process[END_REF][START_REF] Cazacliu | In-mixer measurements for describing mixture evolution during concrete mixing[END_REF][START_REF] Cavinato | Formulation design for optimal high-shear wet granulation using on-line torque measurements[END_REF][START_REF] Chitu | Wet granulation in laboratory scale high shear mixers: Effect of binder properties[END_REF]. The study of the force required to move an intruder through a wet granular mixture has therefore important industrial applications. Torque is a measure of the force needed to rotate the impeller which depends on the resistance of the wet mass against rotation of the impeller blades [START_REF] Bouwman | The relation between granule size, granule stickiness, and torque in the high-shear granulation process[END_REF]. Experiments monitoring the torque in high shear mixers therefore share some similiarity with the present case of the drag force applied by a wet granular system on a probe moving through it at constant velocity. The main difference however is the intensity of mixing, which is much larger in the high shear mixers, with significant centripetal accelerations of the wet material able to generate a toroidal flow pattern [START_REF] Cavinato | Scale-up effects on flow patterns in the high shear mixing of cohesive powders[END_REF] not present in the our pan mixer.

The behavior of wet granular materials can be characterized by a dimensionless number (Bond number), which compares body to capillary forces, and provides the framework to understand why increasing the particle size or decreasing the fluid surface tension yield a mechanical response similar to dry conditions. Bond number is defined as Bo = ρ p gd 2 p /γ, where ρ p is particle density, d p a characteristic particle diameter, γ the surface tension of the liquid and g the gravitational acceleration. Particle shape also plays an important role: for angular particles, capillary bridges, responsible of the attraction between grains, spread on the grain surface and are partially captured by surface roughness.

When focusing on the drag force on an object moving through a granular material, one may imagine that, if cohesion is introduced into the system, for example by wetting the grains, a cohesive stress scale is introduced which may change the drag scaling. It is important to stress that the addition of a liquid not only gives a tensile strength to the material, but also lubricates contacts, which results in a lower macroscopic coefficient of friction [START_REF] Santomaso | Collapse and runout of granular columns in pendular state[END_REF]. It is therefore interesting to study the effect of the addition of liquid (in this case, water) to the drag force on an intruder moving through a granular medium. This is the subject of the present work, which is, to our knowledge, the first study on drag forces on objects moving through wet granular media. It has to be emphasized that the experiments discussed in the paper were performed at a large, industrial scale, in the framework of a global project on the efficiency of industrial mixing of granular materials.

In the following, Section 2 presents the equipment, materials and procedures used; in Section 3 the raw results are briefly discussed, while Section 4 deals with an interpretation of the drag results by renormalization and comparison with shear cell analyses. Section 5 summarizes the main results and concludes the paper.

Materials & Methods

Materials

Experiments were performed on two sands, Fontainebleau sand (FS) and Teillé sand (TS) (Fig. 1). Fontainebleau sand is a fine siliceous sand with a quite narrow PSD. Teillé sand is a coarser and more heterogeneous river sand. Particle size distributions (PSD) for the two materials were determined by sieving and are displayed in Fig. 1. Material parameters are also summarized in Tab. 1. Particle shape was evaluated by computer-assisted microscopy using Fiji software [START_REF] Schindelin | Fiji: an open-source platform for biological-image analysis[END_REF]. The two sands resulted made of rounded particles with similar circularity (evaluated as 4πA/P 2 , where A is the projected area and P the perimeter). Before the tests, materials were air dried to reach a sufficiently low initial moisture level. The wetting liquid added in the subsequent phases was simply tap water.

Mixer

In this work, due to the need of mixing a certain amount of fluid with the granular medium, drag measurements were performed within a modified mixer. This was originally a tractor-mounted 1 m 3 pan mixer, characterized by 5 blades positioned at different distances from the central axis, and rotating in an annular space (see Figure 2). In our case, the mixer was connected to an asynchronous motor equipped with an inverter to easily vary the velocity. One of the blades was replaced by the force measurement system (see Section 2.3). With regard to the dimensions of the mixer, the inner and outer diameter of the mixing zone are respectively 500 mm and 1820 mm, and the maximum height attainable by the material from the bottom of the vessel is 450 mm. The angular velocity of the blades was varied in the range ω = 2.1 -21 rpm.

Force measurement technique

In order to evaluate the drag induced on a body moving through a granular medium, one of the blades of the mixer was replaced by a shaft at the end of which a spherical body was attached. The shaft was made of aluminium and had the length and diameter of 280 mm and 35 mm; the sphere (steel) had the two horizontal poles flattened, and had a diameter of 100 mm and an effective height of 90 mm. The probe axis was mounted at a distance of 620 mm from the center of rotation of the blades (which yields a velocity range of 0.14 to 1.4 m/s ). The height of the sphere with respect to the bottom of the mixer was kept constant, and corresponded to 50 mm (from the bottom part of the sphere), which was chosen based on a balance between the need of maximizing the highest possible level over the sphere and that of minimizing the interaction with the bottom boundary. The shaft was equipped with two pairs of strain gauges glued to its upper part. They allowed to monitor the bending of the shaft in both the streamwise and the transverse direction. The strain gauges were placed 230 mm above the sphere. The acquisition of the signals was performed within the mixer in a box located above the shaft which contained two conditioning cards (18bit Mantracourt DSCUSB OEM), a lightweight single board computer (Raspberry Pi Model 3) for data acquisition, and a battery for power supply. The force measurement system was therefore autonomous. Basic control of the acquisition procedure and data transfer were performed by means of an external laptop connected through a wireless link. Calibration was made by means of known masses attached to the sensing body. It is important to point out that the measurement set-up allows to quantify the moment acting on the shaft, not the force. Given that the material induces a distributed load on the sphere-and-shaft system which profile is not known a priori, it is not possible to deduce the total drag force rigorously from the bending moment. However, even if it is specific to the setup, this quantity can be used to investigate the relative effect of parameters such as velocity, material amount and water addition on the resistance provided by the granular medium on the motion of the intruder. For each configuration, the average moment was obtained by averaging the measured moment in time, while the standard deviation of the signal was taken as a measure of the uncertainty.

Level measurement

The level of material in the mixer, H, was determined after mixing by means of a rigid ruler placed 10 cm downstream the probe. This was chosen due to the need of measuring the level near the probe. The uncertainty estimation for this measurement given alongside with the results (±5 mm) comes from fluctuations with space and time of the height of the free surface. As discussed in the following sections two heights are important in our system: firstly, H, the height of the free surface with respect to the bottom of the mixer and, secondly, H * , the level above the half-sphere plane.

Shear cell

The flow properties of the two sands at different liquid contents were characterized through an annular shear cell (Powder Flow Tester, PFT, Brookfield). Samples of 500g of wet sand were prepared for each tests. Care was taken in order to standardize the procedure before and during the water addition step. Sand was dried in a ventilated oven for 4h at 100 • C and let cool at ambient temperature. The desired amount of tap water was then progressively added to the sample. A short stirring (∼3-4 min), carried out manually, was sufficient to ensure a proper distribution of the liquid in the sand. This procedure was based on visual inspection of the wet samples and on previous analogous experiments assisted by X-ray tomography on calcium carbonate particles [START_REF] Santomaso | Collapse and runout of granular columns in pendular state[END_REF].The sample was then immediately loaded in the shear cell to minimize water evaporation. The average solid fraction ν, internal friction angle φ and cohesion C were measured under application of consolidation stresses σ cons in the range 0.289 and 4.819 kPa. For an ideal material following the Mohr-Coulomb law (which is a reasonable assumption for sand), the cohesion C, the end-point consolidation stress σ cons and the internal friction angle φ are identified in Fig. 4 once the incipient yield locus of the material (τ vs. σ) is measured. In order to quantify the effect of water addition on the mechanical behavior of the granular materials, four parameters were given particular attention: the average solid fraction ν, the effective friction µ = tan φ, the normalized cohesion C/σ cons , and the ratio of the shear stress at incipient yield in wet and dry conditions τ wet /τ dry . These parameters were studied as a function of sand type, water amount and consolidation stress. In order to estimate the uncertainty of the measurements, 3 repetitions were made for both sands and for two water amounts (dry and w = 10%), and the standard deviation of the data was used as measure of the uncertainty.

Experimental procedure

In this work, we studied the effect of (i) the depth of immersion of the probe, (ii) the velocity and (iii) the amount of liquid added on the bending moment on the shaft of the probe resulting from the drag induced by the granular material. In a typical test, dry sand was loaded into the mixer at the beginning. Water was added in successive steps, up to saturation for the fine sand, and up to w ≈13% for the coarse sand (w is the mass ratio of the water to the sand). For the fine sand (FS), saturation was observed around w ≈30%, when a water film started to form at the bed surface. As it will be discussed in the following, for coarse sand we stopped at w ≈13% due to homogeneity issues. For each water amount, the material was first mixed with the liquid at the highest velocity until the mixture appeared to be homogeneous (which took nearly 5 minutes). Homogenization was made with the highest level. Different levels were obtained by partially voiding the vessel. The filling level was set by controlling the mass in the mixer. For the tests with water, three filling levels were studied, corresponding to 600, 800 and 1000 kg of sand. For FS, handling the material close to saturation revealed to be very complex, particularly with respect to partially voiding the mixer for varying the level. Above w ≈ 25%, it was therefore chosen to study only one filling level (corresponding to 600 kg of sand). Close to saturation, it was needed to seal the mixer discharge gate to avoid water leakage. In such cases it was therefore not possible to partially empty the mixer in order to vary the filling level. For this reason, handling FS material close to saturation was very difficult due to the large quantity of material that had to be removed manually from the mixer. Therefore, above w ≈ 25%, it was chosen to study only one filling level (corresponding to 600 kg of sand). For each level, 6 angular velocities of the mixer were tested, and the signal was acquired for more than one minute for each case. The level was manually measured every time.

As discussed in the following sections two heights are important in our system: firstly, H, the height of the free surface with respect to the bottom of the mixer and, secondly, H * , the level above the half-sphere plane. For each water addition the effective water content was measured at the beginning and the end of the tests, by drying and weighing material samples. After the tests for one water amount were completed, the material extracted for varying the level was loaded again into the mixer and new water was added to increase the water content for the further tests.

Results

Figure 5 displays the effect of the mass of sand and velocity on the moment measured by the probe for the materials in the dry state. It is important to note that while Fontainebleau sand was characterized by a negligible initial water content, the different process used for drying the Teillé sand ended with a water content around 0.5%. However, while not perfectly dry, it was verified that the sand was still free flowing at this value of water amount. Fig. 5 is given as a reference for the dry case. We can see that the moment is independent of velocity in the range considered, and it increases when the mass of the solids increases. FS is characterized by slightly larger values of the moment; this is not necessarily due to a direct effect of particle size. As it will be clearer in the following, this is probably more related to a slight difference in the solid fraction of the medium. Fig. 6 displays an example of the effect of the mass of sand and velocity on the moment measured in a wet configuration (around 4%). When adding water to the system, the difference between the two materials increases; rate-independence is also confirmed, as soon as the material remains unsaturated. For example, for water amounts above the saturation value (not shown here), the moment may decrease with velocity given that at high velocity the mixing action partially suspends the sand, while at low velocity sand sediments, and thus requires a higher effort to be displaced. For water amounts above 8%, Teillé sand is not well homogenized by the mixing process and the bed is formed by unstable aggregates. Above 13%, the coarse nature of the grains induces partial percolation of water which migrates towards the bottom of the vessel. For this reason, TS data stop at w ≈ 13%, while for FS it was possible to reach the saturation without observing heterogeneities in the mixture. However, it has to be noted that, for TS, the moment is quite independent of the velocity even in presence of such heterogeneities.

In Figure 7, the average moment measured at an intermediate velocity value is plotted with respect to the effective water content. The effect of the initial mass of sand appears to be similar between wet and dry materials. Then, water addition clearly increases the resistance to the motion of the probe by the material, and we can see how for TS (corresponding to the phenomenology described above), the moment drops down around w ≈ 10%, while for FS it increases up to w ≈ 20%, after which it decreases. Given that for FS, the decrease in the moment is close to saturation, it could be ascribed, for example, to lubrication of interparticle contacts by the surrounding fluid. On the other hand, for TS a sudden drop of the moment is obtained for lower values of w, which seems to be related to the developement of a fractured, unstable, bed which was highlighted previously.

Discussion

Phenomenology

As discussed in the introduction, previous literature has already focused on the force acting on an object moving through a dry granular material. To summarize, the main known features of dry granular drag are that the force is independent of velocity, independent (or weakly dependent) of particle size, but dependent on depth and on object size (and shape). In our case, the independence of drag on velocity is confirmed for dry as well as for wet cases, for all the water amounts tested. This suggests that the frictional, plastic nature of granular materials may remain relevant even for high water amounts. This idea is qualitatively supported by the increase of the measured moment with the initial mass of material, which appears to be similar for wet and dry cases (apart from TS for w > 10%, which will be discussed below). From Figure 7, it is clear that the addition of water into the system yields an increase in the resistance to the motion of the probe, for both TS and FS. It is tempting to directly ascribe this effect to the cohesion induced by water, by thinking that cohesion increases the stress which is needed to shear the material. However, as pointed out by [START_REF] Santomaso | Collapse and runout of granular columns in pendular state[END_REF], the addition of water has also a lubricating effect: if the cohesion increases, internal friction may decrease. Apart from these effects, due to the free surface nature of the flow configuration, another important phenomenon was observed, which is the dilation of the bed with the addition of water. This effect is shown in Figure 8, in which it is clear that cohesion stabilizes loose configurations, and the height of material in the vessel increases when water is added, reaching a value which can be up to 1.4 times the initial level. This increase is observed until a water amount around w = 10%, after which the level decreases. The addition of water implies, for the same solid fraction, an increase of the mass of the material above the probe, and so, the simple fact of adding a liquid could yield an increase in the measured moment due to the extra mass. The dilation induced by cohesion, though reducing the overall bulk density, acts in the same direction.

Given that the position of the probe is constant with respect to the bottom of the vessel, an increase of the level due to dilation means that a higher weight of material interferes with the probe (and this even if the av- erage solid fraction is lower) and will contribute to increasing the measured moment. Therefore four phenomena induced by the liquid addition may contribute to change the resistance to motion provided by the material: cohesion (due to capillary forces), dilatancy (induced by cohesion), extra liquid mass, and internal friction (modulated by contacts lubrication). In the following we try to analyse the data presented in Section 3, with the aim of understanding the effect of these phenomena on drag.

Scaling of the moment with H for dry materials

In order to evaluate the effect of dilatancy, the easiest way is to compare the moment measured in a wet configuration with an equivalent dry configuration characterized by the same level and solid fraction. While we expect from literature that the dependence of forces will be linear on solid fraction, the effect of the level is system-specific (its scaling depends on the probe size and shape). For this reason, a series of measurements on FS was devoted to characterize the dependence of the measured bending moment on the level, for the dry case. In Figure 9, in order to find a scaling, we plot the moment M rescaled by ρg, where ρ is the average bulk density, estimated from the mass and the measured level, as a function of the level (here taken with respect to the half plane of the spherical part of the probe, and noted as H * ). A fitting law can be extracted from these data, characterising the scaling of the moment on solid fraction and depth. Due to the discontinuity related to the passage from the cylinder to the sphere, we expect that a different fitting function should be used when the material is only partially covering the sphere (H * < 0.045 m in Figure 9). In this case, given that wet mixtures were always obtained for levels above the sphere, we fitted only the part of the curve above H * = 0.045 mm. ] sphe e shaft As it is visible in Fig. 9, our data is well fitted by the relation:

M ρg = f dry (H * ) = αH * 2 + βH * + δ, (1)
where α = 1.21 × 10 -1 m 2 , β = 1.55 × 10 -2 m 3 , δ = -9.98 × 10 -5 m 4 . While the fit was obtained on FS, it is surprising to see that the "dry" points of TS fall on the same curve; this may be related to the fact that both TS and FS are rounded siliceous sands, and this supports the idea of a similar macroscopic response. For this reason, we will use the same functional form for the two sands in the following. The form of Eq. 1 is consistent with the work of [START_REF] Albert | Slow drag in a granular medium[END_REF]Albert et al. ( , 2001a)). Note also, that the evolution of the force (and consequently the moment) versus the depth of the probe is far from being trivial due to the discrete nature of sand (Albert et al., 2001a).

Extra mass effects vs mechanical cohesion

From the previous discussion, it is clear that a mass increase in the mixer (due to an addition of solids or liquid), or a level raise due to capillarity-induced dilation, may augment the mass above the probe and therefore increase the resistance to motion of the intruder. It can be interesting to try to separate these "extra mass" effects from the cohesion itself, which, despite having the same micromechanical origin as dilatancy (i.e. capillary forces), can be seen as a separate effect whose importance can depend on the flow configuration. A way to separate the two phenomena is to study the ratio between the moment measured in wet and dry conditions, and to split it by means of the obtained scaling for dry materials:

M wet M dry = Mam Mcoh = ρ wet g f dry (H wet )
M dry M wet ρ wet g f dry (H wet ) ,

(2) where the scaling for dry materials is introduced to mimick the extra mass effect provided by an increase in the mass of the solids or the liquid, or by dilation, and where ρ wet is the average bulk density estimated from the level and the mass of materials, considering also the addition of water. Therefore, the first factor in the RHS of Eq. 2, Mam , can be used as a measure of the importance of the extra mass effects, while the second factor,

Mcoh , is meant to characterize separately the effect of mechanical cohesion. The results of the normalization is displayed in Figure 10, where error bars are computed by propagation of the uncertainties of level and moment measurements. The first observation is that at least for low values of the water content, Mcoh is close to 1, while Mam increases significantly with the water amount. This suggests that for low water contents, for both sands the main mechanism affecting drag is provided by the extra mass effects, and therefore the drag force is determined by the effective depth of immersion of the probe. For a given mass of solids, the strongest effect of the addition of water on the moment recorded by the probe seems to be related to the variation in level due to cohesioninduced dilation. However, if we look in more detail at the materialdependent behavior, the two sands employed behave differently. From Fig. 10a, it is clear that the contribution of extra mass effects is higher for TS than for FS, which means that the dilation induced by the liquid addition is more important for the coarse sand. This is clearly related to the variation in level after water addition, which is stronger for TS than for FS (Fig. 8). Also, the clear decrease of the level observed for the lowest filling of FS results in a decrease of Mam with the water content, for w > 10%. With respect to Fig. 10b, for Teillé sand, the addition of water corresponds to a slight decrease in Mcoh , which may be the signature of lubrication. Teillé sand is the coarser sand used in this work, so it is reasonable to think that lubrication may interfere with cohesion for it. For w ≈ 10%, where the unrescaled data (see Fig. 7) displayed an abrupt decrease of the moment, Mcoh decreases as well, and the points corresponding to the different initial levels do not superpose, which means that the scaling on depth based on dry data does not hold. As it was partially discussed previously, for water amounts larger than 10%, TS starts to be very heterogeneous and forms lumps which are not easily broken by the mixing process. The material is very dilated and fractured, and the probe easily penetrates through this fractured bed. This is probably the phenomenon originating the progressive failure of the scaling obtained for dry data, when w ≥ 10%. For Fontainebleau sand, the picture is clearly different. When increasing the water amount, the effect of the initial mass of the solids becomes increasingly important for this fine sand. Mcoh is larger the lower the mass of the solids, and the difference between the curves increases. This means that the scaling of the moment with depth obtained for dry materials progressively fails as cohesion becomes more important. The order of the curves in Fig. 10b is also important. For the highest initial level, Mcoh is nearly constant with respect to w, and slightly lower than 1. This means that the moment measured by the probe is close to the one corresponding to a dry bed. When the initial level is lower, the moment exceeds the dry reference value by up to 50%.

To summarize, Fig. 10a and Fig. 10b show that the socalled extra mass effects were generally very important in our experiments, apart from the lowest filling level of the fine sand, where for high water amounts mechanical cohesion played a major role. In other words, while the effect of cohesion is important only close to the free surface, the effect of water addition is felt also by the deep layers through the variation in level induced by dilation. This behavior is due to the fact that the importance of cohesion must be considered relatively to body forces. In our case, body forces increase as we move deep inside the bed; moreover, due to the shape of the probe, the frictional nature of the drag and the nature of the quantity measured, the total moment measured does not represent equally all the layers of material touching the immersed part of the probe, but gives more importance to the deepest layers. It is therefore reasonable to think that there is a position in the bed below which cohesive forces are overcome by the mass of the material: this yields the fact that for 1000 kg of sand, the moment measured follows the scaling obtained for the dry material. If the probe is sensing less deeper layers, it is allowed to register the effect of cohesion. We can ascribe the excess of moment with respect to the reference scaling to an increase in cohesion. This marks an important difference with respect to Teillé sand, for which the effect of cohesion was not visible. This difference is clearly related to the particle size, which drives the importance of capillary forces with respect to body forces. As it was stated in the Introduction, a dimensionless number called Bond number can be used to quantify the sensitivity of a granular material to capillary forces.

In our case, Bo, computed with the median diameter, is respectively 0.17 and 0.022 for TS and FS sands. Even if Bond number does not consider the contribution of the confinement stress, and is therefore more relevant close to the bed surface, it gives an indirect proof that for both sands, capillarity can influence the material behavior, but that FS is more sensitive. In addition, one should consider that TS has a broader size distribution than FS, and given that the finest particles in the distribution (d < d 10 ) have particle sizes comparable to FS, this may increase the sensitivity of FS to the water addition.

Shear cell results

In order to understand the origin of the different mechanical behaviors observed with the two types of sand, we discuss in this section the results obtained with the shear cell described in Sec. 2.5. In particular, in Fig. 11 we plot different parameters obtained from the shear cell as a function of water content for the two types of sand and different consolidation stresses. An increase of the consolidation stress mimics an increase of the level in the mixer: in the latter, the material is sheared under its own load and so we can assume that the consolidation stress increases with depth. Given that the probe is fixed, this means that the consolidation stress near the probe end increases with the level of the material. The three values of consolidation stress used in the shear cell are of the same order of magnitude of the hydrostatic pressure near the end of the probe for the different levels in the mixer. However, as it was highlighted in the previous sections, the level in the mixer varied significantly upon the addition of water, while in the shear cell the consolidation stress was kept the same for the sake of simplicity while varying the water content. The shear cell results are not meant to provide a punctual comparison, but were done to obtain an independent measurement of the effect of water addition on the forces needed to shear the material.

The average solid fraction of the materials at incipient yield is displayed in Fig. 11a. TS has generally a higher solid fraction than FS, probably because of the highest polydispersity. This feature causes the formation a more compact packing. It is interesting to note that ν decreases when adding water up to w = 5% and then remains nearly stable. The material loosens when adding water, and this reflects the behavior observed in the mixer (Fig. 8). However, the minimum of solid fraction is obtained for lower values of the water amount than the maximum of the level in the mixer, and above 10% the level in the mixer decreases with w, which is not mirrored by the solid fraction in the cell. This is probably due to the different geometry: in the shear cell, the material can partially expand due to the stresscontrolled configuration, nervertheless it remains a fully confined setup. On the other hand, the mixer is a free surface configuration in which dilation is probably less disturbed. In addition, the bulk solid fraction can be assumed as uniform in the shear cell while it will probably depend on depth in the mixer.

In Fig. 11b, the effective friction µ (the slope of the incipient yield locus) is displayed. TS is characterized by higher values of the friction coefficient, which can be related to a higher level of interlocking between the particles due to the higher solid fraction. It is evident that the addition of water decreases µ at least until w < 5%. An increase of µ is then apparent for w > 10%. With respect to effective friction, the two types of sand respond slightly differently to the addition of water, with the finer (FS) displaying a larger variation of µ. A decrease of friction with w for w < 5% was already discussed by [START_REF] Santomaso | Collapse and runout of granular columns in pendular state[END_REF], who postulated that the micromechanical origin of such phenomenon could be related to the lubrication of interparticle contacts by the wetting liquid. While this is confirmed by the present data, we also observe a non-monotonic variation of µ for higher water contents which would require further investigation.

Fig. 11c displays the evolution of cohesion parameter C, normalized by the consolidation stress σ cons , with the water content. As one may expect, cohesion generally increases when the liquid is added to the granular material, and this increase is nonlinear, with a larger variation between w = 0 and 2%. The coarse sand appears to be less sensitive to cohesion (C/σ cons for FS is roughly 4/3 of C/σ cons for TS), which is the expected effect of particle size, and the importance of cohesion decreases when increasing the consolidation stress. This behavior is analogous to what was observed in the mixer (Fig. 10), where the effect of cohesion on the measured moment decreases with the filling level of the mixer. Together, Figs. 11b and11c show the importance of both internal friction and cohesion phenomena for the two types of sand studied in this work.

In order to evaluate the effect of water addition on the global behavior of the mixture, we plot in Fig. 11d, the ratio of the shear stress at incipient yield in wet conditions to the one in the dry case. This quantity is related to the normalized moment shown in Fig. 10, because it compares the effort needed to shear the material in wet and dry conditions. From the figure it is evident that τ wet /τ dry increases with the water addition: this means that the decrease in the slope of the yield locus is always compensated and exceeded by an increase in cohesion. This is different from drag data, in which the coarser TS displayed a decrease of the normalized moment with w, and therefore a global balance in favour of internal friction, modulated by contact lubrication. The opposite variation of τ wet /τ dry in the shear cell and M in the mixer have probably to be ascribed to the different flow configuration, and particularly to the more pronounced confinement provided by the shear cell, as it was discussed above with respect to solid fraction. However, apart from the variation, it has to be noted that τ wet /τ dry for TS seems to be insensitive to the consolidation stress; this peculiar behavior finds a correspondence in Fig. 10 and can be seen as a signature of the reduced importance of cohesion for the coarser sand. On the other hand, for the finer sand (FS), the effect of water visible in Fig. 11d is stronger and the curves are sensitive to the consolidation stress σ cons . Even in this case, the shear cell provides an overestimation of the effect of water with respect to the results obtained in the mixer, but it is important to notice that the order of the curves is the same that in Fig. 10, and the curve corresponding to the highest consolidation stress is close to the curves obtained for TS. The origin of the increase in the yield stress upon the addition of water in conditions for which we would normally suppose that the liquid is negligible (coarse particles, high consolidation stress) needs further investigation. Still, shear cell results confirm that the drag force on an object moving through a wet granular material is a complex phenomenon which involves the effect of cohesion, dilatancy and internal friction. These three phenomena may all be relevant at the same time, as the present results suggest.

Conclusions

The present paper reports an experimental study of the drag force applied on an object moving through a wet granular medium with the aim to compare our results to those obtained with the well-known dry case. We used both fine (Fontainebleau sand, d 50 ≈ 250 µm) and coarse (Teillé sand, d 50 ≈ 690 µm) sands. First, the moment of the drag force is found to be independent on the velocity of the probe and to increase with the amount of granular medium, similar to what has been observed for dry granular materials. Yet, strong differences between the dry and wet cases exist, consequences of (i) cohesion due to capillary forces, (ii) pos-sible contact lubrication and (iii) dilation induced by cohesion. We have shown that, among the three aforementioned mechanisms, the effect of dilation is the most important, at least in our experimental configuration. The effect of cohesion at the system scale depends on the relative influence of the cohesive forces with respect to the grains' weight. This implies that, for equal material densities, a coarse sand is less sensitive to cohesion than a fine one. Consequently coarse grains may experience a significant contact lubrication. At important water contents, the granular system may become strongly heterogeneous. We also carried out shear cell experiments to obtain an independent measurement of the effect of water content on the stress required to shear the granular media used in the industry-scale mixer. Our results confirm that both contact lubrication (which prevails at low water content) and cohesion influence the mechanical response of wet granular systems. They also shed light on the importance of the geometry of the system which by its characteristics (e.g. the existence of a free surface) influences the dilation of the system and thus its behavior.

Figure 1 :

 1 Figure 1: Microscope pictures of (a) Fontainebleau and (b) Teillé sand. (c) The cumulative particle size distributions of the two types of sand show that the former is coarser and its size distribution broader.

Figure 2 :

 2 Figure 2: Internal view of the mixer adapted for this study. The four blades placed at different radial positions allow to mix the sand and the water. The upper part of the force measurement system is visible in the foreground.

Figure 3 :

 3 Figure 3: Top-view sketch of the mixer configuration with the probe location and detail of the system used for measuring the drag induced by the granular material. The heights H and H * are respectively the position of the free surface with respect to the the bottom of the mixer and the level above the half-sphere plane H * .

Figure 4 :

 4 Figure 4: Mohr-Coulomb representation of the material parameters obtained with the shear cell.

Figure 5 :

 5 Figure5: Average moment recorded by the probe as a function of probe velocity and amount of sand in the mixer for the "dry" materials, for Fontainebleau sand (dashed lines, w < 0.01%) and Teillé sand (solid lines, w = 0.46%).

Figure 6 :

 6 Figure6: Average moment recorded by the probe as a function of probe velocity and amount of sand in the mixer for wet materials with a water amount around 4%, for Fontainebleau (dashed lines) and Teillé sand (solid lines).

Figure 7 :

 7 Figure7: Average moment recorded by the probe as a function of water amount and amount of sand in the mixer, for Fontainebleau (dashed lines) and Teillé sand (solid lines). The data correspond to a probe velocity of 0.68m/s.

Figure 8 :

 8 Figure 8: Evolution of the height of material in the mixer as a function of water amount and amount of sand in the mixer. Dashed lines correspond to FS and solid lines to TS. Water induces dilatancy up to w = 10% due to the free surface nature of the flow configuration.

Figure 9 :

 9 Figure 9: Average moment recorded by the probe for dry materials (squares: FS, circles: TS), rescaled by ρg, as a function of the level above the half-sphere plane H * .

Figure 10 :

 10 Figure 10: Contribution of (a) the extra mass effect induced by addition of solids, liquid and by dilation and (b) mechanical cohesion on the average moment recorded by the probe as a function of water amount and mass of sand in the mixer, for Fontainebleau (dashed lines) and Teillé sand (solid lines). Same data as in Fig. 7.

Figure 11 :

 11 Figure 11: Effect of sand type, water content and consolidation stress on bulk material parameters estimated with a shear cell : (a) average solid fraction ν at incipient yield, (b) effective friction coefficient µ = tan φ, (c) normalized cohesion C/σ cons , and (d) ratio of the shear stress at incipient yield in wet and dry conditions τ wet /τ dry .

Table 1 :

 1 Characteristic properties of the sands used in this work.

Acknowledgments

We thank D. Hamon, A. Cothenet, F. Buisson for the help with the sample preparation and mixer operation and J. M. Paul for the development of the sensor.


Arthur Le Rhun
email: arthur.le-rhun@inria.fr

Frédéric Bonnans
email: frederic.bonnans@inria.fr

Giovanni De Nunzio
email: giovanni.de-nunzio@ifpen.fr

Thomas Leroy
email: thomas.leroy@ifpen.fr

Pierre Martinon
email: pierre.martinon@inria.fr

A stochastic data-based traffic model applied to vehicles energy consumption estimation

des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

A stochastic data-based traffic model applied to vehicles energy consumption estimation Arthur Le Rhun, Frédéric Bonnans, Giovanni De Nunzio, Thomas Leroy, Pierre Martinon

Abstract-A new approach to estimate traffic energy consumption via traffic data aggregation in (speed,acceleration) probability distributions is proposed. The aggregation is done on each segment composing the road network. In order to reduce data occupancy, clustering techniques are used to obtain meaningful classes of traffic conditions. Different times of the day with similar speed patterns and traffic behavior are thus grouped together in a single cluster. Different energy consumption models based on the aggregated data are proposed to estimate the energy consumption of the vehicles in the road network. For validation purposes, a microscopic traffic simulator is used to generate the data and compare the estimated energy consumption to the measured one. A thorough sensitivity analysis with respect to the parameters of the proposed method (i.e. number of clusters, size of the distributions support, etc.) is also conducted in simulation. Finally, a real-life scenario using floating car data is analyzed to evaluate the applicability and the robustness of the proposed method.

Index Terms-Traffic modeling, Clustering, Energy Consumption

I. INTRODUCTION

I N 2015, according to data from the European Environment Agency, road transportation contributed to 21% of total EU-28 greenhouse gas emissions. In order to meet the long-term emissions reduction target, emissions from transportation need to fall by more than two thirds by 2050 [START_REF]European Strategies -White Paper[END_REF]. These emissions are essentially a function of the vehicle propulsion technology and the driving style [START_REF] Evans | Multivariate Analysis of Traffic Factors Related to Fuel Consumption in Urban Driving[END_REF].

Estimating energy consumption of the vehicles is a great challenge in the objective of improving global transportation efficiency, since this information is used in energy management, eco-routing, eco-driving, traffic management, ... Traffic congestion has a major impact on the driving behavior, and thus plays a key role in the level of fuel consumption [START_REF] Vlieger | Environmental effects of driving behaviour and congestion related to passenger cars[END_REF].

Therefore, accurate predictions of vehicles energy consumption must take traffic conditions into account. To perform this objective, faithful modeling of traffic behavior is of primary importance. Energy-oriented modeling approaches can be divided in two main categories. On the one hand, several mathematical traffic models are available nowadays, see for instance [START_REF] Treiber | Traffic flow dynamics[END_REF]. Such models typically depict the reality either from a macroscopic point of view, based on the road vehicular density [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF], or from a microscopic perspective, based on the description of the instantaneous behavior of each vehicle [START_REF] Panwai | Comparative evaluation of microscopic carfollowing behavior[END_REF]. Both approaches have limitations in providing an accurate energy consumption estimation. Macroscopic models typically provide average traffic speeds to compute energy consumption [START_REF] Zegeye | Integrated Macroscopic Traffic Flow, Emission, and Fuel Consumption Model for Control Purposes[END_REF], thus neglecting the impact of speed fluctuations due to congestion. Higher precision of the energy consumption estimation could only be obtained at the expense of a denser discretization of the road network, therefore compromising scalability. Microscopic models could achieve precise energy consumption estimation, but they require a significant calibration and validation effort. Also, the computational burden and the amount of collected data grows rapidly with the size of the network, therefore these models are more suitable for off-line use. On the other hand, data-based models rely on collected traffic information to estimate traffic behavior and energy consumption. Instantaneous models are able to precisely estimate energy consumption by using large amounts of data, generally the measured driving profile of each vehicle. To tackle this drawback, aggregated models use the average value of the measured speed profiles to compute energy consumption, but they suffer from the same accuracy problems previously discussed for the macroscopic traffic models. Furthermore, the data sparsity and availability is an issue [START_REF] Guo | EcoMark 2.0: Empowering Eco-Routing with Vehicular Environmental Models and Actual Vehicle Fuel Consumption Data[END_REF]. Other approaches try to solve the problem of the data sparsity by simply classifying road segments by category (e.g. urban, arterial, freeway, etc.), in order to associate each category with a typical energy use. This type of models may lead to inaccuracy in energy consumption estimation, as road segments belonging to the same category may show very different traffic patterns [START_REF] Ericsson | Variability in urban driving patterns[END_REF].

In this work, a new way to represent traffic behavior on large road networks is proposed. The objective of this model is to accurately depict the effect of traffic conditions on the vehicles energy consumption in each road segment. The key idea is to use a statistical approach based on vehicle speed and acceleration data, measured from real vehicles. In particular, the entire observation time during which speed and acceleration data are collected is subdivided into time-frames. During each time-frame a (speed,acceleration) distribution is generated, and such a distribution is then used as an input for an energy consumption model to estimate the traffic energy consumption on the analyzed road segment during the specific time-frame. Therefore, each road segment is defined by its own collection of (speed,acceleration) distributions.

The proposed model also includes the possibility of reducing the dimension of the traffic data and increasing scalability, by applying clustering techniques to the probability distributions of each road segment. For instance, the different distributions representing traffic in one road segment over different hours of the day may be aggregated in clusters modeling only significant traffic conditions (e.g. peak, offpeak, etc.).

The paper is organized as follows. The clustering technique and the proposed energy consumption model are presented in Section II. The traffic data collection and the model validation procedure are discussed in Section III. Using traffic data from the simulator SUMO, Section IV illustrates the method first on a single road segment, then for a larger set of segments. Finally, section V presents an application of the method to actual traffic data, as well as a comparison to a more basic approach, using only mean speeds.

II. PROPOSED METHOD

A. Road segments

In the following, the road network is assumed to be subdivided into a collection of segments. A segment is a portion of road with homogeneous topographic characteristics. Segments are typically delimited by network elements such as traffic lights, crossroads, or roundabouts. With a traffic simulator such as SUMO we can retrieve the segments from the simulator model, see IV-A. When working on actual traffic data, segments can be obtained from Geographic Information Systems (GIS) data (e.g. OpenStreetMap, HERE Maps, etc.). Segments length typically range between a few meters and a few hundreds meters. In the following we denote S the set of road segments, of size N S .

B. Probability distribution for (speed,acceleration)

The consumption of an engine depends on its operating point, which can be determined by the speed and acceleration of the vehicle. The key point in our method is to assume that an accurate (with respect to energy consumption) description of the traffic can be derived from the probability distributions of measured speeds v and accelerations a for each road segment. Note that these distributions do not retain the temporality of the speed profiles.

segments and time-frames, we obtain a total of N S N T discrete distributions of support size N V N A . Fig. 1 shows an example of such a distribution.

C. K-means clustering with strong patterns

For each segment s and time-frame t i,s we have a probability distribution µ i,s . When considering large road networks with several thousands of segments, the data size grows rapidly to several gigabytes. A natural idea is to reduce the information through clustering techniques applied to these distributions. Then, for each cluster j, we can take its barycenter μj as representative of the traffic conditions for all the segment/time-frame pairs (s, t i,s ) that belong to this cluster. This way, we only have to store the distributions corresponding to the barycenters of the clusters.

Since the elements to be clustered are probability distributions, we use the 1-Wasserstein distance ( [START_REF] Vasershtein | Markov processes over denumerable products of spaces describing large system of automata[END_REF], [START_REF] Dobrushin | Prescribing a System of Random Variables by Conditional Distributions[END_REF]). This distance based on optimal transport theory tends to preserve the geometrical aspects (shapes) of the distributions. We recall the definition of the Wasserstein distance:

W 1 (µ, ν) = min π∈Π(µ,ν) d • π (1)
where Π(µ, ν) is the set of transportation plans from µ to ν, i.e. the set of nonnegative matrices π of size

(|N V | × |N A |) 2
with marginals π 1 = µ and π1 = ν. As speed and acceleration have comparable magnitude orders when expressed in SI units, we simply define the displacement cost d by the Euclidean norm between the origin and destination points.

The notion of barycenter for distributions has been extended to the Wasserstein distance in [START_REF] Agueh | Barycenters in the Wasserstein space[END_REF], where the barycenter is defined as the distribution that minimizes the sum of the squared Wasserstein distances to all distributions of the set. For the practical aspects we follow [START_REF] Benamou | Iterative Bregman projections for regularized transportation problems[END_REF] in which the barycenter minimizes the sum of W distances, with W an entropic regularization of the Wasserstein distance, computed by a Sinkhorn-type algorithm.

We use the K-means algorithm [START_REF] Macqueen | Some methods for classification and analysis of multivariate observations[END_REF] to compute the clusters, and note k the (fixed) number of clusters. To reduce the sensitivity of the K-means algorithm to the initial guess, we use the strong patterns method [START_REF] Lefebvre | Dynamic clustering and strong patterns recognition: new tools in automatic classification[END_REF]. Strong patterns are subsets whose elements always end up in the same cluster regardless of the K-means starting point. In practice, we run a first batch of K-means with random initializations to identify the strong patterns, and then perform a final K-means which is initialized by taking one element in each of the k largest strong patterns. The principle of the K-means is recalled in Algorithm 1.

Algorithm 1: K-means algorithm

Input : Distributions to be clustered µ 1

D. Computing energy consumption

The energy consumption computed is the energy at the wheel, neglecting the losses due to the powertrain. The instantaneous power at the wheel is denoted as a general function P (v, a), which can be for instance of the form presented in [START_REF] Lighthill | On kinematic waves. II. A theory of traffic flow on long crowded roads[END_REF]. The ultimate objective of the proposed method is to estimate energy consumption by using only the information extracted from the (v, a) probability distributions µ. More accurately, we seek to obtain the energy consumption from the barycenter μ of the cluster containing µ. In the following we introduce two methods to compute the consumption of a generic vehicle passing through a segment. These will be referred to as "Average Consumption Method" and "Memoryless Sampling Method".

1) Average Consumption method:

The first idea is to use the average power P (in the probabilistic sense)

P ( μt ) = V×A μt (v, a)P (v, a) (2)
where μt is the barycenter of the cluster containing the current segment at time t. The barycenter may indeed change if t crosses different time-frames while the vehicle is on the segment. This average power is integrated over the time interval [t i , t f ] spent by the vehicle on the segment, thus

C Avg = t f ti P ( μt )dt (3)
Note that knowledge of the time interval is required in this method, in addition to the (speed,acceleration) distributions.

Indeed, using here the average time would give identical consumption for every vehicle. Therefore, we need some more statistically significant time information in order to capture the deviation of the consumption distribution. A resulting drawback of this method is that a faster vehicle has a shorter travel time and thus a lower energy consumption, which may seem unrealistic.

2) Memoryless Sampling method: In the second method the energy consumption is still obtained by integrating the instantaneous power, but we do not use the average power. Instead, we implement the idea that the vehicle must follow the traffic at every time, in a statistical sense. More precisely, its speed and acceleration should follow the probability distribution of the barycenter μ of the cluster for the current pair (s, t). Another difference is that the integration is performed over the segment length L s instead of travel time.

So the Memoryless Sampling method generates a sequence of (v n , a n ), independent samples according to the probability distribution μ. Setting a time step δt, we use this sequence to integrate both the traveled distance and the instantaneous power. We assume δt = 1s, in order to have the same order of magnitude as the reaction time of a driver. We stop the generation of (v n , a n ) when the vehicle reaches the end of the segment2 . Since the distance will be covered in a finite time, we obtain the finite set of samples (v n , a n ) n=1,...,n f and the consumption writes as

C M SM = n f n=1 P (v n , a n )δt (4)

III. VALIDATION APPROACH A. Power and Reference energy consumption

In order to assess the accuracy of our two methods (Average and Memoryless Sampling), both based on the statistical representation of traffic, we introduce a "reference" energy consumption for the sake of comparison. First we choose a simple equation for the instantaneous power P (v, a), neglecting all road slope effects:

P (v, a) = (ma + a 2 v 2 + a 1 v + a 0 )v ( 5
)
where m is the vehicle mass and a 0 , a 1 , a 2 define a vehicledependent polynomial approximation of the road-load force.

We use the numerical values m = 1400kg, a 0 = 185.4, a 1 = 0, a 2 = 0.3, corresponding to a passenger vehicle [START_REF] De Nunzio | Energy-optimal driving range prediction for electric vehicles[END_REF]. We would like to point out that energy consumption models are a wide topic, and we refer interested readers for instance to [START_REF] Hayes | Simplified electric vehicle powertrain model for range and energy consumption based on epa coast-down parameters and test validation by argonne national lab data on the nissan leaf[END_REF].

In our case we essentially want to compare the consumptions from raw traffic data and statistically processed data, so the key point is to use the same power expression for all consumption formulas. In all the following, these energy consumptions C are computed by integrating3 the instantaneous power P along the (speed,acceleration) profiles, without taking into account any regenerative braking, namely

C = t f ti max(0, P (v(t), a(t)))dt (6)
In order to compare the different methods, we compute the distribution of the energy consumptions for a given (segment,timeframe): -the Reference energy distribution C Ref is obtained by plugging into (5) the recorded speed and acceleration (v(t), a(t)) of all vehicles passing through the segment during the timeframe, according to the traffic simulator.

-for the Average and Memoryless Sampling methods, we recall that the consumption distributions C Avg , C M SM are obtained using the speed/acceleration probability laws μ from the barycenters, see ( 3) and (4).

B. Indicators

Since our aim is to compare distributions of energy consumptions, we study several indicators.

• Mean and Standard deviation errors are classical indicators. We compute these relative errors as follows, with C method denoting either the Average consumption or the Memoryless Sampling consumption:

ε mean (s) = C method (s)-C ref (s) C ref (s) ε σ (s) = σ C method (s)-σ C ref (s) σ C ref (s)
• Kullback-Leibler divergence [START_REF] Kullback | On information and sufficiency[END_REF], also called 'relative entropy', is a particular case of ϕ-divergence. KL divergence can be used to measure distances between two probability distributions P and Q, however it is not a metric (no triangular inequality or symmetry). Another drawback is that it cannot be computed for instance when the probability of the model q is 0 while the probability of the reference p is not.

KL(P |Q) = i p i log( p i q i ) (7)
Instead, we use the Jensen-Shannon divergence, which is a symmetrized version of KL divergence, sometimes referred to as 'total divergence to the average' [START_REF] Dagan | Similarity-based Methods for Word Sense Disambiguation[END_REF]. Note that the square root of the JS divergence is a metric called JS distance [START_REF] Endres | A new metric for probability distributions[END_REF], [START_REF] Österreicher | A new class of metric divergences on probability spaces and its applicability in statistics[END_REF]. Fig. 2 illustrates the JS divergence on a Gaussian with noised parameters.

JS(P |Q) = 1 2 KL(P |M ) + 1 2 KL(Q|M ) (8)
with

M = 1 2 P + 1 2 Q (9)
Here, we have P = C method and Q = C Ref .

IV. RESULTS WITH SIMULATED TRAFFIC DATA A. Traffic Data from simulation

We illustrate our approach with data obtained from the traffic simulator SUMO [START_REF] Krajzewicz | Recent Development and Applications of SUMO -Simulation of Urban MObility[END_REF]. The simulation runs the scenario LUST [START_REF] Codecá | Luxembourg SUMO Traffic (LuST) Scenario: Traffic Demand Evaluation[END_REF], which models a 24h traffic in the city of Luxembourg. The time step is set to one second, i.e. 1 Hz sampling frequency for the variables of the vehicles. For the road segments we take the subdivision from the scenario, consisting of roughly 24000 elements. We aggregate contiguous road lanes together to obtain more data per segment, and end up with 18322 segments. The length of these road segments is shown on Fig. 3.

From the simulation raw data (5.7GB), we extract the segment, speed, and acceleration of each vehicle in the network at all time steps. We aggregate the records for a fixed time-frame length ∆t s,i (in practice we use a constant frame length ∆t for all segments). This is done both to gather sufficient data on the segments, and to decrease the number of traffic distributions for the clustering phase. For instance a time-frame of 1h will give 24 distributions per segment, while 10 min gives 144 distributions per segment. In the following all codes are written in Python/NumPy and run on a standard desktop computer.

B. Numerical results on one segment

To begin with, we compare the different consumptions on a single road segment. The objects compared are therefore the consumptions distributions of all the vehicles that went through the segment during each time-frame. Results are shown as the cumulative distribution function of the consumptions for all time-frames.

We analyze in particular the influence of the (speed,acceleration) discretization, the clustering, and the choice of time-frame duration. Unless specified otherwise, the distributions are shown for a discretization N V = N A = 10, a time-frame ∆t = 10min, taking the full set of distributions without clustering.

1) Influence of the (speed,acceleration) discretization: First, Fig. 4 shows the consumption distribution for both the Average method and Memoryless Sampling method. We test the discretizations N V = N A = 10, 20 and 30 and compare to the reference consumption. We see that for the Memoryless Sampling method: i) the general shape of the distribution is similar to the reference and ii) finer (v, a) discretizations give distributions closer to the reference. On the other hand, for the Average method: i) we observe some linearization of the consumption and ii) the effect of discretization is much less significant. 2) Influence of the time-frame: Next we study the effect of the length of the time-frame ∆t, which is the time interval over which we aggregate the vehicles data. Longer time-frames may cause some over-averaging and loss of specific traffic information. On the other hand, shorter time-frames may lead to insufficient vehicle data (for statistical relevance), and also increase the number of (v, a) distributions to handle. Fig. 5 shows the energy consumptions obtained for time-frames of 5s, 1min and 10min. On this segment, for both the Average and the Memoryless Sampling methods, the influence of the time-frame duration seems rather small.

3) Influence of the clustering: To conclude this first batch of results, we examine the information loss due to the clustering stage. The consumptions obtained using only the barycenters (green and cyan curves in Fig. 4) are almost identical to the ones obtained without the clustering stage. This indicates that the data reduction performed by the clustering comes with a negligible loss of information.

C. Numerical results on 500 segments

Now, we consider the entire road network of the city of Luxembourg, as implemented in LUST, with 18000 road segments. In order to study the influence of the different parameters while keeping reasonable computational times, we pick a test set of 500 random segments, for which we perform the clustering and compute the energy consumption with the Average and Memoryless methods. For each segment, we use the indicators defined in III-B to compare the computed energy consumptions to the Reference consumption. We discuss the relevance of the models based on the distributions of these indicators on the set of segments. More precisely, we investigate the influence of the discretization of speed and acceleration, and the number of clusters. Unless specified otherwise, the simulations use a speed and acceleration discretizations of 20 steps, a time-frame of 10 minutes, with 2 clusters.

1) Influence of speed discretization: The speed discretization has a direct influence over the barycenters computed by the K-means, since it changes the support of the distribution obtained. We expect a finer discretization to give computed consumptions (Average and Memoryless Sampling methods) closer to the Reference ones. The obvious drawbacks are an increased cost of the barycenter computation and size of the distributions. We test N V = 5, 10, 20 steps for the discretization, the speed interval being [0, 20] in m/s.

• Relative mean and standard deviation errors

We begin with the distribution over the 500 test segments of the mean and standard error (both relative). Fig. 6a and6b (upper graphs) show the errors between the Reference energy consumption and the Memoryless Sampling method. The mean and standard error both appear to be reasonably well centered around 0. We also observe that finer discretizations of the speed clearly improve the standard error, possibly due to a better reconstruction of the travel times. On the other hand the mean error is shifted towards positive values, and is indeed smaller for 10 steps than 20.

Fig. 6a and 6b (lower graphs) compare the Average method consumptions to the Reference. Here the mean error is almost always negative, and the standard error is also negative for 80% of the segments. This strong unbalance towards the negative indicates that the Average method tends to underestimate the consumption. Increasing the speed discretization reduces the mean error, but does not really improve the standard error. A possible explanation is the fact that the travel time for each vehicle is taken from the simulation, thus the speed discretization has no effect on it.

• Jensen-Shannon Divergence Fig. 6c shows the distribution over the 500 segments of the Jensen-Shannon divergence. Upper graph is for the Memoryless Sampling versus Reference, and lower graph is for Average method versus Reference. We observe that the JS divergence is much smaller overall for the Memoryless Sampling than the Average method, with distributions more concentrated towards zero. For both methods, increasing the speed discretizations from 5 to 10 steps improves the JS divergence, while 20 steps yield very little additional gain.

2) Influence of acceleration discretization: Now we study the discretization of the acceleration. As for the speed, this parameter influences the support of the (v, a) distributions, and therefore the K-means clustering. We want to know if finer discretizations of a give more accurate energy consumptions for the Average and Memoryless Sampling methods. We test N A = 5, 10, 20 steps for the discretization, the acceleration interval being [-4.5, 4.5] in ms -2 . • Relative mean and standard deviation errors Fig. 7a shows the mean error distribution for the two methods, with different discretizations of a. We observe that for both the Memoryless and Average methods, going from 5 to 10 steps gives a significant improvement, while 20 steps is similar to 10. With a sufficient discretization, the mean error is extremely good for the Average method. The Memoryless method, on the other hand, tends to slightly overestimate the consumption.

Fig. 7b shows the standard error distributions. As observed for the speed discretization, error for the Memoryless is well balanced while Average method has mostly negative standard errors. Finer discretizations of a seem to give no improvement for the standard error. This may be due to the fact that acceleration has no influence on the travel time for either method, unlike speed which is used to reconstruct the travel times in the Memoryless method.

• Jensen-Shannon Divergence Fig. 7c shows the distribution of the Jensen-Shannon divergence when varying the discretization of a. Like in the speed discretization study, we observe that the JS divergence is much smaller overall for the Memoryless Sampling method. Once again, for both methods increasing from 5 to 10 steps improves the indicator, while 20 steps give no additional benefit.

3) Influence of the number of clusters: Our traffic model uses clustering techniques to reduce the size of the traffic data, while retaining the useful information. One would expect some kind of trade-off, where using a larger number of clusters would keep more information at the expense of data size.

For this particular segment, simulations with k = 2, 3, 4 clusters give almost identical consumption distributions. Note that these simulations were run for a time-frame ∆t = 1h, due to the increased computational cost for 3 and 4 clusters. It remains to be seen whether a smaller time-frame would benefit more from a higher number of clusters. Also, the ideal number of clusters is most likely segment-dependent. Experiments with real data in Sec. V give more insight into the impact of the number of clusters.

4) Summary and performance analysis: Generally speaking, increasing the speed, acceleration, time discretization and the number of clusters will improve the accuracy of the whole approach. However, this comes at a cost, and some parameters have more impact than others. The detailed influence of the method parameters is as follows: -speed/acceleration discretization: greater N V improves the Memoryless method indicators except mean error, and the Average method except std error. Greater N A improves all indicators except std error for Average method.

-time-frame ∆t: seems to bring little improvement to both methods with this data set.

-number of clusters k: for this data set, seems to have little effect overall as well.

The drawbacks of each of these parameters are: -finer speed/acceleration discretizations come at the cost of increased data occupancy (before and after clustering) and higher CPU time for the clustering.

-smaller time-frames increase data occupancy before clustering, clustering CPU time, and the risk of not having enough data on the segments for the distributions to be meaningful.

-a larger number of clusters leads to a higher data occupancy after clustering, and clustering CPU time.

Overall, the validation indicates that the Memoryless method has an acceptable mean error and a good std error, while the Average method has a very good mean error but a bad std error. The JS divergence is harder to interpret but is lower for the Memoryless method (with half of the distribution below 0.025) than the Average method (half of the distribution below 0.125).

At the core, we are interested in making use of the statistical information from the traffic distributions. In this respect, the Average method appears too limited, with a bad std error regardless of the parameters choice. Therefore, in the following sections we will focus on the Memoryless method.

For the 1-day simulation on the whole city of Luxembourg, the raw data corresponding to the 18322 segments take up to 5.7GB. Table I indicates the data size before ('HIST') and after ('BARY') clustering, as well as the CPU time for the clustering, for a 10% subset of the 18322 segments. For the whole dataset, with a discretization N V = N A = 10 and a time-frame ∆t = 10min the histograms for the (speed,acceleration) distributions are computed in less than 15min and amount to 1.2GB. Setting k = 4, the clustering step takes 19h and the distributions for the barycenters occupy 59MB. Taking the histograms already reduces data occupancy by 79%, and clustering pushes it to 99% in total. More generally, the data reduction done by the histograms is related to the speed/acceleration discretization, while the clustering further reduces data size according to the ratio k/∆t. Memory usage during the computations was about 100MB and posed no difficulties. It should be noted that the clustering is an offline step that only has to be done once. Also, it is independent for all segments and may benefit heavily from parallelization.

In the end, we observe that the range of parameters that give a reasonable trade-off between accuracy, data size and computation times is rather small, which limits the interest of trying sophisticated techniques to deduce some 'optimal' settings. The number of clusters k is discussed again in V. As for the N V , N A discretization, possible improvements could be for instance the kind of adaptive exploration described in [START_REF] Yang | Pace: a path-centric paradigm for stochastic path finding[END_REF], provided the overall cost of the exploration does not outweight the cost of simply taking a sufficient fixed discretization.

V. RESULTS USING REAL DATA

Thanks to floating car data collected by the smartphone application Geco air [START_REF]Geco air[END_REF], we were able to test our method on a real-life scenario. We focused on a portion of the A7 highway near Lyon, France, which is known to be regularly used by commuters. For our analysis, the traffic data collected during the working days of the last two years were aggregated as they were recorded over one day, reasonably assuming that the data share similar traffic patterns. The speed measurements were then divided into 10-minutes time-frames. The Memoryless method is applied to the 1632 collected speed profiles. The discretization of the (speed,acceleration) space is limited at a 10 × 10 grid to reduce the computation time. One of the main goal of our method is to represent traffic with a small number of (speed,acceleration) distributions, thus we choose to explore the traffic representation with a number of clusters between 1 and 10. As we did with the data from SUMO, we will compare the consumption distribution of the Memoryless method C M SM to the reference consumptions C Ref .

A. Clustering impact

The choice of the number of clusters k is related to a tradeoff between data size, computational time and accuracy of the consumption estimate. A small number of clusters will reduce the data occupancy and CPU time, with the risk of a coarser consumption estimate. A larger k will give consumptions closer to the unclustered raw data, at the cost of increased data size and CPU time. Fig. 8 shows the mean and standard deviation errors between the consumption distributions C Ref and C M SM , for a number of clusters k ∈ 1, • • • , 10, as well as the unclustered case (dashed line). We see that the mean error is below 5%, and further decreases below 2% for k = 2, 3, 4 clusters, while larger k have larger errors. Concerning the standard deviation, it is below 0.4 and appears globally decreasing with k. The few quirks in the curve may be due to some lingering sensitivity of the barycenter computation to the initialization. In the following, we set k = 4 clusters since this value gives a good trade-off. The discussion in V-C shows the link between the clusters and the traffic conditions, and sheds some light on why a small number of clusters seems optimal, since typical traffic conditions likely encompass only a limited number of situations (fluid, jammed, etc). Errors for the unclustered data and the mean speed method are indicated as well. The value k = 4 seems to be a good choice in terms of errors and computational cost.

B. Mean speed method

Furthermore, we decide to compare our statistical model based on (speed,acceleration) distributions, to a simpler approach using only mean speeds, available from tools such as HERE or Google Maps. Noting vs,t the mean speed for a given pair (segment,timeframe), we compute a basic estimate of the energy consumption C AV GS with a constant speed and null acceleration

C AV GS = P (v s,t , 0)L vs,t . (10)
Taking our set of (speed,acceleration) profiles from section V, we recompute the mean speed information. We show on Fig. 9 the probability distributions for the reference consumption, memoryless method and mean speed method. We observe that the memoryless results are closer to the reference, with relative errors (ε mean = 4.6%, ε σ = 24.4%) better than (ε mean = 22.7%, ε σ = 43.6%) for the mean speed method. ENERGY CONSUMPTION (J)

C. Analysis of the 4-cluster case

In the following we focus on the 4-cluster case that seems to offer a good trade-off between model size and accuracy. On Fig. 10 we plot the subset of speed profiles belonging to each of the 4 clusters. We observe that these speed profiles do appear rather similar in each cluster. Clusters 1 to 3 correspond to relatively smooth traffic conditions, with little speed variations, 2 being the fastest, followed by 1 then 3. Cluster 4, on the other hand, obviously corresponds to a traffic jam situation, with large variations in speed and frequent drops to null speeds. Fig. 11 represents the same 4 subsets of profiles for each cluster in the (speed,acceleration) space. This representation confirms that the (speed,acceleration) distributions are quite distinct for each cluster. The level sets on each graph correspond to the distribution of the barycenter of the cluster, and we observe that the barycenters coincide rather well with the profile subsets. Furthermore, we can interpret the clusters in terms of traffic conditions, as summarized on Tab. II. In order to see whether the traffic interpretation is coherent with reality, the clusters associated with each time-frame were plotted in Fig. 12, in which the grey portions correspond to the time-frames with not enough data for the analysis. We can observe that the cluster 3, corresponding to a dense traffic, appears in the morning between 7am and 9:30am. The cluster 4, traffic jam, is present essentially between 4:30pm and 8pm. Cluster 1, normal traffic conditions, during the day and cluster 2, fluid conditions, during the night. Overall, this traffic pattern seems very consistent with reality, exhibiting two peak hours in the morning and evening, typical of commuting behavior.

D. Summary

The analysis conducted on real data confirms that representing traffic by means of (speed,acceleration) distributions is effective not only to estimate the energy consumption distributions in different portions of the road network but also to identify different traffic conditions. Also, clustering proves to be an effective method to reduce data occupancy.

The cluster stage fits its role, allowing us to keep a reasonable data size while retaining most of the useful information from the original set of (v, a) profiles.

VI. CONCLUSIONS

In this paper, we have presented a new approach to use traffic data to predict the energy consumption of vehicles. The key point is to consider the (speed,acceleration) data in a statistical sense without the temporal aspect, coupled with a decomposition of the road network into a collection of small segments, based on topological aspects.

Numerical experiments carried out with traffic data generated by the traffic simulator SUMO indicate that our approach is able to reconstruct the distribution of the energy consumption over a set of vehicles. We introduce two methods to compute the energy consumption, called Average and Memoryless Sampling methods. The Memoryless Sampling method gives a more accurate estimate of the distribution of energy consumptions, according to indicators such as std error and Jensen-Shannon divergence.

We also investigate the influence of several parameters such as the (speed,acceleration) discretization, length of time-frame for data aggregation, and number of clusters for the data reduction.

The analysis on real data shows that the Memoryless Sampling method performs better than a more basic approach based on mean speed only. Another interesting point is that the clusters are consistent with the traffic conditions.

Possible future works include a second level of clustering, creating clusters of road segments with close traffic conditions, and variants of the Memoryless Sampling method that would retain some temporality (e.g. Markov). Another direction could be to investigate the possible coupling of the presented model with approaches in flow management such as the SS-CTM in [START_REF] Jin | Throughput-improving control of highways facing stochastic perturbations[END_REF], since both use a spatial decomposition and account for stochastic aspects of the traffic.

5 Fig. 1 :

 51 Fig. 1: Example of (speed,acceleration) distribution.

Fig. 2 :Fig. 3 :

 23 Fig. 2: Jensen-Shannon divergence for a noised Gaussian distribution with base parameters µ = σ = 1.

Fig. 4 :

 4 Fig. 4: Cumulated consumption distribution for one segment -effect of (v, a) discretization and clustering. Reference consumption (blue line) is recomputed along each vehicle speed/acceleration profile. The two graphs show the consumptions estimated by the memoryless sample and the average consumption methods. We compare several settings for the speed/acceleration discretization, as well as the number of clusters.

Fig. 5 :

 5 Fig. 5: Cumulated consumption distribution for one segment -effect of the time-frame duration.

  Cumulated distribution of relative Standard deviation error. Cumulated distribution of Jensen-Shannon divergence.

Fig. 6 :

 6 Fig. 6: Speed discretization.

  (a) Cumulated distribution of relative Mean error.

  Cumulated distribution of Jensen-Shannon divergence.

Fig. 7 :

 7 Fig. 7: Acceleration discretization.

Fig. 8 :

 8 Fig.8: Memoryless sampling method. Relative mean error and standard deviation for different numbers of clusters k. Errors for the unclustered data and the mean speed method are indicated as well. The value k = 4 seems to be a good choice in terms of errors and computational cost.

Fig. 9 :

 9 Fig.9: Probability distributions of energy consumption. Reference (blue) vs Memoryless sampling method (red) and mean speed method (yellow). The mean and std are indicated for the 3 discrete distributions.

4 Fig. 10 :

 410 Fig. 10: Speed profiles associated with the clusters.

4 Fig. 11 :

 411 Fig. 11: Real distributions compared to their cluster.

Fig. 12 :

 12 Fig. 12: Traffic clusters according to the time of day.

Frederic

  Bonnans Got an engineering degree from the Ecole Centrale Paris, a PhD (docteur ingenieur)from the Université Technique de Compiègne, and the Habilitation from Université Paris IX Dauphine in 1979, 1982 and 1993 resp. He is currently a senior researcher at Inria Saclay and Ecole polytechnique. His research theme is dynamic optimization for deterministic and stochastic systems, and applications to energy management and biological models. Giovanni De Nunzio received the B.Sc. and M.Sc. degrees in information and automation engineering from the University of L'Aquila, Italy, in 2007 and 2010. He received the PhD in automatic control from the Grenoble Institute of Technology, France, in 2015. He is currently a research engineer at IFP Energies nouvelles, Lyon, France. His research activities focus on modeling, control, and simulation of traffic and transportation systems, vehicle powertrain modeling, advanced driver-assistance systems, energy consumption minimization, and graph theory. Thomas Leroy graduated from ESSTIN and received the M.Sc. degree in control theory from the University Henri Poincaré in Nancy, France, in 2006. He received the doctorate degree in Control theory and Mathematics from the Ecole des Mines ParisTech, France, in 2010. He is currently a research engineer at IFP Energies nouvelles, Rueil-Malamaison, France. His research activities focus on modeling, control, and simulation of internal combustion engines and hybrid electric powertrains. Pierre Martinon received his M.Sc. and Ph.D. degrees in applied mathematics from the ENSEEIHT (Toulouse, France) in 2001 and 2005. He is currently a researcher at Inria Saclay and Ecole polytechnique. His research field is optimization for dynamical systems, more specifically the numerical aspects.

  . . . , k}, μj <> μj do Backup all barycenters μj ← μj , ∀j ∈ {1, . . . , k} for each distribution µ i do Find the closest barycenter μc Label µ i as part of cluster c Update all cluster barycenters μj , ∀{j ∈ 1 . . . k}

	• • • µ I
	Output: Clusters barycenters μ1 • • • μk
	Initialization
	for j ∈ {1, . . . , k} do
	μj ← choose randomly in µ 1 • • • µ I
	Iterate until convergence
	Initialize barycenter backups μj ← 0
	while ∃j ∈ {1,

TABLE I :

 I Influence of parameters on CPU time and data size (HISTograms before clustering and BARYcenters after). Dataset is 10% of the total 18322 segments of the Luxembourg scenario, with raw size about 480MB.

TABLE II :

 II Traffic Interpretation

	Cluster	Mean Speed Speed Spread	Acceleration Interpretation
	1	Medium	Important	Small	Normal
	2	High	Low	Small	Fluid
	3	Low	Low	Small	Dense
	4	Very low	High	Well spread	Traffic jam

For the sake of simplicity, we assume the same number of time-frames for all the segments

which happens with probability 1 since vehicles never stop indefinitely.

in practice integration is done by the Euler scheme.

We set N T the number of time-frames for each segment 1 . We denote the family of time-frames (t i,s ) (1...Nt)×(1...N S ) , and their length (∆T i,s ) (1...Nt)×(1...N S ) . For each segment s, we record the (speed,acceleration) of all the vehicles passing through the segment during all the time-frames t i,s .

In the following, we work with discrete distributions in the (v, a) space. We denote V and A the sets of feasible speed and acceleration. To simplify, these sets are taken identical for all segments and time-frames, thus all (speed,acceleration) discrete probability distributions have the same support in V × A. We denote N V and N A the discretization size of V and A. Recalling that N S and N T are the number of road