Sergio Pirozzoli

Matteo Bernardini

Simon Marié

Francesco Grasso

Early evolution of the compressible mixing layer issued from two turbulent streams

Keywords: compressible turbulence, shear layer turbulence, turbulence simulation

Direct numerical simulation of the spatially developing mixing layer issuing from two turbulent streams past a splitter plate is carried out under mild compressibility conditions. The study mainly focuses on the early evolution of the mixing region, where transition occurs from a wake-like to a canonical mixing-layer-like behaviour, corresponding to the filling-up of the initial momentum deficit. The mixing layer is found to initially grow faster than linearly, and then at a sub-linear rate further downstream. The Reynolds stress components are in close agreement with reference experiments and follow a continued slow decay till the end of the computational domain. These observations are suggestive of the occurrence of incomplete similarity in the developing turbulent mixing layer. Coherent eddies are found to form in the close proximity of the splitter plate trailing edge, that are mainly organized in bands, initially skewed and then parallel to the spanwise direction. Dynamic mode decomposition is used to educe the dynamically relevant features, and it is found to be capable of singling out the coherent eddies responsible for mixing layer development.

Introduction

Turbulent mixing phenomena are ubiquitous in nature, from small-scale mixing in the bloodstream to the large scales of oceanic and atmospheric streams. This research subject has been extensively studied in the past and is also of major interest for industrial engineering systems, in which strongly heterogeneous flows can mix in a wide variety of situations. Turbulence is known to be a major mixing catalyst, and its dynamical behaviour is a crucial issue for the understanding of mixing processes [START_REF] Dimotakis | Turbulent mixing[END_REF]. Flow momentum redistribution is also responsible for the radiation of acoustic waves, and the study of mixing layers is a key building block for the understanding of jet-radiated noise [START_REF] Tam | Supersonic jet noise[END_REF].

Turbulent mixing is typically realized by merging of streams having different properties in terms of speed, temperature, and chemical composition, past a separating surface (in experiments, a splitter plate). Since early studies [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF], it has been known that large coherent structures form past the trailing edge of the splitter plate, which are mainly organized as rollers with significant spanwise coherence, and being the main cause of the entrainment of fresh fluid into the mixing zone. Later, it became clear that the spatial organization of the coherent structures in the mixing layer is significantly affected by compressibility effects [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF], by the state of the upstream boundary layers [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF][START_REF] Slessor | Turbulent shear-layer mixing at high Reynolds numbers: effects of inflow conditions[END_REF], and by the shape of the trailing edge of the splitter plate [START_REF] Laizet | Direct numerical simulation of a mixing layer downstream a thick splitter plate[END_REF]. Regarding the first issue, [START_REF] Papamoschou | The compressible turbulent shear layer: an experimental study[END_REF] observed that the growth rate of the mixing layer is significantly reduced in the presence of mild compressibility effects, and the controlling parameter is the convective Mach number, defined as

Ma c = u 1 -u 2 c 1 + c 2 , (1.1)
where u and c are the undisturbed velocity and sound speed associated with the high-(subscript 1) and the low-speed (subscript 2) stream, respectively. Regarding the second issue, it should be noted that the state of the two merging streams may be laminar, turbulent, or transitional. In this respect, [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF] noticed significant history effects, with an influence of the upstream state even on the asymptotic structure of the developing mixing layer. Specifically, it was found that tripped upstream boundary layers yield lower growth rates than untripped boundary layers, although the mixing layer structure at about 150 boundary layer thicknesses downstream of the splitter plate is similar to within 10 % for the Reynolds stress components. In most analyses presented in the literature both boundary layers are laminar, or one is laminar and one is turbulent [START_REF] Sandham | Direct numerical simulation of the early development of a turbulent mixing layer downstream of a splitter plate[END_REF]. A notable exception is the study of [START_REF] Goebel | Mean and turbulent velocity measurements of supersonic mixing layers[END_REF], [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF], who considered a wide range of flow conditions, in terms of speed and density ratios of the two streams, with two turbulent incoming boundary layers. Although the (alleged) asymptotic developed state of turbulent mixing layers has been extensively studied [START_REF] Bradshaw | The effect of initial conditions on the development of a free shear layer[END_REF], much less information is available from experiments regarding the non-similar region just past the trailing edge of the splitter plate, mainly because of technical difficulties in obtaining the necessary resolution in experimental facilities.

In one of the few dedicated studies, [START_REF] Ghaemi | Counter-hairpin vortices in the turbulent wake of a sharp trailing edge[END_REF], used tomographic particle image velocimetry to investigate the evolution of the wake forming past an airfoil trailing edge (hence, for the special case of unit velocity ratio), under low-speed conditions. Those authors highlighted the formation of counter-hairpin vortices, and concluded that these structures were formed as a result of the generation of spanwise vortices along the wake centreline, upon the disappearance of the solid wall. Direct numerical simulation (DNS) of turbulent mixing layers has been possible only in relatively recent times. To keep the computational cost within acceptable bounds, most numerical studies have been carried out in a temporal setting, by monitoring the time evolution of spatially periodic mixing layers. Temporal simulations allow some insight into the early evolution of mixing layers, and in fact in the pioneering DNS by [START_REF] Rogers | Direct simulation of a self-similar turbulent mixing layer[END_REF] the two streams were initialized with the same velocity field obtained from a DNS of a turbulent boundary layer, to realistically reproduce the flow conditions downstream of a splitter plate. Temporal simulations have also shown that it is possible that the ultimate self-similar state may be dependent on details of the initial conditions [START_REF] George | Some new ideas for similarity of turbulent shear flows[END_REF], and provided insight into the growth reduction mechanism associated with compressibility effects [START_REF] Vreman | Compressible mixing layer growth rate and turbulence characteristics[END_REF][START_REF] Foysi | The compressible mixing layer: an LES study[END_REF]. Spatial mixing layer simulations have been rather limited, so far, and they have been typically carried out by assuming a laminar mixing layer profile at the inflow with superposed small perturbations, consisting of a combination of the most unstable modes to accelerate transition to a turbulent state [START_REF] Fu | Numerical simulation of compressible mixing layers[END_REF][START_REF] Wang | Coherent fine scale eddies in turbulence transition of spatially-developing mixing layer[END_REF][START_REF] Zhou | Direct numerical simulation of a spatially developing compressible plane mixing layer: flow structures and mean flow properties[END_REF]. This approach allows the capture of some particular features of the spatial shear layer, including the non-symmetric spreading, but it suffers from the same limitations as temporal studies, in terms of unrealistic initial development. Notable exceptions include the study of [START_REF] Laizet | Direct numerical simulation of a mixing layer downstream a thick splitter plate[END_REF], who analysed the effect of trailing edge shape on the process of mixing layer formation and that of [START_REF] Sandham | Direct numerical simulation of the early development of a turbulent mixing layer downstream of a splitter plate[END_REF], who investigated the mixing of a laminar and a turbulent stream in the low-Ma c regime by large eddy simulation (LES), including a splitter plate with zero thickness in the computation. These studies showed significant differences with respect to the case of mixing layers initiated through synthetic disturbances, and highlighted the significant influence of the initial velocity deficit on the fully developed mixing layer region. [START_REF] Sharma | Large-eddy simulation of supersonic, turbulent mixing layers downstream of a splitter plate[END_REF] further observed that the merging of a laminar and a turbulent boundary layer results in faster development of self-similar velocity profiles as compared to the case of two incoming laminar boundary layers. They also observed that in the former case the flow adjustment from the wake-dominated regime to the mixing layer regime occurs without any major instability, resulting in a significantly weaker radiated acoustic near field.

The goal of the present paper is to study the process of mixing layer formation in a flow configuration with sizeable compressibility effects, by numerically reproducing the flow conditions of a well-documented flow case with two turbulent streams, namely flow case 2 of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF]. For that purpose we perform a large-scale DNS in a wide computational domain which includes the two upstream turbulent boundary layers developing on the two sides of a zero-thickness splitter plate, and their early merging region. To our knowledge, this type of computation has not been attempted before. The paper is organized as follows. The numerical setup for the DNS is described in § 2; the main velocity statistics (including those related to the early mixing layer formation and growth) are reported in § 3, where the quality of the DNS is assessed through detailed comparison with experimental data; the structural organization of the flow field, including a study of the scale change resulting from mixing layer formation, is analysed in § 4; flow snapshots are processed using dynamic mode decomposition (DMD) in § 5, to bring out dynamically relevant flow features; and concluding remarks are given in § 6.

Computational setup

The computational arrangement for the DNS is sketched in figure 1. Two turbulent boundary layers having the same thickness (δ 0 ) at the inflow of the computational domain, develop on the two sides of a splitter plate with zero thickness and merge past the trailing edge, forming a turbulent mixing layer. The mesh is uniform along the streamwise (x) and spanwise (z) direction, and stretched along the wall-normal direction to have y + 1 at the wall. The computational box has a size of L x × L y × L z = 120δ 0 × 47δ 0 × 12δ 0 , to be able to have fully developed upstream boundary layers, and follow the mixing layer evolution into the developed state. The computational domain has been selected based on previous experience [START_REF] Foysi | The compressible mixing layer: an LES study[END_REF][START_REF] Sharma | Large-eddy simulation of supersonic, turbulent mixing layers downstream of a splitter plate[END_REF], to minimize numerical confinement effects, within the constraints of computational cost. The origin of the coordinate system is set at the trailing edge of the splitter plate, at x 0 ≈ 40δ 0 downstream of the inflow. The Stream 1, where subscripts 1 and 2 are used to denote the upper and lower streams, respectively. While the free-stream conditions are made to fit those of case 2 of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF], the properties of the two boundary layers are somewhat different, mainly as we had to reduce the Reynolds numbers by about a factor of three (in the experiment of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF], Re δ 1 ≈ 51 765, Re δ 2 ≈ 56 286) to make the simulation feasible. Furthermore, in the experiment the two streams were brought together to form the mixing layer at an angle of 2.5 • across a thin splitter plate with a base height of approximately 0.5 mm (18 % of the averaged boundary layer thickness). As shown in the discussion below, these differences seem to have a very minor impact on the mixing layer development. The computational mesh includes N x × N y × N z = 6144 × 301 × 768 points, and the spacing in the wall-parallel directions (taken just upstream of the trailing edge) is x + = 4.42, z + = 3.52 for the upper stream, and x + = 9.83, z + = 7.68 for the lower stream, where the superscript + indicates normalization by local wall units, based on the friction velocity u τ and the viscous legthscale δ v = ν w /u τ , being ν w the kinematic viscosity evaluated at the wall. In our experience [START_REF] Pirozzoli | Probing high-Reynolds-number effects in numerical boundary layers[END_REF], and based on the results of the next section, the mesh resolution is sufficient to have well-behaved boundary layer statistics.

Ma ∞ T ∞ (K) T w (K) Re δ Re θ Re τ r = U 2 /U 1 s = ρ 2 /ρ 1 δ 1 /δ 2 θ 1 /θ 2 Ma c Π c 1
The numerical simulation has been carried out by means of a conservative finite-difference algorithm. The flow solver relies on central sixth-order discretization of the convective terms of the Navier-Stokes equations cast in fully split form [START_REF] Pirozzoli | Numerical methods for high-speed flows[END_REF], which guarantees discrete conservation of the total kinetic energy in the limit case of inviscid, incompressible flow, also in the presence of grid stretching in the coordinate directions. This approach allows robust spatial discretization of the convective terms without the addition of any spurious numerical dissipation in the form of upwinding or filtering. Being cast in locally conservative form, the method also guarantees excellent computational efficiency, facilitating hybridization with shock-capturing methods. A fifth-order weighted-essentially-non-oscillatory (WENO) scheme is used to replace central differencing in the presence of discontinuities, with switch controlled by the Ducros sensor [START_REF] Bernardini | The wall pressure signature of transonic shock/boundary layer interaction[END_REF]. It is worth noting that in the present simulation WENO is only activated at about 0.36 % of the total number of mesh points, hence artificial viscosity has very little impact on the computed solution. The diffusive terms in the Navier-Stokes equations are expanded to Laplacian form for improved numerical stability, and also discretized with sixth-order central-difference formulae, to guarantee proper action of molecular diffusion at the smallest scales resolved on the computational mesh. The resulting semi-discrete system of equations is advanced in time by means of an optimized, fully explicit fourth-stage Runge-Kutta algorithm [START_REF] Bernardini | A general strategy for the optimization of Runge-Kutta schemes for wave propagation phenomena[END_REF]. Particular attention is devoted to the correct prescription of the inflow boundary conditions, which is a key ingredient in the simulation of spatially developing turbulent flows. In the present simulation the inflow conditions are prescribed through a recycling-rescaling procedure, suitably adapted to compressible flows [START_REF] Bernardini | The wall pressure signature of transonic shock/boundary layer interaction[END_REF], and the recycling station is placed sufficiently far from the inlet (about 30δ 0 ) to achieve full decorrelation from the inflow, as seen from inspection of the streamwise and spanwise two-point velocity correlations. In the absence of detailed information from experiments, the two walls are assumed to be adiabatic.

In the following, unless otherwise specified, the velocity difference between the two streams, U = U 1 -U 2 , is used to normalize the velocity statistics, and lengths are made non-dimensional with respect to the arithmetic average of the 99 % thicknesses of the upper and lower boundary layers (δ = (δ 1 + δ 2 )/2) evaluated just upstream of the trailing edge. The following scaled coordinates will then be used:

ξ = (x -x 0 )/δ, η = y/δ, ζ = z/δ.

Velocity statistics

The statistical properties of the incoming boundary layers are provided in figure 2, where we show the mean velocity profile in van Driest scaling, u VD = u 0 (ρ/ρ w ) 1/2 du, and the density-scaled Reynolds stress components, R ij = (ρ/ρ w )u i u j (here the overbar denotes Reynolds-averaged values, and the primes fluctuations thereof). In figure 2(a) the linear viscous sublayer law and the logarithmic overlap layer law are also shown for reference, with standard values of the log-law constants, k = 0.41, C = 5. The wake region is found to be stronger for the lower stream, which has higher friction Reynolds number, and consistently the streamwise velocity variance is also found to be somewhat higher than for the upper stream. Overall, the figure suggests that the incoming boundary layers are fully developed, thus ensuring proper initiation of the mixing process of the two streams.

Mean vorticity contours are shown in a side view in figure 3. A weak compression wave is observed to form at the trailing edge on the high-speed (upper) side, whereas a weak expansion wave forms on the low-speed (lower) side, with subsequent upward deflection of the mixing layer. It is known in the incompressible case [START_REF] Pope | Turbulent Flows[END_REF] that the deflection of the mixing layer axis occurs towards the side with less entrainment, which is the low-speed side, under ordinary circumstances. The slight upward deflection observed here is thus likely to be an effect of the non-unity density ratio. To visually appreciate the mixing layer evolution, in figure 3 we show the mixing layer mid-line (y 0 , marked as a thick white line), and the position of the lower and upper boundary of the mixing layer (y L , y U , thick black lines). These are defined by vertically scanning the mean velocity profiles as suggested by [START_REF] Goebel | Mean and turbulent velocity measurements of supersonic mixing layers[END_REF], as follows. For those profiles for which a location where u < U 2 -0.1 U does not exist (i.e. the momentum deficit is low), the upper and lower boundaries of the mixing layer are defined as the transverse locations where the mean streamwise velocity equals U 1 -0.1 U and U 2 + 0.1 U, respectively, and the mid-line is defined as the transverse position where the mean velocity equals (U 1 + U 2 )/2. For those profiles with strong momentum deficit, the lower edge of the mixing layer is defined to be the transverse location where the mean velocity equals U 2 -0.1 U, and the mixing layer centreline is defined as the average of y U and y L . This definition explains the jump in y L occurring at ξ ≈ 20.

A comparison of the evolution of the mixing layer apparent thickness, b = y U -y L , with the experimental data of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF] is presented in figure 4(a). Very good agreement is recovered, even in the early development region. It is noteworthy that the position with a sudden drop of b, which according to the previously given definition corresponds to a change from a strong velocity deficit state to a 'canonical' mixing layer state, is well captured. In figure 4 we also show the evolution of the vorticity thickness, defined as ∆ ω = U/(max y ∂u/∂y), which is the most frequently used indicator for mixing layer thickness in numerical simulations. As expected from its definition, ∆ ω shows fast growth past the splitter plate trailing edge, associated with a sudden reduction of shear, and it attains a nearly linear behaviour at a distance of about 20δ past the trailing edge. In this region, which is typically associated with a self-similar state of the mixing layer [START_REF] Pope | Turbulent Flows[END_REF], we find ∆ ω ≈ b, in agreement with the claims of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF].

Clearer scrutiny of the mixing layer thickness evolution is provided in figure 4(b), where the vorticity data are shown in a logarithmic scale, to highlight power-law behaviour. It is found that in the very early stages of mixing layer formation (ξ 3), the vorticity thickness increases faster than linearly, as

∆ ω ∼ ξ 1.36 . (3.1)
This behaviour can be tentatively explained as follows. Referring to figure 5, we assume that: (i) the merging boundary layers have the same thickness δ (see table 1), and have a power-law mean velocity profile, hence u i (y) = U i (y/δ) α for y δ; (ii) the mixing layer spreads symmetrically with respect to the geometric mid-plane, and ∆(ξ ) is its local thickness; (iii) for |y| ∆/2, the boundary layers are not affected by the entrainment process. From these hypotheses, taking as mixing length the local mixing layer thickness, and the local velocity difference across the mixing layer, u 1 (∆/2) -u 2 (∆/2) as the typical eddy velocity scale, it follows that the relevant eddy viscosity at a given streamwise station is

ν t ∼ ∆ • (u 1 (∆/2) -u 2 (∆/2)) = ∆ • U • (∆/(2δ)) α . (3.2)
Dimensional arguments easily yield the mixing layer thickness evolution equation

d∆ dξ ∼ 1 U 0 ν t ∆ = U U 0 ∆ 2δ α , (3.3)
where U 0 = (U 1 + U 2 )/2 is an average convection velocity, whence it readily follows

∆(ξ ) ∼ ξ 1/(1-α) . (3.4)
Incidentally, after the mixing layer has grown to a thickness greater that the boundary layer thickness, the above analysis returns the classical linear growth law [START_REF] Pope | Turbulent Flows[END_REF]. However, until saturation is complete, (3.4) suggests that the growth rate should be faster than linear. At the low Reynolds number of the boundary layers under consideration, power-law exponents are expected to be in the range 0.182 α 0.196 [START_REF] Pope | Turbulent Flows[END_REF], which we have verified by fitting the DNS velocity profiles. The growth rate exponent predicted by (3.4) is then in the range between 1.22 and 1.24, not too far from the DNS data fit of (3.1).

It is also noteworthy that, although figure 4(a) seems to suggest that linear growth is achieved past ξ ≈ 20, figure 4(b) in fact favours power-law growth, with ∆ ω ∼ ξ 0.587 .

(3.5)

In retrospect, a similar behaviour can be traced in the numerical simulation of [START_REF] Sandham | Direct numerical simulation of the early development of a turbulent mixing layer downstream of a splitter plate[END_REF], for the merging of a laminar and a turbulent stream. This non-canonical behaviour will be further commented on below. A comparison of the velocity statistics with the experimental data of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF] is presented in figure 6 at various stations across the developing mixing layer, where the cross-stream coordinate is scaled with respect to b, and the reference velocity is U. The scaled mean velocity profiles (left-most column) highlight the presence of significant momentum defect in the proximity of the plate trailing edge caused by flow retardation within the issuing boundary layers, which is gradually filled up as the flow evolves, although remnants are still visible (see also the analysis below) near the outflow of the computational domain. Overall, agreement with the experiments is quite remarkable, also given the previously noticed discrepancy in the Reynolds numbers and the assumption of a splitter plate with zero thickness. Similar considerations also apply to the streamwise turbulence intensity (second column from left), whose main peak is well predicted throughout, occurring in all cases near the mixing layer mid-line, where turbulence kinetic energy production is maximum. A shoulder in the streamwise velocity variance distribution is observed in the lower part of the mixing layer, approximately corresponding to the inflectional point of the momentum defect zone. Poorer agreement is found for the vertical Reynolds stress component (third column), which is overestimated with respect to experiments by up to 50 %, and which is likely to be due to difficulties in measuring the cross-stream velocity with laser Doppler velocimetry. This also implies, although to a lesser extent, discrepancies in the prediction of the shear stress peak (last column). The turbulence levels away from the mixing layer core are found to be substantially higher in experiments, presumably because of sound background wind tunnel turbulence. Similar results for the fully developed state are obtained from spatial simulations with synthetic inflow conditions [START_REF] Salvadore | GPU accelerated flow solver for direct numerical simulation of turbulent flows[END_REF].

0.5 0 -0.5 -1.0 -1.5 -2.0 1.0 1.5 2.0 0.5 0 -0.5 -1.0 -1.5 -2.0 1.0 1.5 2.0 0.5 0 -0.5 -1.0 -1.5 -2.0 1.0 1.5 2.0 0.5 0 -0.5 -1.0 -1.5 -2.0 1.0 1.5 2.0 0.5 0 -0.5 -1.0 -1.5 -2.0 1.0 1.5 2.0 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 -1 0 1 2 -2 0 0.
The evolution of the maximum velocity deficit across the mixing layer, defined as U 2 -min y u(y), is shown in figure 7. Based on its definition, this property is expected to tend asymptotically to zero as x → ∞, where a canonical mixing layer behaviour should be recovered. Figure 7 shows a two-stage evolution whereby the momentum defect initially scales as ξ -1/2 , and then decays exponentially past ξ ≈ 20. This behaviour can be explained by the prevalence of a wake-like behaviour in the trailing edge near field, recalling the scaling of momentum defect in a planar wake [START_REF] Pope | Turbulent Flows[END_REF], which is then followed by rapid fill-up in the far field. It is noteworthy that linear stability analysis of mixing layers in the presence of momentum defect [START_REF] Zhuang | Instability of wake-dominated compressible mixing layers[END_REF] suggested the existence of two modes of instability, namely a wake mode which is dominant in the presence of strong deficit, superposed onto a canonical mixing layer mode. The results of the stability analysis appear to be consistent with the observed two-stage evolution of the mixing layer properties.

The evolution of the Reynolds stress peaks with the distance from the trailing edge is depicted in figure 8. A sudden rise of all Reynolds stress components is observed past the trailing edge up to ξ ≈ 3, which is associated with flow rearrangement, and then followed by monotonic decay. This behaviour is overall consistent with the findings of [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF] for a low-speed tripped mixing layer. However, those authors were unable to measure close enough to the trailing edge to clearly probe the near-field peaks. Most previous studies, including [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF], [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF], as well as the theory of self-similar mixing layers [START_REF] Pope | Turbulent Flows[END_REF] suggest saturation of the Reynolds stresses with the trailing edge distance, when scaled by the mixing layer velocity difference. While tendency to saturation here may be inferred from a linear representation as in figure 8(a), the same data cast in a log-log representation (figure 8b) do in fact suggest continuing decrease of the Reynolds stresses in the far-field region, according to a weak inverse power law of the streamwise coordinate, and no clear tendency to saturation. In particular, the inferred behaviour for the turbulent shear stress, which is the only component entering the mean streamwise momentum equation, is -u v ∼ ξ -0.33 , for ξ 10. This finding is consistent with the previously noticed nonlinear growth of the mixing layer thickness (recalling figure 4b). In fact, analysis of the mean streamwise momentum equation [START_REF] Townsend | The Structure of Turbulent Shear Flow[END_REF][START_REF] Pope | Turbulent Flows[END_REF] indicates that a necessary condition for self-similarity is d∆ dξ

∼ u 2 0 , (3.6)
where u 0 is the typical velocity scale of turbulent fluctuations. The power-law exponents found in DNS for the turbulent shear stress and for the mixing layer thickness (such as from (3.5), which yields d∆ ω /dξ ∼ ξ -0.413 ) are in satisfactory agreement with (3.6). However, full self-similarity also requires similarity of the Reynolds stress balance equation, which can only be achieved if ∆ ∼ ξ , u 0 ∼ const. [START_REF] Townsend | The Structure of Turbulent Shear Flow[END_REF]. The present data then point to the attainment of a state of incomplete similarity, whereby equilibrium of the mean velocity field is reached prior to equilibration of the velocity fluctuations field. While the existence of such a state does not seem to be very surprising, to our knowledge it has not been theoretically predicted nor numerically studied before, and probably deserves some further investigation. Full similarity, hence saturation of the Reynolds stresses, is not reached in the present simulation as it probably requires a much longer fetch, and indeed the experiments of [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF] showed evidence for saturation only at ξ 150. Analysis of the Favre-averaged stresses shows differences of no more than 1 % with respect to the Reynolds averages, hence density variations do not seem to be responsible for the observed scalings. In any case, it is noteworthy (see table 2) that the mixing layer state (in terms of the Reynolds stress peaks) at the right-most station available is comparable with the asymptotic structure quoted in experiments [START_REF] Bell | Development of a two-stream mixing layer from tripped and untripped boundary layers[END_REF][START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF].

Although the mixing layer growth rate is found to be significantly different from linear, we have nevertheless attempted to compare with existing experimental data, by locally fitting the vorticity thickness distribution in figure 4 with a straight line towards the end of the computational domain. The results of this exercise are summarized in table 2. As probably to be expected, the table highlights significant differences in the alleged asymptotic growth rates from incompressible experiments, as well as from experiments with the same flow conditions. In this respect we note that the generally accepted correlation for incompressible mixing layers [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF]

d∆ dx i = C δ (1 -r)(1 + s 1/2 ) 2(1 + rs 1/2 ) , (3.7)
by taking C δ = 0.165 [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF], would yield d∆/dx| i = 0.0465. Hence, the present DNS would yield a growth rate reduction factor with respect to the incompressible case (d∆/dx)/(d∆/dx) i ≈ 0.601, whereas the data of [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF] yield (d∆/dx)/(d∆/dx) i ≈ 0.816. Several (more or less empirical) formulae have been proposed to relate compressible growth rates to incompressible ones, one of the most well established being by [START_REF] Slessor | Turbulent shear-layer mixing: growth-rate compressibility scaling[END_REF], who stipulated the correction factor to be

(d∆/dx)/(d∆/dx) i = (1 + αΠ 2 c ) -β , (3.8)
where α ≈ 4, β ≈ 0.5, and Π c is an 'effective' convective Mach number

Π c = max √ γ i -1 c i U,
(3.9) thus also covering the case that the two streams have different specific heat ratios, γ i . Applied to the flow conditions under study here, Π c = 0.52, which translates into a growth rate reduction factor (d∆/dx)/(d∆/dx) i = 0.69, which is not too far from the present DNS. It should also be noted that the growth reduction factor observed here is also entirely consistent with the findings of [START_REF] Clemens | Large-scale structure and entrainment in the supersonic mixing layer[END_REF], who reported (d∆/dx)/(d∆/dx) i = 0.63 at Ma c = 0.42, and (d∆/dx)/(d∆/dx) i = 0.56 at Ma c = 0.50.

Flow organization

Contours of velocity, temperature and pressure from instantaneous flow fields are shown in figure 9, in a side view. The figure highlights well the slow process of mixing of the two streams which leads to strongly convoluted interfaces with the surrounding irrotational fluid. However, no clear sign of flow organization into coherent spanwise vortices (rollers) is observed, as in classical visualizations [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF]. This is a likely effect of flow compressibility which acts to suppress the Kelvin-Helmholtz instability. Indeed, in the experiments of [START_REF] Clemens | Large-scale structure and entrainment in the supersonic mixing layer[END_REF], the formation of well-organized rollers at a similar convective Mach number was only observed at ξ ≈ 200, well downstream of the present computational outflow. The pressure field is found to be organized into small-sized blobs with positive/negative deviation from the free-stream value in the proximity of the trailing edge, which is a clear imprint of the small-sized vortices issued from the incoming boundary layers. Large-scale organization of the pressure field into alternating, positive and negative, delta-scaled blobs is observed further downstream, which is probably the precursor of the formation of rollers.

Top view visualizations are shown in figure 10, which indeed indicates inception of large-scale organization in the developing mixing layer. Specifically, the streamwise velocity organizes itself (as in any shear flow) into low-and high-speed streaks with spanwise spacing which monotonically grows from the trailing edge proximity (as inherited from the upstream near-wall streaks) to the computational outflow. On the other hand, the temperature and the pressure fields exhibit the signature of very wide flow structures, which eventually fill up the whole computational domain. Hints of spanwise organization are clear even at ξ ≈ 7-8, and large eddies inclined by about 70 • with respect to the free stream are observed until ξ ≈ 40. Similar diagonal patterns were also observed in the flow visualizations of [START_REF] Clemens | Large-scale structure and entrainment in the supersonic mixing layer[END_REF], although at slightly higher Mach number (Ma c 0.62). At ξ ≈ 40, the eddies become more or less parallel to the spanwise direction. It is worth recalling that temporal stability of compressible laminar mixing layers [START_REF] Sandham | Three-dimensional simulations of large eddies in the compressible mixing layer[END_REF] indicates the presence of both two-dimensional and three-dimensional unstable modes, the former expected to be dominant at Ma c 0.6, and the latter at higher convective Mach numbers. The results of the present simulation (at Ma c ≈ 0.45) apparently support a similar scenario in which two-dimensional disturbances eventually occur, although strong amplification of oblique disturbances is observed in a transitory stage. Since a box sensitivity study is not feasible, these observations should be taken with some caution, being potentially affected by the spanwise size of the computational domain (here, L z /∆ ω ≈ 3.7 at the end of the computational domain).

To quantify the growth of the size of the turbulent eddies, in figure 11 we provide the spanwise spectral densities of the streamwise velocity and pressure, taken at 24, 6.37, 14.91, 25.16, 35.42, 45.67, 55.92 (the arrow indicating increasing x). The dot-dashed line in (a) denotes the trend E(k) ∼ k -5/3 , and in (b) the trend E(k) ∼ k -7/3 . the mixing layer mid-plane. As already observed from the flow visualizations, the spectra show an overall shift to lower wavenumbers (i.e. larger scales) on progressing downstream, and the formation of an inertial layer in the u spectrum, with exponent not far from -5/3. On the other hand, the pressure spectra do not seem to exhibit the k -7/3 inertial layer which would be expected in isotropic turbulence at sufficiently high Reynolds number [START_REF] Batchelor | Pressure fluctuations in isotropic turbulence[END_REF]. This is presumably due to the action of the mean shear which activates 'rapid' source terms related to mean shear/turbulence interaction in the Poisson equation for pressure [START_REF] Blake | Mechanics of Flow-Induced Sound and Vibration[END_REF]), or to low-Reynolds-number effects, as the Taylor-scale Reynolds number here is no larger than 80, all the way to the end of the computational domain.

The spectral densities are plotted in pre-multiplied form as a function of the distance from the trailing edge in figure 12. This type of representation has the advantage of clearly highlighting the size of the energy-containing eddies, and it has the property that equal areas underneath the curves correspond to equal energy, when a semi-logarithmic representation is used. Further, in order to scale out the effect of variation of the amplitude of the disturbances with ξ , the spectra are scaled with respect to the total energy, i.e. Êu (k) = E u (k)/u 2 . The normalized velocity and pressure spectra are shown in figure 12(a,b) as a function of kδ, and confirm the scale change effect previously noticed in figure 11, in that the spectral peak shifts to larger scales with the downstream coordinate. Apparently, collapse of the pre-multiplied spectra is found if the wavenumber is scaled with respect to the local vorticity thickness, as shown in figure 12(c,d). In this representation, the spectral ridge of u is very nearly constant starting from ξ ≈ 10, whereas the peak of p is slightly shifted to higher wavenumbers. This is further confirmed from inspection of the integral length scales, as shown in figure 13, which highlight well the proportionality of the velocity length scale to the vorticity thickness (with Λ z /∆ ω ≈ 0.2), and to a lesser extent of the pressure scale (with Λ z /∆ ω ≈ 0.6, still decreasing). It is noteworthy that the observed integral length scale of the u-containing eddies is comparable to that observed in the outer part of turbulent boundary layers, in which Λ z /δ ≈ 0.2-0.3 [START_REF] Pirozzoli | On the size of the energy-containing eddies in the outer turbulent wall layer[END_REF], which probably points to a universal scaling law for the momentum streaks in all shear flows, either free or wall bounded. Finally, the enhanced spanwise coherence of the pressure field as compared to the velocity field is consistent with the flow visualizations of figure 10.

Dynamic mode decomposition

The DMD [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]) is a technique of spectral analysis which operates on a sequence of snapshots u(t n ) obtained from either experiments or, as in our case, simulations to extract information based on Koopman's analysis of dynamical systems [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF]. Introducing a linear time advancement operator A as u(t n+1 ) = Au(t n ), the DMD algorithm extracts approximate eigenvalues and eigenmodes of a reduced-order representation of A. While having some connection with proper orthogonal decomposition (POD), DMD produces modes associated with a single frequency, whereas POD modes may in general have broadband frequency content. Furthermore, for linear systems, the DMD modes coincide with the global modes, and for perfectly time-periodic data the DMD modes coincide with the Fourier modes. The essence of the DMD can be viewed as follows, with more complete and rigorous descriptions provided in [START_REF] Rowley | Spectral analysis of nonlinear flows[END_REF] and [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF]. Snapshots of data are organized as columns into a matrix V N 1 = {v 1 , v 2 , . . . , v N } (of size M × N, say), where two subsequent snapshots v j and v j+1 are separated by a constant time interval t. Following the idea underlying Krylov subspace techniques, and in particular the Arnoldi method, it is possible to approximate the linear mapping A from one snapshot to the next using an (N -1)-dimensional snapshot basis (5.1) or in matrix form

{v 2 , . . . , v N } = A{v 1 , v 2 , . . . , v N-1 } ≈ {v 1 , v 2 , . . . , v N-1 }S,
V N 2 = AV N-1 1 ≈ V N-1 1 S, (5.2)
where the reduced-order matrix S is of companion type. Its eigenvalues (Ritz values) approximate some of the eigenvalues of A and the associated eigenvectors provide the coefficients of the linear combination that is needed to express the modal structure within the snapshot basis. S is easily estimated through a least-square procedure, by performing a QR decomposition of V N-1 1 to solve (5.2), thus obtaining S = R -1 Q T V N 2 . The DMD algorithm has been initially applied to temporal snapshots of the pressure field taken in an x-y plane (at z = L z /2), in a subdomain containing the final part of the splitter plate and the initial part of the mixing layer, extending from ξ = -2 to ξ = 11.6. A total of 300 snapshots, with time separation t = 0.2δ/ U, have been processed. The Ritz values (λ) associated with the resulting reduced-order system S are shown in figure 14 in the complex plane. Essentially all eigenvalues lie on the unit circle, which is to be expected as the system is statistically stationary. The eigenvalues of S are transformed via the mapping ω = log(λ/ t), and shown in figure 15, where the real and imaginary parts represent the exponential growth or decay rate and the temporal frequency, respectively. Since the processed data are real, all eigenvalues are real or in complex-conjugate pairs, and the associated DMD spectrum is symmetric about the zero frequency. The size and colouring of the symbols reported in the figure are indicative of the coherence of the associated modes, which is useful to distinguish large-scale energetic structures from smaller-scale, less-energetic ones. Coherence is gauged by projecting each dynamic mode onto the POD basis (U) computed from the singular value decomposition of the data sequence V N-1 1 = UΣ W . This decomposition provides the matrices U and W containing the spatial and temporal structures respectively, whose energy ranking is given by the singular values Σ [START_REF] Schmid | Dynamic mode decomposition of numerical and experimental data[END_REF].

The pressure fields associated with three representative dynamic modes (labelled in figure 15) are displayed in figure 16. Besides the one stationary mode corresponding to the mean field (not shown), the most coherent mode (a) is characterized by large blobs with alternate positive and negative pressure, covering approximately three oscillatory periods over the streamwise extent of the interrogation window, and growing in the streamwise direction. Interestingly, these organized modes arise immediately past the end of the splitter plate, and they are apparently not inherited from any large-scale structure embedded in the upstream boundary layers. A similar pattern is also retrieved in mode (b), which is characterized by structures with comparable scale, but leaning in the streamwise direction. Mode (c) in figure 16 has a similar shape as mode (a); however its associated eigenvalue has lower frequency, and the associated spatial scale is correspondingly larger. Furthermore, it is characterized by a spatial extent in the η direction which is less than for mode (a). The same analysis has been repeated for snapshots taken in the geometric mid-plane in the same streamwise interval, and the shape of the leading modes is shown in figure 17. It should be noted that, since the DMD algorithm has been applied plane-wise, the modes in figure 17 do not necessarily correspond to those of figure 16. Here the most salient feature, clearly noticeable in the leading modes, is the appearance of substantial spanwise coherence of the pressure field emerging without apparent connection with the upstream turbulent boundary layers. At this stage, it appears that the leading pressure-containing eddies are strongly undulated in the spanwise direction. These eddies may be interpreted as the precursors of the formation of mixing layer rollers, which have been observed in mixing layers at low/moderate values of the convective Mach number [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF][START_REF] Clemens | Large-scale structure and entrainment in the supersonic mixing layer[END_REF]. It is noteworthy that, while their presence cannot immediately be inferred from the instantaneous visualizations reported in figures 9 and 10, they are clearly detected by the DMD algorithm, which identifies these structures as the most significant from the dynamical point of view. To analyse the downstream evolution of these large-scale structures, the DMD analysis in the mid-plane has been repeated for snapshots in the intermediate region 12 ξ 32 (figure 18) and in the far-field region 32 ξ 56 (figure 19). It becomes evident that in the intermediate region the pressure eddies gain further spanwise coherence, organizing themselves into diagonal bands filling the whole computational box. Further downstream, the bands become almost aligned along the spanwise direction, and the spanwise modulation is much reduced. Correspondingly, the modal pressure in the x-y plane exhibits a pattern typical of mixing layer rollers. The DMD modes here appear to provide a condensed form of the information contained in the instantaneous flow visualization of figure 10.

Conclusions

We have carried out a DNS study of a mixing layer developing from two turbulent boundary layers, at relatively high convective Mach number. The study has highlighted the predictive capabilities of DNS, in that the mean velocity profile and the main velocity statistics agree very well with reference experimental results [START_REF] Goebel | Experimental study of compressible turbulent mixing layers[END_REF]. The flow is found to be organized into three parts: (i) a near field (extending up to about 2-3 boundary layer thicknesses past the end of the splitter plate), which is substantially affected by the momentum deficit inherited from the parent boundary layers, and in which the Reynolds stress components attain strong peaks; (ii) an intermediate field (extending up to about 20 boundary layer thicknesses past the splitter plate), in which the momentum deficit decays according to an inverse-square-root law, as in plane turbulent wakes; and (iii) a far field in which the momentum deficit is re-absorbed, and the flow converges to the classical mixing layer structure. The analysis shows that the mixing layer growth in region (i) is quite fast, occurring according to a power law of the trailing edge distance significantly greater than unity. This observation (common in previous spatial studies) is interpreted as the consequence of incomplete entrainment of the boundary layers into the developing mixing layer. A qualitative change in the mixing layer growth takes place in regions (ii) and (iii), in which we observe slower, sub-linear growth. Correspondingly, we find decay of the Reynolds stress components according to a weakly negative power-law, without evidence of saturation as expected in fully self-similar mixing layers. This observation has been found to be the likely result of incomplete similarity, whereby mean momentum attains an equilibrium state prior to the turbulent stress distribution. The lack of full self-similarity in the present study would then be the consequence of the limited streamwise extent of the computational domain. Establishing the spatial development scale for the achievement of full self-similarity would then be an interesting topic for a follow-up study.

In terms of flow structures, we find very little hint of spatial organization in region (i), whereas spanwise coherence is observed in region (ii), where the typical eddies appear to be skewed with respect to the flow direction, and in region (iii), where they tend to be aligned along the spanwise direction. Notably, no distinct roller of the type found in experiments [START_REF] Brown | On density effects and large structure in turbulent mixing layers[END_REF] has been observed, which is due either to compressibility effects, or again to limitation in the computational domain, or both. As expected, the size of the velocity-and pressure-containing eddies is found to grow in the streamwise direction, apparently in approximately proportional fashion to the local vorticity thickness. In terms of the latter, the normalized spanwise integral scales are found to have values comparable to those of boundary layers (scaled by the boundary layer thickness), pointing to a universal behaviour of the eddy size in all shear flows. Further insight into the flow organization is provided by the DMD of the flow, which usefully isolates the process of mixing layer formation and saturation. In particular, DMD seems to provide hints on the formation of rollers, which are not (yet) clearly visible in the instantaneous flow realizations.

  FIGURE 1. (Colour online) Sketch of the computational domain; δ 0 is the boundary layer thickness of the two streams at the inflow station. The subscript 1 is used to denote the upper (high-speed) stream, and 2 for the lower (low-speed) stream.

  FIGURE 2. (Colour online) Mean velocity profile (a) and Reynolds stresses (b) on the two sides of the splitter plate, upstream of the trailing edge (ξ = -1). Solid lines indicate the upper stream, and dashed lines the lower stream. The dot-dashed lines in (a) are the viscous sublayer linear law (u + = y + ), and the logarithmic overlap layer (u + = C + 1/k log y + , with C = 5, k = 0.41).

  FIGURE 4. (Colour online) Evolution of mixing layer thicknesses past the splitter plate, scaled by the reference upstream boundary layer thickness, δ = (δ 1 + δ 2 )/2, in linear scale (a) and logarithmic scale (b). Open symbols: apparent thickness b (squares), vorticity thickness ∆ ω (circles). Solid squares: Goebel & Dutton (1991). Dashed line: ∆ ω ∼ ξ 1.36 ; solid line: ∆ ω ∼ ξ 0.587 .

FIGURE 5 .

 5 FIGURE 5. Idealized sketch of mixing layer early development; δ is the initial thickness of the boundary layers (assumed to be the same), and ∆ is the local mixing layer thickness.

  FIGURE 6. (Colour online) Comparison of mean velocity and Reynolds stress components at various streamwise stations past the splitter plate: ξ = 3.65 (a), ξ = 9.11 (b), ξ = 19.69 (c), ξ = 36.46 (d), ξ = 54.68 (e). Lines: present DNS; symbols: (Goebel & Dutton 1991).

  FIGURE 7. (Colour online) Evolution of maximum velocity defect with trailing edge distance. The dashed line denotes the trend ξ -1/2 .

  FIGURE 9. (Colour online) Instantaneous contours of (a) streamwise velocity, (b) temperature and (c) pressure in the x-y plane. 32 contour levels are shown in the respective ranges: 0.45 < u/U 1 < 1.1, 0.65 < T/T 1 < 1.2, 0.75 < p/p 1 < 1.07. The colour scale ranges from blue (lower end) to red (higher end).

  FIGURE 10. (Colour online) Instantaneous contours of (a) streamwise velocity, (b) temperature and (c) pressure in the x-z geometric symmetry plane. 32 contour levels are shown in the respective ranges: 0.45 < u/U 1 < 1.1, 0.65 < T/T 1 < 1.2, 0.75 < p/p 1 < 1.07. The colour scale ranges from blue (lower end) to red (higher end).

  FIGURE 11. (Colour online) Spanwise spectral densities of streamwise velocity (a) and pressure (b) in the mixing layer symmetry plane, at ξ = 1.24, 6.37, 14.91, 25.16, 35.42, 45.67, 55.92 (the arrow indicating increasing x). The dot-dashed line in (a) denotes the trend E(k) ∼ k -5/3 , and in (b) the trend E(k) ∼ k -7/3 .

FIGURE 13

 13 FIGURE 12. Pre-multiplied and normalized spanwise spectral densities of streamwise velocity (k z Êu , a,c) and pressure (k z Êp , b,d), as a function of the trailing edge distance. Wavenumbers are scaled with respect to δ in (a,b), and with respect to the local vorticity thickness, ∆ ω in (c,d). Contour levels are shown from 0.01 to 0.24, in intervals of 0.01, from light to dark grey shades.

FIGURE 14

 14 FIGURE 14. (Colour online) Ritz values of the reduced system detected by DMD for temporal snapshots of the pressure field in an x-y plane between ξ = -2 and ξ = 11.6.

  FIGURE 16. (Colour online) Leading dynamic modes in an x-y plane, with snapshots taken in the interval -2 ξ 11.6, corresponding to the labelled eigenvalues in figure 15. Contours of modal pressure in colour scale from blue (negative values) to red (positive values).

  FIGURE 17. (Colour online) Leading dynamic modes in the y = 0 plane, with snapshots taken in the interval -2 ξ 11.6. Contours of modal pressure in colour scale from blue (negative values) to red (positive values).

  FIGURE 18. (Colour online) Leading dynamic mode in the y = 0 plane, with snapshots taken in the interval 12 ξ 32. Contours of modal pressure in colour scale from blue (negative values) to red (positive values). The solid lines are inclined by 70 • with respect to the streamwise direction.

  FIGURE 19. (Colour online) Leading dynamic mode in the x-y plane (a) and in the y = 0 plane (b), with snapshots taken in the interval 32 ξ 56. Contours of modal pressure in colour scale from blue (negative values) to red (positive values).

TABLE 1 .

 1 Boundary layer parameters at (x -x 0 )/δ 0 = -1. Ma ∞ and T ∞ indicate the free-stream Mach number and absolute temperature, and T w the wall temperature; δ is the 99 % boundary layer thickness, and θ is the momentum thickness. Re τ = u τ δ/ν

	w is

TABLE 2 .

 2 Asymptotic mixing layer parameters from the present DNS and from experiments. In the DNS, the growth rate is estimated by linearly fitting the vorticity thickness distribution for ξ 50, and the Reynolds stresses are taken at ξ ≈ 58.

		DNS	Bell & Mehta (1990) Goebel & Dutton (1991)
	r = U 2 /U 1 s = ρ 2 /ρ 1 d∆/dx	0.57 1.55 0.028	0.6 1 0.046	0.57 1.55 0.038
	max(u 2 )/( U) 2	0.0292	0.03	0.0289
	max(v 2 )/( U) 2	0.0143	0.015	0.0098
	max(w 2 )/( U) 2	0.0195	0.022	NA
	max(-u v )/( U) 2 0.0090	0.011	0.0086

† Email address for correspondence: sergio.pirozzoli@uniroma1.it

Acknowledgements

We acknowledge that most results reported in this paper have been achieved using the PRACE Research Infrastructure resource JUQUEEN based at the Forschungszentrum Jülich (FZJ) in Jülich, Germany. This research also used resources of the Argonne Leadership Computing Facility, which is a DOE Office of Science User Facility supported under Contract DE-AC02-06CH11357. We would like to thank J. C. Dutton for making available experimental data.


R Le Goff Latimier
email: roman.legoff-latimier@ens-rennes.fr

T Kovaltchouk
email: thibaut.kovaltchouk@ens-rennes.fr

H Ben Ahmed

B Multon

Preliminary Sizing of a Collaborative System: Photovoltaic Power Plant and Electric Vehicle Fleet

Keywords: photovoltaic energy, electric vehicle, collaborative system, production forecast, smart charging, sizing, power generation planning, production commitment

Simultaneous upcoming of photovoltaic generation and electric vehicles increases constraints on electric power system. This paper explores the possible synergy between these players so as to jointly improve the production predictability while ensuring a low carbon mobility. First a context is defined for this collaboration. It consists in the association of a photovoltaic producer and some electric vehicles owners so as to both manage the EV recharge and meet a day-ahead production commitment. Several currently studied questions such as commitment strategies or optimal charging can be transposed into the proposed context called collaborative system. Here, we mainly focus on its sizing in terms of PV rated power and number of vehicles. A simplified model of the system is thus realised, including a day ahead commitment and an optimal vehicle charging planning, based on deterministic vehicle characteristics.First results show a strong influence of the sizing on the potential added value of vehicles in this association. Then, we assess the impact of day-ahead production forecast quality by comparing persistence forecast with some meteorological data. It appears that other things remaining equal, an imprecise forecast will increase the optimal number of vehicles that are supposed to get into the collaborative system. Finally, the robustness of the charging planning is investigated.

Growing concerns around renewable electricity production compel us to look beyond its advantages in terms of low carbon content. We will here focus our attention on photovoltaic (PV) devices but similar investigations can be carried on others energy sources such as wind or ocean waves. As a matter of fact, the electricity production by means of photovoltaic (PV) panels presents a high variability and relatively low predictability. Under such characteristics, the spread of PV plants cannot but be limited as it brings about some additional stress on distribution and transportation networks, while increasing the need for spinning reserves. Canova et al. [START_REF] Canova | Electrical Impact of Photovoltaic Plant in Distributed Network[END_REF] were thus able to show that PV producers could cause voltage fluctuations, fostered by large line impedance and important solar capacity. Then increasing heavily the penetration rate would first require an important strengthening of the grid if no precautionary measure was introduced.

Various proposals are currently being put forward to take into account this poor predictability. For instance, the call for tenders [START_REF]et l'exploitation d'installations de production d'électricité à partir de l'énergie solaire d'une puissance supérieure à 250 kWc[END_REF] which is in force on the French island territories for PV plants above 250 kW c requires that the photovoltaic operator commits himself some time ahead on its production profile. This profile should moreover respect a trapezoidal shape with a period of constant production during mid-day. As a consequence of such constraints, the addition of a storage device which would be coupled to the PV power plant is required. Studies like those of Ru et al. [START_REF] Ru | Storage Size Determination for Grid-Connected Photovoltaic Systems[END_REF] suggest some methods to design a storage capacity associated with PV plants in the context of exchanges with the electrical network. The possibility to take advantage of fluctuations in the electricity price is also taken into account.

Furthermore, the development of electric vehicles (EV) causes additional constraints on the transmission and distribution network. Indeed, these vehicles represent a modification of the total energy demand distribution. Transportations which were previously ensured by liquid fuels are supposed to be more and more electricity seeker. It is therefore necessary to be able to convert more energy into the electricity form and to distribute this energy according to the charging requests of users. Drovtar and Landsberg [START_REF] Drovtar | Large Scale Electric Vehicle Integration and its Impact on the Estonian Power System[END_REF] describe the impact that the emergence of these new consumers can have on a grid, using the Estonian example, just as Turker et al. [START_REF] Turker | Impact of Plug-in Hybrid Electric Vehicles (PHEVs) on the French electric grid[END_REF] do on the French one. It is concluded that in order to integrate a significant share of the vehicle fleet, it is necessary to set up some tariff incentives so as to shift during the night the vehicles charging consumption [START_REF] Rautiainen | Case studies on impacts of plug -in vehicle charging load on the planning of urban electricity distribution networks[END_REF].

However, concurrently to these constraints, electric vehicles do have the potential to bring new services to energy systems. Indeed, they represent a storage capacity which is connected to the electrical network most of the time, as mobility often represent a little share of the vehicle life. Opportunities associated to this scattered storage were highlighted by Kempton and Tomic [START_REF] Kempton | Vehicle-to-grid power implementation: From stabilizing the grid to supporting large-scale renewable energy[END_REF]. First of all, even if the stored energy is limited in comparison with some others grid players, interconnection power and fast response of electric vehicle batteries can enable them to play a part within spinning reserves [START_REF] Han | Development of an optimal vehicle-to-grid aggregator for frequency regulation[END_REF]. But Kempton and Tomic had already considered a particularly fruitful combination of electric vehicles with renewable energy sources. Indeed, the environmental performance of transportation must be considered from well to wheel. As a consequence trying to improve this performance through the use of electric vehicles involves to use some primary energy sources with a low environmental impact. The direct association of electric vehicles with renewable energy sources is therefore very relevant. On the one hand, photovoltaic production -as well as any intermittent production -could use vehicle batteries as a storage device spread out over the grid, so as to mitigate their production variability and improve their predictability. On the other hand, the vehicles would be charged with a production of sustainable energy. Many studies address this collaboration. For instance, Traube et al. [START_REF] Traube | Mitigation of Solar Irradiance Intermittency in Photovoltaic Power Systems With Integrated Electric-Vehicle Charging Functionality[END_REF] investigate the possibility of using vehicles batteries to compensate some fast variations of PV production during cloudy days so as to only send slow power slopes into the grid. Such a compensation will result in a high frequency variation of the charging power that might hasten batteries ageing. Guillou et al. [START_REF] Guillou | Energy management Strategies for optima charging of vehicles with photovoltaic production[END_REF] investigates the collaboration of a PV plant and electric vehicles fleet, but to maximize the self-consumption of the solar production. This minimisation of power exchanges with the grid could be seen as a zero commitment. This paper then proposes as a first step to present a context for the direct collaboration of photovoltaic producers and electric vehicles users, in order to comply with a production commitment which would be communicated one day-ahead to the grid manager (section II). Within this context, different problems must be solved to achieve the operational management of the system. A breakdown of these problems and of their interactions will be proposed. Afterwards, section III will examine more specifically the sizing issue in terms of peak power and number of vehicles. The influence of the quality of the day-ahead production forecast will be specifically investigated.

II. DESCRIPTION OF THE INVESTIGATED SYSTEM

A. A collaborative system Some contexts aiming at structuring the electric vehicles charging in order to allow their integration into the electrical network have already been proposed [START_REF] Kempton | Vehicle-to-grid power implementation: From stabilizing the grid to supporting large-scale renewable energy[END_REF], [START_REF] Guille | A conceptual framework for the vehicle-to-grid (V2G) implementation in Energy Policy[END_REF]. The framework which is here proposed relies on the merging of the usual functions of electric vehicles aggregator and of intermittent electricity producer. The goal is to build an operator who could be considered from a grid manager point of view as a single actor, to some extent similar to a virtual power plant, but with a distributed storage capacity. In this framework, photovoltaic producer is subjected to a constraint of dayahead production commitment. Thus, the manager of the system has to announce today the power profile he might exchange with the grid during the following day. He suffers some penalties according to the gap between the profile he has committed and the profile he actually achieves to exchange. Therefore, he gets associated with some electric vehicles owners who entrust him the control of the battery charging power, subject to the requirement of a fully charged battery at the time they indicated when dropping their vehicle. The commitment constraint does not any more affect the single PV power profile, but is shifted to the power injected into the grid, which is the difference between the PV production and vehicles consumption. The possibility to shed some of the PV power is also taken into account. Such an association between a power plant and vehicles which stands as a mobile storage capacity is afterwards called collaborative system.

It may be noticed that the questions of physical interconnection of vehicles and power plant is not into the scope of this study. That means vehicles can either be in the direct neighbourhood of the PV plant or scattered over a large area but connected through the grid and collaborating with a distant plant. Although these two situations are very different, the management rules are in both cases similar. That is why the theoretical study that is here described does not exclude some remote collaboration, using the electrical network to gather the spread vehicles into a huge and virtual storage capacity. This second possibility does imply a complex infrastructure of communication and measurement, as well as a suitable regulatory framework and the absence of grid constraints. Each of these conditions represents an entire topic by itself and we only focus here on management investigations. Figure 1 presents the main actors in this system.

P pv stands for the maximum power that could be produced by the PV plant at a given instant. A share of this maximum achievable power can be shed P shed if there is no suitable solution to deliver it. Vehicles are overall consumers, but power flow can occasionally inverse in this specific context. Once again, this suppose to use some bidirectional battery chargers. P ev can thus be positive or negative. In the case where several vehicles are involved in the collaborative system, it becomes a vector quantity: P ev = {P ev,1 , ..., P ev,nev } with n ev the number of involved vehicles. A priori, the collaborative system is not supposed to strictly be an energy producer: power exchanged with the grid P grid can be bidirectional. Within this context, the objective is to obtain a power exchanged with the grid which is as close as possible to the commitment profile P * grid taken out the day before. At any time, the system is ruled by the Kirchhoff law:

P pv = P grid + P shed + nev i=1 P ev,i (1)

B. Breaking down of the system management into subproblems

In order to manage a collaborative system, it is necessary to get some predictions of its behaviour, both in terms of maximum power output from the PV panels and of vehicle power needs. The latter is actually a forecast of the energy needed to refill the battery -how much energy ? -and of presence hours of the vehicles -when are they available ? These forecasts are written P X pv and P X ev where X is the value of the prediction horizon. A dayahead forecast is written D+1 and an intraday forecast is written D. Based on day-ahead forecasts, the manager of the collaborative system must then be able to compute its commitment for grid exchanged power P * grid . real-time management of the vehicle charging power depends on both commitment -the objective to be strived towards -and real-time gathered information dealing with the solar irradiance and vehicles user behaviours. Finally a grid operator -a balance responsible party for instance -monitors the difference between commitment and actually achieved profile to determine what should be the payouts and penalties to be applied to the collaborative system. Several problems must be alternately solved to study the system. Figure 2 summarizes these problems and the interactions between them.

We can break down the collaborative system management into three entities :

• strategy of real-time vehicle charging. The goal is to find a way to determine the charging power depending on electricity production, availability and energy needs of vehicles so as to comply as much as possible with the grid commitment. It is therefore necessary to have decided upstream the commitment profile P * grid as well as the collaborative system composition in terms of number of vehicles and peak power.

• commitment strategy. It aims to compute the grid commitment profile for tomorrow based on day-ahead production forecasts P D+1 pv and charging demand forecast P D+1 ev . Knowledge of the system sizing is obviously necessary for such predictions, which can either come from an external specialized service or be produced internally by numerical methods.

• design of collaborative system. It appears than the decision which has to be made far upstream of the previous problems is the sizing of the collaborative system in terms of the peak power and number of involved vehicles.

Each one of the aforementioned issues is a complex question which is currently being investigated in other contexts and other systems. Among others, Peralta et al. [START_REF] Peralta | Unit Commitment with Load Uncertainty by Joint Chance-Constrained Programming[END_REF] propose to determine production commitments within a framework where the consumption of domestic users is uncertain, using a joint chance-constrained programming method. Such an approach could be transposed to a situation where uncertainty is shifted from consumption to production. Furthermore, Vaya and Andersson [START_REF] Vaya | Integrating renewable energy forecast uncertainty in smart-charging approaches for plug-in electric vehicles[END_REF]develop commitment strategies and vehicles charging programs that will allow to compensate forecasting errors of wind generation. In there study, the objective is to move the complexity towards the commitment strategy so as to get the real-time charging strategy simpler. On the contrary, Sortomme and El-Sharkawi [START_REF] Sortomme | Optimal Charging Strategies for Unidirectional Vehicle-to-Grid[END_REF] carry forward the complex task on the instantaneous charging management with an additional coupling with energy spot price. The model which describes vehicles is then deterministic.

C. Scope of the further investigation

This section has presented a context which involves electric vehicles and intermittent electricity producer into a collaborative system. In addition, a breakdown of issues that are to be solved in order to manage such a system 2. Outline of sub-problems occurring along the collaborative system management and connections among them.

was presented. One could notice that several ongoing works consistently fit into the suggested context. The next section of this article tackles more specifically the issue of collaborative system sizing. As stated earlier, it is the very first problem on which both depend commitment and instantaneous charging power strategies.

III. SYSTEM SIZING

As mentioned in Section II, the design of the collaborative systems is crucial to build commitment and real-time charging strategies. The term design here stands for choosing the rated peak power of the PV system P peak and the number of cars it will collaborate with n ev , which means the equivalent storage capacity as well as the interconnection power. However, to achieve this design study, it is necessary to get a relevant model which will take into account the entire system. One can hardly consider to model some of the complex management strategies that can be implemented for solving dayahead (D-1) commitment and real-time charging issues. Simplifying assumptions must be made so as to build a model which takes into account the whole system while maintaining a low level of complexity, consistent with the exploration of many design possibilities.

A. System modelling and assumptions

Most of the vehicles are used for a daily journey. Since the state of energy of an electric vehicle when it is dropped into the collaborative system mainly depends on its utilisation for mobility, it does not seem relevant to intent to develop charging strategies with an horizon farther than the coming evening. This model thus only considers optimisation throughout a day.

In order to determine the number of vehicles that might be associated to a PV plant of a given rated power P peak , we seek to evaluate the added value of a given number of vehicles for the collaborative system. This added value is basically the difference between how much vehicles cost to the system and how much they make.

• how much vehicles cost: the energy needed to charge vehicles could otherwise had been injected into the grid. So as to assess the vehicles potential in terms of improving the overall system predictability, we consider that the daily charging of vehicles is ensured for free. The income loss is thus C out = n ev .a.(E ev -E ev , where E 0 ev the initial battery energy, E ev the rated battery energy and a the selling price to the grid manager. Here we assume a = 0.20 e/kW h. This is a quite • how much vehicles make: the control of vehicles charging power allows the system to strive towards its commitment power profile P * grid and thus avoid some penalties. Income is then the following reduction of penalties:

C inc = C 0 pen - C nev
pen where C 0 pen is the penalty amount for a single PV plant and C nev pen for a collaborative system with n ev vehicles.

Our goal will thus be for the rest of the study to assess the added value of n ev vehicle within a collaborative system. This added value AV is obviously the difference:

AV = C inc -C out (2)
Computing penalties requires to choose a penalisation law. We adopt the shapes presented Figure 3. This law has as a parameter a tolerance threshold ε, expressed as a fraction of the plant peak power. Below this threshold, commitment gap ε grid = P grid -P * grid is not penalized. Above, a linearly increasing penalty is applied to the exchanged power. An error corresponding to twice the threshold means that the penalty and the selling price are equal. I.e. if the commitment gap is greater than 2ε, injecting power into the grid would be a burden for the collaborative system.

We then need to proceed to a phenomenological modelling of the day-ahead commitment strategy and the instantaneous charging strategy. During a real operational management of the collaborative system, they would both use some complex methods and real-time information. As we here intent to achieve a sizing study, we will use time series of past records describing the production of a PV plant and the forecast of its production. The considered asset is a 2.64 M W plant in Corsica island (France), monitored from May 2012 to April 2013 with a ∆T = 15 min time step. Corresponding meteorological data are also available.

1) Day ahead commitment modelling:

The first concern is the day-ahead commitment strategy. It has to be computed on the base of a production forecast P D+1 pv and a prediction of the charging power profile of the vehicles P D+1 ev . These forecasts are the ones available one dayahead (D + 1).

Assumption: deterministic behaviour of electric vehicles. All vehicles are considered identical and have the same driving patterns. Their individual storage capacity is E ev = 25 kW h, the interconnection rated power is P ev = 3 kW . They are assumed to be available for the collaborative system between H start = 9 a.m. and H end = 6 p.m. and have an initial state of energy E 0 ev = 10 kW h. While this assumption may seem very restrictive, it roughly fits to the macroscopic behaviour of many fleets obeying working hours [START_REF] Sortomme | Optimal Charging Strategies for Unidirectional Vehicle-to-Grid[END_REF].

The P D+1 ev forecast is thus a constant power profile from H start to H end which value is n ev . E ev -E 0 ev Hend-Hstart . The P D+1 pv is in a first attempt done by a persistence forecast: the profile measured today is the forecast for tomorrow. Paragraph III-D will investigate the impact of this forecast technique. Grid commitment is then:

P * grid = P D+1 pv -P D+1 ev (3)
2) instantaneous charging modelling: The next concern is the instantaneous charging strategy. An operational strategy might be based on tools like short terms models of the PV production, stochastic approaches or dynamic programming so as to take at each instant the best decision. Here, as we have to keep a low computational cost, it is not possible to implement such a real time strategy. We would instead try to determine a charging planning for the entire day, thanks to the PV production forecast P D pv which is available at the morning of the considered day.

Assumption: in a first attempt, it is supposed that P D pv is perfect. From the morning of a considered day, a perfect forecast of the coming production is available until the evening. This perfect prediction assumption from the morning may seem very strong. The paragraph III-E will aim to study how robust the chosen charging strategy is with respect to the forecast error.

B. Description of an optimisation problem

Under the previously described modelling, the charging strategy can be computed as a linear optimisation problem: min Pev,Pshed C nev P en (P * grid -P grid )

such that

E ev (t) = E 0 ev + ∆T t Hstart P ev (τ ) (5) ∀t ∈ [H start , H end ],
-n ev .P ev ≤ P ev (t) ≤ n ev .P ev (6)

0 ≤ E ev (t) ≤ n ev .E ev (7) E ev (H end ) ≥ 90%n ev .E ev (8) ∀t / ∈ [H start , H end ], P ev (t) = 0 (9)
I.e. the aim is to find the best charging planning P ev and production shedding planning P shed for the considered day so as to minimize penalties coming from the commitment gap. Constraints are the dynamic equation of the storage (5), maximum interconnection power of vehicles [START_REF] Rautiainen | Case studies on impacts of plug -in vehicle charging load on the planning of urban electricity distribution networks[END_REF], extreme battery states of energy [START_REF] Kempton | Vehicle-to-grid power implementation: From stabilizing the grid to supporting large-scale renewable energy[END_REF], charged battery at the outgoing hour (above 90% of full charge) (8) and zero power exchanged while batteries are not connected [START_REF] Traube | Mitigation of Solar Irradiance Intermittency in Photovoltaic Power Systems With Integrated Electric-Vehicle Charging Functionality[END_REF].

It may be noticed that the presented formulation is generic and could in the future be used to take into account some other system modelling including for instance random vehicles availability or imperfect forecast of future production. In the deterministic form which is here adopted, this optimization problem is solved with an interior point algorithm. Optimisation results over a 4 day sample is shown Figure 4. The time step is 15 min. This sample illustrates several different situations. During the first and fourth days, forecast has both underestimated the production and made temporary overestimations. Thus power shedding is used and batteries are punctually discharged to compensate the production lack. The second day is when the production has been quite well predicted by persistence. Thus the charging planning is very close to its day-ahead forecast P D+1 ev , with a rectangular shape. Finally, the third day is a heavy overestimation case, when the grid commitment can not be fulfilled. Overall adequacy to the commitment profile is thus improved and therefore penalties lowered using vehicles charging power control. However, introducing vehicles into the system results in a diminution of the energy injected into the electrical network, and thus income losses. An optimum participation of electric vehicles is likely to exist.

C. Sizing results

The design of a collaborative system such as we here consider requires a priori that we determine the value of the PV plant peak power and the number of vehicles involved into the association. Because of the way the problem has been set and of its characteristics, it can be noticed this design is only based on the ratio of these two quantities. Therefore the resolution of the previous optimization problem is iterated over a one year period and for different values of the n ev .E ev /P peak ratio. The results depend on many parameters, not limited to the used time series describing the photovoltaic production and forecast, the supposed behaviour for vehicles, as well as the modelling assumptions. That is why the aim of the here presented results is not to propose predictions of profitability values, but rather to point out some tendencies that may be similar whatever the renewable sources and the storage system. Results are shown Figure 5. The impact of penalisation law is investigated for three different values of the tolerance threshold, as represented Figure 3. The y-axis represent the added value AV = C inc -C out over a one year period divided by the peak power P peak . It is thus an annual benefit or loss by photovoltaic watt.

These results are indeed highlighting two very different sizing areas. In a first sizing region, there are not enough vehicles to make the most from the potential of the collaboration. They do not represent a sufficient storage capacity to absorb the production forecast errors. On the contrary, if the vehicles are too numerous, an important share of the produced energy is dedicated to the batteries charging and thus is directly consumed. The residue that still can be sold to the grid manager is then no longer sufficient to compensate for the shortfall, even though compliance with the commitment were excellent. Two significant points can be discussed: the maximum profitability and the zero profitability. The position of these points necessarily depends on how the system is ruled by the penalisation law. On the first hand, if the law is stringent with a low tolerance threshold, collaboration is very relevant as a small storage capacity is able to save a large amount of penalties. On the other hand, if the penalty law is more tolerant, it becomes difficult to compensate for the energy absorbed by vehicles charging. In the median case, where the threshold is set at 5% of the peak power, sizing with maximum profitability is 1.7W h/W p , i.e. in a 2.6 M W PV plant as the one used for production records, 180 vehicles with a 25 kW h capacity. This represents a total storage capacity of 4.5 M W h or a roughly 2 hours at P peak . The energy needed for the daily charging of these vehicles is 2.7 M W h which means an hour of full power production. Concerning the zero profitability sizing point of 4.2 W h/W p , it is very close to a total self consumption of the entire production by vehicles.

We will now question in the two coming paragraphs the influence of the two forecast techniques used so far. The day-ahead production forecast P D+1 pv has been made by persistence and the forecast for the current day from the morning P D pv was a perfect prediction.

D. Influence of the day-ahead production forecast precision

So as to assess the influence of the P D+1 pv day-ahead production forecast quality on the optimal sizing of the collaborative system, a comparison is carried between the three following cases, all other things remaining equal:

• the day-ahead production forecast is obtained by persistence (section III-C).

• the day-ahead production forecast is deduced from irradiance forecast which has been computed by Meteo France meteorologists, using the AROME numerical weather model.

• the time series used for day-ahead production forecast is a perfect prediction.

In particular, the forecast for the current day P D pv is still in this paragraph a perfect prediction. For the considered PV plant in Corsica, while persistence forecast has a mean squared error equal to 25% of the mean power output, the forecast coming from meteorologists lowers the error to 15% of the mean output. Figure 6 draws the comparison between computed profitability of the collaborative system for the three described cases.

It therefore appears that a better prediction greatly reduces the predicted added value of vehicles into a collaborative system. Indeed, if the manager of the plant can use a better forecast, he is thus able to compute a safer commitment. Penalties related to the commitment gap will then be necessarily lowered. Having a storage capacity to compensate for this gap will be less useful than in the case where the prediction and the realisation are very different. Added value is therefore generally reduced and the optimum sizing is shifted to a smaller number of vehicles, as it represents smaller energy need and a sufficient capacity storage regarding Impact of the day-ahead production forecast quality over the vehicles added value. The day-ahead commitment profile is computed using either a persistence forecast, a second one done by meteorologists and a perfect prediction.

to low forecast errors. The extreme case is a perfect day-ahead production forecast. In such a situation, it is theoretically possible to get committed without any error. Providing a storage capacity to compensate for commitment gaps is no longer of any use. Moreover, as this storage capacity is also an overall energy consumer, it becomes a burden. The profitability of the collaborative system is then necessarily negative since it is always mandatory to charge the vehicles even when controlling their charging power is no longer useful. Profitability then decreases proportionally to the number of vehicles.

We here only focus on the system predictability issue. That is why if perfect forecast are available, storage device is of no interest. If we had tackled also the PV variability, a storage would have always been useful, even in perfect forecast situations.

It may be noticed that for heavily oversized systems, each one of the three investigated forecast techniques reaches the same asymptotic behaviour: a linearly decreasing profitability. It means that when the vehicles storage capacity is sufficient, the commitment gap has already been perfectly mitigated. There is thus no more penalty to save. But vehicles still cost some energy to the collaborative system.

E. Robustness of the charging strategy with respect to current day forecast error

We aims in this last paragraph to assess the consequences of the assumption of production perfect forecast from the morning and for the entire current day P D pv . Whereas the charging power planning was previously computed with a perfect production prediction for the current day, we will now once again use a persistence forecast. The charging power profile is thus optimised on the basis of a forecast and not any more of the real achievement. It is then blindly applied to the real production profile. As a consequence, the profile of power exchanged with the grid is not the one guessed while optimising. This sensitivity study is made possible by the deterministic context. Even if the blindly applied charging planning does not suit the real PV production, it will not cause some inconsistency such as for instance extreme state of energy violations.

As the persistence forecast is here used for the intraday forecast P D pv , the grid commitment strategy, based on a day-ahead forecast P D-1 pv , has to be computed with another prediction technique. Moreover, the P D-1 pv should be of a poorer quality compared to P D pv , as available information is supposed to increase as the clock is ticking. That is why a two days ahead persistence is here pv . The production profile measured today becomes the forecast for the day after tomorrow. Figure 7 illustrates this modification over the same sample as figure 4.

The profile of power exchanged with the grid is here necessarily downgraded in comparison with the perfect forecast case, regarding to commitment fulfilment. The degradation of performance is however highly dependent on the considered PV plant location. Corsica has a significant share of very sunny days when production forecast is very easy to achieve using a clear sky model. This characteristic reduces the impact of poor quality forecast. Performance is also dependent on the chosen modelling of the system. Figure 8 presents variations of vehicles added value depending on the forecast technique.

As vehicles charging is not optimal towards the real production profile, their added value is heavily lowered and the optimal sizing is shifted to fewer vehicles.

IV. CONCLUSION

In a first time, this study has presented a context for the collaboration between electric vehicles and photovoltaic plants. The management of such a collaborative system requires to solve several problems: a production commitment issue and a real-time charging strategy issue. Several independent and current research studies can fit into the context proposed here. In order to highlight the influence of the system sizing in terms of peak power and of electric vehicles number, a simplified modelling of the collaborative system behaviour is carried out. Under these assumptions, it is possible to show that the value of the ratio n ev E ev /P peak -number of vehicles over PV peak power -set the overall system profitability. The influence of the production forecast quality over the profitability has been investigated by a comparison between performance achieved on a meteorologist based forecast and a persistence based forecast. It appeared that a poor quality day ahead forecast increases the added value of vehicles into the collaborative system. Indeed, the commitment gap -that has to be mitigated thanks to vehicles charging power planning -is more important.

The sizing which has been obtained under several modelling assumptions may be subject to further investigations. The impact of many of the parameters describing vehicles behaviour have not been been studied here. Moreover, it is possible to take into account more complex management rules for commitment and real time control and then to quantify optimal sizing modifications it would imply. Furthermore, a real time modelling would enable to implement some stochastic vehicles behaviours.

The strongest assumption made throughout this study aiming at a sizing proposal is that a perfect knowledge of the coming production is available from the morning and for the current day. It has be checked out that the profitability of collaborative system is strongly affected when production is not the one that was expected. Attempting to establish some robust charging policies regarding to the forecast error is thus one of the main investigation paths to yield a better performance.

V. ACKNOWLEDGEMENT This study has been made possible thanks to the support of Langa Solar for the monitoring of its PV plant and Frank Baraer (Météo France) for the extraction of the corresponding meteorological data.

Figure 1 .

 1 Figure 1. Outline schematic of the collaborative system players and main parameters.

Figure 3 .

 3 Figure 3. Penalisation law CP en for three tolerance threshold values: 3%P peak , 5%P peak and 10%P peak

Figure 4 .

 4 Figure 4. 4 day sample for a {P peak = 2.64 M Wp, nev = 184} collaborative system.

Figure 5 .

 5 Figure 5. Evolution of the collaborative system profitability against nev.E ev /P peak ratio, over one year and for different tolerance thresholds.

  Figure 6.Impact of the day-ahead production forecast quality over the vehicles added value. The day-ahead commitment profile is computed using either a persistence forecast, a second one done by meteorologists and a perfect prediction.

Figure 7 .

 7 Figure 7. 4 day sample for a P peak = 2.64 M W and nev = 184 collaborative system. Impact of a persistence forecast for the current day production.

Figure 8 .

 8 Figure 8. Impact of the intraday forecast quality over the vehicles added value.