S Carpy

R Manceau
email: remi.manceau@lea.univ-poitiers.fr

Turbulence modelling of statistically periodic flows: synthetic jet into quiescent air

Keywords: Turbulence modelling, URANS, Synthetic Jet, Statistical Unsteadiness, Control

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

Nomenclature

.

Long-time averaging operator

< . >

Phase averaging operator

u * i
Instantaneous velocity p * Instantaneous pressure

U i = u * i
Long-time-averaged velocity

U i =< u * i >
Phase-averaged velocity Although statistical turbulence modelling (RANS) represents the current industrial standard, there is a considerable interest in the possibility of obtaining informations on the large-scale unsteadiness. As long as Large Eddy Simulation will remain too expensive for overnight simulations, some room exists for the development of less expensive methods able to predict the very large scales of the flow. This type of approaches has important applications in different domains: among others, aeronautics [START_REF] Gatski | DNS/LES for NASA aerodynamic needs and engineering applications[END_REF][START_REF] Jin | Two-equation turbulence model for unsteady separated flows around airfoils[END_REF], control [START_REF] Getin | Simulation numérique du contrôle actif par jet pulsés de l'écoulement turbulent autour d'un cylindre circulaire[END_REF][START_REF] Hassan | Effects of zero-mass synthetic jets on the aerodynamics of the NACA-0012 airfoil[END_REF], aeroacoustics [START_REF] Bastin | Computation of jet mixing noise due to coherent structures: the plane jet case[END_REF], power generation [START_REF] Benhamadouche | LES, coarse LES, and transient RANS. Comparisons on the flow across a tube bundle[END_REF] and environment [START_REF] Kenjereš | Contribution to elliptic relaxation modelling of turbulent natural and mixed convection[END_REF][START_REF] Kenjereš | Combined effets of terrain orography and thermal stratification on pollutant distribution in a town valley: a T-RANS simulation[END_REF]) can be cited.

u ′′ i = u * i -U i Fluctuation around phase-averaged veloci- ties < u ′′ i u ′′ j > Reynolds stresses k ′′ = 1 2 < u ′′ i u ′′ i >
In many statistically steady flows, as soon as the time derivatives are included in the equations, unsteady RANS solutions appear naturally [START_REF] Johansson | Numerical simulation of vortex shedding past triangular cylinders at high Reynolds number using a k-ε turbulence model[END_REF][START_REF] Bosch | Simulation of vortex shedding past a square cylinder near a wall[END_REF][START_REF] Benhamadouche | LES, coarse LES, and transient RANS. Comparisons on the flow across a tube bundle[END_REF][START_REF] Jin | Two-equation turbulence model for unsteady separated flows around airfoils[END_REF][START_REF] Lasher | On the computation of turbulent backstep flow[END_REF][START_REF] Iaccarino | Unsteady 3D RANS simulations using the v 2 -f model[END_REF]. In such cases, long-timeaveraged solutions are often better than solutions obtained by steady-state computations [START_REF] Durbin | Separated flow computations with the k-ε-v 2 model[END_REF][START_REF] Iaccarino | Unsteady 3D RANS simulations using the v 2 -f model[END_REF]. However, there is no clear consensus on the definition of the mean operator for this type of computations: the definition of the turbulent scales resolved by the method is not a priori known. Recent approaches, like Semi-Deterministic Modelling (SDM: Ha Minh and [START_REF] Ha Minh | Semi-deterministic modelling for flows dominated by strong organized structures[END_REF][START_REF] Aubrun | The semi-deterministic approach as way to study coherent structures. Case of a turbulent flow behind a backward-facing step[END_REF][START_REF] Bastin | Computation of jet mixing noise due to coherent structures: the plane jet case[END_REF], Very-Large-Eddy Simulation (VLES: [START_REF] Speziale | Turbulence modeling for time-dependent RANS and VLES: a review[END_REF], Detached-Eddy Simulation (DES: [START_REF] Spalart | Comments on the feasibility of LES for wings, and on a hybrid RANS/LES approach[END_REF], Limited Scales simulations (LNS: [START_REF] Batten | LNS-An approach towards embedded LES[END_REF], Scale-Adaptative Simulation (SAS: [START_REF] Menter | SAS turbulence modelling of technical flows[END_REF], Partially-averaged Navier-Stokes modelling (PANS: [START_REF] Girimaji | PANS turbulence model for seamless transition between RANS and LES: Fixed-point analysis and preliminary results[END_REF] or Partially-integrated Transport modelling (PITM: [START_REF] Chaouat | A new partially integrated transport model for subgrid-scale stresses and dissipation rate for turbulent developing flows[END_REF][START_REF] Schiestel | Towards a new partially integrated transport model for coarse grid and unsteady turbulent flow simulations[END_REF] are designed in order to ensure unsteady solutions.

In the present study, the flow belongs to the class of statistically periodic flows, i.e., flows in which the boundary conditions are periodic in time. The use of RANS models does not raise any theoretical issue in this case, since Reynolds average is equivalent to phase average. However, obtaining a periodic solution is not guaranteed: some frequencies can appear that are not harmonics of the forcing frequency. Moreover, while the relative performance of the numerous turbulence models available have been investigated extensively in statistically steady flows, their ability to close successfully the phase-averaged equations in periodic flows have received much less attention, and generally for modified versions of the models (Tardu and Da Costa, 2005;Stawiarski and Hanjalić, 2002;[START_REF] Bremhorst | Two-time-scale turbulence modelling of a fully-pulsed axisymmetric air jet[END_REF][START_REF] Lardeau | Unsteady RANS modelling of wake-blade interaction: computational requirements and limitations[END_REF]. Therefore, a workshop dedicated to the evaluation of modelling strategies in statistically periodic flows was recently organized [START_REF] Rumsey | Summary of 2004 CFD validation workshop on synthetic jets and turbulent separation control[END_REF][START_REF] Sellers | Langley Research Center workshop on CFD validation of synthetic jets and turbulent separation control[END_REF]. The present study focuses on the comparison of the performance of first and second moment closures: the standard k-ε model (Laun-der and Spalding, 1974) and the Rotta+IP second moment closure [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF] are used, without any case-specific modification. Standard wall functions are applied since near-wall regions do not play a significant role.

The flow studied in the present paper is one of the test cases selected for the workshop. The flow consists in an oscillatory jet with zero-net-mass flow issuing into quiescent air. This flow is of high industrial interest, in particular in the frame of turbulence control, since it is representative of MEMS piezoelectric technology, as developed by [START_REF] Seifert | Oscillatory blowing. a tool to delay boundary-layer separation[END_REF] and [START_REF] Amitay | Modification of the aerodynamics characteristics of an unconventional airfoil using synthetic jet actuators[END_REF].

Separation control through the use of suction, blowing and periodic excitation has lead to a number of studies (see, e.g., the review by [START_REF] Greenblatt | Control of flow separation by periodic excitation[END_REF]. The present case was designed in order to evaluate the capabilities of CFD and modelling techniques in this type of problems, concentrating on the synthetic jet region, without interaction with an outer flow. The experiments used as reference were performed at NASA LaRC by Yao et al. (2004) and are available on the web site of the workshop [START_REF] Sellers | Langley Research Center workshop on CFD validation of synthetic jets and turbulent separation control[END_REF].

In such a problem, it is convenient to introduce two averaging operators: longtime averaging, denoted by . and phase (or ensemble) averaging, denoted by < . >. Instantaneous quantities are denoted by stars (u * i , p * ) and the following definitions are used:

Long-time-averaged velocity:

U i = u * i
Phase-averaged velocity:

U i =< u * i >
In order to avoid confusion with usual notation in RANS modelling, where u ′ i stand for the total fluctuations (around the long-time-averaged velocities), the fluctuations around phase-averaged velocities are denoted by u ′′ i :

u ′′ i = u * i -< u * i > (1)
Thus, in the present periodic case, the transport equations to be solved are those of the phase-averaged velocity components U i :

∂ U i ∂t + U j ∂ U i ∂x j = - 1 ρ ∂ P ∂x i + ν ∂ 2 U i ∂x j ∂x j - ∂ < u ′′ i u ′′ j > ∂x j (2)
and those of the Reynolds stresses

< u ′′ i u ′′ j > (or of k ′′ = 1 2 < u ′′ i u ′′ i >
in the frame of the eddy-viscosity model), and of ε ′′ , the dissipation rate of k ′′ : all these equations are formally identical to the standard steady RANS equations, except for the presence of the time derivative.

In order to close the equations, two different standard turbulence models are used: the standard k-ε model [START_REF] Launder | The numerical computation of turbulent flows[END_REF], denoted k ′′ -ε ′′ in the present case, and the Rotta+IP Reynolds-stress model [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF].

Description of the geometry and available experimental data

The oscillating jet of air emanates from a slot in the floor. The slot is h = 1.27 mm wide and 28 h long, in order to produce a 2-D plane jet. The slot is at the bottom of a cubic enclosed box 480h = 609.6 mm large, at the centre of the floor and parallel to the outer walls. The geometry and the coordinate system are shown in Fig. 1.

The jet is produced by a piezo-electric diaphragm on the side of a narrow cavity under the floor, which is driven at 444.7 Hz [START_REF] Sellers | Langley Research Center workshop on CFD validation of synthetic jets and turbulent separation control[END_REF]. The maximum jet velocity generated could reach 28 m s -1 .

Although the actuator operating parameters including diaphragm displacement, internal cavity pressure, and internal cavity temperature were monitored to provide boundary conditions for CFD simulations, we chose not to solve the flow in the cavity: the present work focuses on the performance of turbulence models, rather than on the other difficult issue of representing a vibrating diaphragm with simple boundary conditions.

Computational results in the present paper are compared with the PIV measurements performed by Yao et al. (2004). Phase-averaged (every 5 deg) and long-timeaveraged velocities and Reynolds stresses (only components in the x-y plane) are available in the window -2.15 < x/h < 3.1, 0.1 < y/h < 6.6, z = 0, corresponding, in actual size, to -2.73 mm < x < 3.94 mm, 0.127 mm < y < 8.38 mm.

Numerical method and sensitivity studies

Computations are performed using Code Saturne, a finite volume solver on unstructured grids developed at EDF [START_REF] Archambeau | Code Saturne: A finite volume code for the computation of turbulent incompressible flows -Industrial applications[END_REF], for vectorial and parallel computing.

The code solves time-dependent RANS or LES equations for Newtonian incompressible flows. Space discretization is based on the collocation of all the variables at the centre of gravity of the cells. Velocity/pressure coupling is ensured by the SIM-PLEC algorithm, with the Rhie & Chow interpolation in the pressure-correction step. The Poisson equation is solved using a conjugate gradient method, with diagonal preconditioning. The present study is performed using the usual standard k-ε model [START_REF] Launder | The numerical computation of turbulent flows[END_REF] and the Rotta+IP second moment closure [START_REF] Launder | Progress in the development of a Reynolds-stress turbulence closure[END_REF]. Central differencing is used for convection. For timemarching, a second order Crank-Nicholson scheme is used for the k-ε model and, for stability reasons, a first order Euler implicit scheme for the Reynolds-stress model. However, as shown below, time-step-independent solutions are reached in all cases. Computations are performed on 8 processors of a PC cluster.

Boundary conditions

Since the flow inside the cavity is not solved, unsteady Dirichlet boundary conditions are specified at the slot exit (Fig. 2), which are functions of x and phase: profiles at y/h = 0.079 of the available phase-averaged velocities U i and Reynolds stresses < u ′′ i u ′′ j > are extracted from PIV data and used as inlet conditions. The missing component < w ′′ 2 > is reconstructed by assuming < w ′′ 2 >=< u ′′ 2 >.

Concerning the dissipation rate ε ′′ , prescribing an inlet value is already an issue. Indeed, ε ′′ obviously varies with time, as well as all the other usual turbulent scales (time scale, length scale, eddy-viscosity). A possible way is to arbitrarily assume that one of them is constant, e.g., the length scale L. A periodic ε ′′ can thus be obtained from ε ′′ = k ′′ 3/2 /L. Preliminary computations (not shown here) have indicated that a value of L of the order of magnitude of what can be expected in a fully developed turbulent channel flow (typically one-tenth of the half-width of the slot) gives too much turbulent dissipation. Therefore, the value of the length scale prescribed in the computations presented here is half the slot width. The necessity of using such a large value is to be linked to the fact that, despite a turbulence rate of 5%, the flow issuing from the cavity is far from being a fully developed turbulent flow: the Reynolds number based on the slot width and the maximum inlet velocity is only 2400, and the geometry of the cavity does not allow the turbulence to reach a spectral equilibrium before issuing from the slot.

Solid-wall boundary conditions cannot be used at the remaining boundaries of the domain, since the fluid is considered incompressible: when fluid is blown from the slot, it is necessary to enable an equivalent volume of fluid to go out of the calculation domain in order to satisfy incompressibility. Since the experimental box is huge compared to the size of the slot (480h), a smaller domain was used, as shown in Fig. 1. Outlet boundary conditions (homogeneous Neuman conditions) are applied at the open boundaries of the domain: as shown in section 3.3, the results are not influenced by this restriction of the domain. However, in this type of flows, fluid enters the domain through the "outlet boundaries" in different regions, at different periods of the cycle, in particular during suction: everywhere the mass flux enters the domain, Dirichlet boundary conditions are used for turbulent quantities, using values identical to those used for the initial conditions described just below.

Initial conditions

Computations are started from rest state ( U i = 0), with a low residual turbulent energy (k ′′1/2 / V max =0.05%) and a dissipation rate corresponding to the ratio ν t /ν = 10. Computations are directly started with central differencing with large time steps (∆t corresponding to ∆φ = 2.5 deg), which shows the stability of the code (for convenience, time steps are given in phase degrees ∆φ throughout the paper, i.e., the value corresponds to ∆φ = 360f ∆t, where f is the frequency of the jet, such that ∆φ = 360 deg corresponds to a cycle). When a satisfactory periodic state is reached (after about 15 cycles), the time step is reduced and the computation is run until a new periodic state is obtained. Note that periodicity was only evaluated in the region where the flow is of interest (basically the region in which the results can be compared with experiments): indeed, it would take hundreds of cycles to reach a periodic solution in the entire domain.

Sensitivity studies

In order to investigate the influence of the size of the domain, the time-step and the grid, seven different computations were performed with the standard k ′′ -ε ′′ model and the Rotta+IP second moment closure. Conclusions are given below.

Influence of the domain size

k ′′ -ε ′′ computations were performed using two different domain sizes (small and large, Fig. 1) with the same grid resolution (i.e., for the larger domain, cells were added around the grid used for the smaller domain). The grid covering the smaller domain is shown in Fig. 3a. Computations were carried out with the same time step corresponding to ∆φ = 2.5 deg, i.e., 144 time steps per cycle. Fig. 4a clearly shows that the restriction of the domain has no influence on the solution. Therefore, the smaller domain in used in the following of the paper.

Influence of the time step

The influence of the time-step on the solution has been carefully investigated. Using the coarse grid in the smaller domain, the time step was progressively reduced: for the k ′′ -ε ′′ model (second-order time-marching scheme), the solution becomes independent of the time step below a value of 0.5 deg (see Fig. 4b). This conclusion was confirmed using the fine grid. For the Reynolds-stress model (first-order timemarching scheme), the time-step independence is reached at 0.25 deg, i.e., 1440 time steps per cycle. It is worth pointing out that reaching a time-step-independent solution requires a large number of time steps per cycle. In a recent study, [START_REF] Lardeau | Unsteady RANS modelling of wake-blade interaction: computational requirements and limitations[END_REF] showed that the solution obtained by URANS in the case of a flow around a cylinder is qualitatively modified (from a quasi-steady solution to a vortex street) when the time step is reduced: at least 125 time steps per cycle are necessary to properly obtain the vortex shedding. The present study indicates that a true time-step independence is actually much more demanding.

Grid refinement

Three different grids (coarse grid: 32900 cells, middle grid: 47610 cells, fine grid: 130950 cells), shown in Fig. 3, are used. Fig. 5 shows that the influence of the grid refinement is small (the discrepancies between the middle grid and the fine grid are of the order of 1%). Therefore, all the results presented in the following of the article were obtained using the middle grid.

Results and discussion

Global characteristics of the flow

Figs. 7 and8 show comparisons of the isocontours of the vertical phase-averaged velocity V from PIV measurements, k ′′ -ε ′′ and RSM computations, for both blowing and suction phases of the cycle (see Fig. 2). It can be seen that the global topology of the phase-averaged flow is well reproduced by both models. It is first worth pointing out that the turbulence models play the role that is expected in this type of flow: the resolved velocity field V i is free from turbulent eddies and V i only represents the oscillating motion induced by the diaphragm.

The pair of large vortices (Fig. 7) that appears at the beginning of the blowing phase is directly inherited from inviscid dynamics: it is induced by the vorticity injected into the domain through the slot (see section 4.4). However, even though the turbulence models do not play a role in their creation, the turbulence field has a significant influence on their evolution into the domain, and the Reynolds-stress model produces a much more realistic picture. The k ′′ -ε ′′ model is much too diffusive as can be observed in Fig. 6. The vortices are dissipated very rapidly, as seen in the movies that can be found on a dedicated website [START_REF] Carpy | Movies associated to the paper: Turbulence modelling of statistically periodic flows: synthetic jet into quiescent air[END_REF]. Fig. 6 shows profiles in the symmetry plane of the phase-averaged velocity V . It appears that the k ′′ -ε ′′ model predicts a much too rapid decrease of the peak velocity, and, in particular, is unable to reproduce the acceleration above the slot which is observed in the PIV data: the velocity at the slot exit reaches its maximum 28.3 m s -1 at phase 45 (see Fig. 2) and, at phase 90, the peak velocity at about y = 4 mm (3.15h) is close to 35 m s -1 . The "fluid particle" with the max-imum velocity at phase 90 has been experiencing an acceleration since it issued from the slot. Fig. 6 shows that this behaviour is reproduced, but underestimated, by the Reynolds-stress model. In both cases, the convection velocity of the peak is underestimated (Figs. 6 and7). The misprediction of the peak velocity and of the convection velocity are both due to the overestimation of the turbulence intensity, which is particularly severe with the k ′′ -ε ′′ model, as will be shown in following sections.

During the suction phase of the cycle, Fig. 8 shows that the fluid is mainly sucked from the sides. This is the main explanation (apart from experimental errors) of the fact that, in Fig. 2, the integral of the curve is not zero: what is plotted is the velocity at the centre of the slot, at the elevation y = 0.08h. Since fluid is mainly blown at the centre of the slot and sucked from the sides, a very significant part of the fluid that enters the cavity during the suction phase is not seen at that point. Therefore, computations cannot be performed using uniform profiles (independent of x): if the flow in the internal cavity is not solved, it is necessary to base the inlet conditions on profiles extracted from the PIV velocity fields, as is done in the present study.

This characteristic behaviour of the flow is confirmed by Fig. 9, in which long-timeaveraged velocity fields are shown. It clearly appears that the fluid tends to issue from the cavity at the centre of the slot and to enter from the sides. Both models (eddy-viscosity and Reynolds-stress models) correctly reproduce this behaviour.

Momentum budget

In order to investigate the discrepancies between the two models, it is convenient to focus on the material derivative of the velocity in the symmetry plane x = 0:

d V dt = - 1 ρ ∂ P ∂y - ∂ < v ′′2 > ∂y - ∂ < u ′′ v ′′ > ∂x + ν ∂ 2 V ∂x k ∂x k (3)
The acceleration of the fluid particle described above results from a competition between the first three terms on the right hand side (pressure gradient and Reynolds stresses), since the molecular diffusion is negligible. The main difference between the two models obviously lies in the Reynolds-stress tensor, and in particular in the prediction of the sign of the contribution of normal stress < v ′′2 >: the shear stress < u ′′ v ′′ > is globally well reproduced by both models, and always contributes as a sink in Eq. (3).

Figs. 10 and 14 show that the Reynolds-stress model predicts a region of negative production right above the slot: in the present 2-D case, using the incompressibility condition, the production of k ′′ reduces to:

P k ′′ =< u ′′2 > ∂V ∂y -< v ′′2 > ∂V ∂y -< u ′′ v ′′ > ∂U ∂y + ∂V ∂x (4)
In the symmetry plane, the last term is zero because < u ′′ v ′′ >= 0. Fig. 10 shows that the first term is positive and the second term is negative. This is a direct consequence of the unsteadiness of the boundary conditions. Indeed, since the velocity at the slot exit decreases rapidly after phase 45 (Fig. 11a→b), the partial derivative ∂V /∂y at phases > 45 deg is positive between the slot exit (y = 0 mm) and the location of the peak velocity (Fig. 11c→d): for instance, it is seen in Fig. 11 that at phase 90, ∂V /∂y is positive between y = 0 mm and y ≃ 4 mm. Therefore, in or close to the symmetry plane, the second term in the RHS of Eq. ( 4) is necessarily negative, and leads to a decrease of < v ′′2 >, since this term is actually 1 2 P 22 . This mechanism, due to a locally positive streamwise velocity gradient, similarly to the case of the flow in a convergent nozzle, can be expected in many unsteady flows, in particular in pulsed jets.

On the contrary, the first term, which is 1 2 P 11 , is positive. The production of k ′′ is negative because < v ′′2 > is larger than < u ′′2 > in the jet. This phenomenon was observed in the far field of a transient axisymmetric turbulent jet submitted to a large and sudden decrease in its ejection velocity [START_REF] Borée | Measurements and image analysis of the turbulent field in an axisymmetric jet subject to a sudden velocity decrease[END_REF]. During the velocity decrease, analysis of the transient turbulent field shows a reduction of the Reynolds-stress tensor anisotropy. It is thus clear that in this two different flows, the combination of anisotropy in the jet and decrease of the velocity at the slot exit after the maximum necessarily implies negative production. Kinetic energy is thus transfered from turbulence to the oscillating motion through the component < v ′′2 >. The other components < u ′′2 > and < w ′′2 > also lose their energy via redistribution to < v ′′2 >.

The acceleration of the fluid particle observed on Fig. 6 is thus due to the fact that the third term in the RHS of Eq. ( 3), the shear stress, is not sufficient to balance the first two terms, in which the pressure gradient is dominant.

With the standard k ′′ -ε ′′ model, which is unable to predict negative production, the picture is reversed: instead of providing energy to the oscillating motion, < v ′′2 > increases rapidly and this misprediction is sufficient to changes the sign of the momentum budget (Eq. 3). The previously quoted increase of the peak velocity between phases 45 and 90 is then missed: the peak velocity is much too weak, as seen in Fig. 6, which leads to too short a penetration of the jet into the quiescent air. Fig. 12 shows profiles of long-time averaged velocity V . Surprisingly, despite the strong overestimation of k ′′ , shown in Fig. 15, the k ′′ -ε ′′ model reproduces the spreading of the jet better than the RSM. This is a good example of the beneficial effect of the compensation of errors: basing the comparison of the relative performance of the models on long-time-averaged velocity profiles is clearly misleading.

The RSM also predicts regions of negative production at the periphery of the large eddies (Fig. 14a). A detailed analysis (not shown here) indicates that production in these regions are dominated by the first two terms in Eq. ( 4): < u ′′2 > -< v ′′2 > is negative and ∂V /∂y is positive.

Misalignment of the strain and anisotropy tensors

It is well known that negative production results from a misalignment of the proper axes of the strain and anisotropy tensors: in particular, in flows subjected to periodic irrotational strains, a time lag is observed between these two tensors, which is at the origin of a turbulence decay (see, e.g., [START_REF] Hadžić | Modeling the response of turbulence subjected to cyclic irrotational strain[END_REF]. The strain and anisotropy tensors are both 3 × 3 symmetrical tensors: the eigenvectors form orthogonal bases. Since the flow field is 2 dimensional, the third eigenvector is the same for both tensors and aligned with the z-axis. Let θ be the angle between the eigenvectors of the anisotropy tensor a ′′ and those of the strain tensor S in the x-y-plane, as shown in Fig. 16.

The anisotropy tensor a ′′ , defined by a ′′ ij =< u ′′ i u ′′ j > /k ′′ -2 3 δ ij , can be written:

a ′′ =        λ 1 0 0 0 λ 2 0 0 0 λ 3       
(5) in its own eigenaxes, where λ i 's are the three eigenvalues (with, arbitrarily, λ 1 > λ 2 > λ 3 ). In the same basis, the strain tensor reads:

S =        S 11 S 12 0 S 12 S 22 0 0 0 0        =        -β cos 2θ -2β cos θ sin θ 0 -2β cos θ sin θ β cos 2θ 0 0 0 0        ( 6
)
where β is the positive eigenvalue of S. Unlike eddy-viscosity models, the production term in RSM does not need modelling, and reads:

P k ′′ = -ka ′′ : S = kβ(λ 1 -λ 2 ) cos 2θ (7)
For instance, in the symmetry plane, just above the slot exit, at phase 60 (Fig. 14), the eigenvalues of the anisotropy tensor are λ 1 = a 22 and λ 2 = a 11 , and the eigenvectors are aligned with the y-and the x-axes, respectively. The positive eigenvalue of S is β = ∂V /∂y, and its eigenvectors are aligned with x-and the y-axes, respectively, such that θ = 90 deg, and the production is negative. Fig. 17 illustrates the correspondence between θ and the regions of negative production at phase 60.

The k ′′ -ε ′′ model is of course unable to capture such a mechanism because the anisotropy and strain tensors are forced to be directly in phase, leading to an overestimation of turbulent kinetic energy, as shown in Fig. 15, which is in turn at the origin of the wrong prediction of the peak velocities, as shown in section 4.2. Therefore, the stress-strain lag, whose role was emphasized by [START_REF] Hadžić | Modeling the response of turbulence subjected to cyclic irrotational strain[END_REF] in the case of homogeneous turbulence subjected to periodic compression-dilatation strains, also plays a fundamental role in the present flow.

Vortex generation

Using standard vortex generation analysis, it can be shown that the vortex dipole created during a cycle is mainly driven by inviscid dynamics. Indeed, let Σ be the part of the x-y-plane delimited by a the line ∂Σ = ABCDE, as shown in Fig. 13. The time-derivative of the circulation Γ along ∂Σ is given by:

∂Γ ∂t = ∂Σ ∂ U ∂t • dℓ = Σ ∂ ω ∂t dS (8)
where the vorticity ω is the component in the z-direction of the curl of the velocity vector ∇ × U . In an inviscid 2D flow, the material derivative of vorticity is zero, such that:

∂Γ ∂t = - Σ ∂ ω U k ∂x k dS = - ∂Σ ω U k n k dℓ = - h/2 0 ω V dx (9)
where the third equation holds because BC is a wall, CD and DE are supposed sufficiently far from the blob of vorticity ( ω ≃ 0), and ω = 0 on the symmetry line EA. Just at the outlet of the slot, ∂U/∂y is small compared to ∂V /∂x, such that Eq. ( 9) reduces to:

∂Γ ∂t = - V 2 0 2 (10)
where V 0 is the velocity at the centre of the slot (x = 0).

The circulation along ∂Σ predicted by a purely inviscid dynamics can thus be estimated at every time using Eq. ( 10), from the experimental value of V 0 . Now, using the Lamb-Oseen [START_REF] Lamb | Hydrodynamics[END_REF] model for a vortex:

V θ (r) = Γ 2πr 1 -e -r 2 /R 2 (11)
an estimation of the velocity field is obtained (r and θ are the cylindrical coordinates with respect to the vortex centre). In Eq. ( 11), the only parameters that cannot be predicted by the inviscid theory are the location of the vortex centre and the radius of the vortex core R. In Fig. 18, the locations (x 0 , y 0 ) and (-x 0 , y 0 ) of the centres of the vortices were determined from the PIV velocity field, and the velocity profiles were extracted at y = y 0 . The radius of the vortex core R is then obtained from the velocity gradient ∂V /∂x at the centre of the vortices.

The velocity induced by the vortex dipole is the vectorial sum of the velocities induced by each vortex:

V tot θ = V θ (r) + V θ (r * ) ( 12
)
where r * 2 = (x + x 0 ) 2 + (yy 0 ) 2 . Taking the x-derivative of Eq. ( 12) at (x = x 0 , y = y 0 ), shows that R can be evaluated by:

1 R 2 = 2π Γ    ∂ V tot θ ∂x x 0 ,y 0 + Γ 8πx 2 0    (13)
The excellent agreement in Fig. 18a of the Lamb-Oseen model with the experimental profile shows that, up to this phase, the vortices are mainly driven by an inviscid dynamics: despite a significant level of velocity fluctuations (above 5 %), the flow is more transitional than fully turbulent when it issues from the slot. It is moreover difficult to estimate the relative contribution in the velocity fluctuations of the deviation from periodicity in the diaphragm motion. In Fig. 18b, it can be seen that the vorticity injected through the slot has partly been dissipated at phase 135. Since the radius of the vortex cores is not given by the inviscid theory but rather chosen to fit the profile, the Lamb-Oseen vortices still reproduce accurately the experiment in the region -0.002 m < x < 0.002 m. However, outside this region, it is observed that the inviscid theory predicts too strong a backflow, which shows that the circulation Γ is overestimated, since in this region, Eq. ( 11) goes asymptotically to V θ (r) = Γ/(2πr). This conclusion is completly independant of the particular choice of R.

In such a transitional region, what is expected from turbulence models is to influence the flow as little as possible. As shown in section 4.2, the k ′′ -ε ′′ model strongly overestimates the turbulent energy. Therefore, a strong viscous term is artificially added to the vorticity equation, and, in turn, introduced in Eq. ( 9). The circulation is then dissipated, which explains at the same time why the acceleration of the fluid particle on the centreline is missed and why the convection velocity is underestimated (Fig. 18 shows that the convection velocity is a consequence of mutual induction of the vortex dipole).

Conclusions

The turbulence rate of the synthetic jet studied in the present work is above 5% at the slot exit and can reach 20% around the vortex core, which could question the use of turbulence models. Indeed, despite this significant level of fluctuations, the analysis of the dynamics of the vortex dipole shows that the flow is more transitional than fully turbulent. Nevertheless, in industrial applications [START_REF] Getin | Simulation numérique du contrôle actif par jet pulsés de l'écoulement turbulent autour d'un cylindre circulaire[END_REF][START_REF] Ben Chiekh | Synthetic jet control for flows in a diffuser: vectoring, spreading and mixing enhancement[END_REF][START_REF] Amitay | Modification of the aerodynamics characteristics of an unconventional airfoil using synthetic jet actuators[END_REF]) synthetic jets interact with the developed turbulence of the flow to be controlled: therefore, a turbulence model must be able to reproduce the dynamics of the synthetic jet region as well as the fully turbulent region. This is very challenging for statistical modelling.

The present work was dedicated to the investigation of the relative performance of two closure levels: eddy-viscosity and Reynolds-stress models. The solutions are free from turbulent eddies (i.e., represent only the oscillating motion), which shows that the standard models derived for statistically steady flows can also be suitable for flows treated in phase-averaging.

However, the results clearly demonstrate the superiority of the Reynolds-stress model. Indeed, the k ′′ -ε ′′ model severely alters the evolution of the vortex dipole: the generation of the large-scale periodic contra-rotating vortices are directly inherited from the inviscid dynamics, but their dissipation and their convection velocity are driven by the level of viscosity (molecular and turbulent). The k ′′ -ε ′′ model predicts a much too high turbulence energy, a rapid decrease of the peak velocity, and, consequently, too short a penetration of the jet into quiescent air. The Reynoldsstress model gives a much more realistic solution, but slightly underestimates the peak velocities, the convection velocity and the spreading rate of the jet.

The analysis of the turbulence production mechanisms shows that the discrepancies between the two models can be traced to the presence, during the deceleration phase, of a region of negative production that plays a crucial role in the global dynamics of the flow. Such a negative-production region is due to the combination of the anisotropy of turbulence in the jet and the positive gradient of velocity in the streamwise direction, which are both present in every pulsed jet. Since the k ′′ -ε ′′ model is not able to account for the stress-strain lag at the origin of the negative production, it predicts an erroneous, intense increase of turbulent energy, which is sufficient to change the sign of the momentum budget.

The present study clearly pleads against the use of linear eddy-viscosity models for this type of flows. [START_REF] Rumsey | Summary of 2004 CFD validation workshop on synthetic jets and turbulent separation control[END_REF] arrived at the opposite conclusion by synthetizing the results of the 2004 CFD validation workshop on synthetic jets and turbulent separation control : they concluded that linear eddy-viscosity models produce the best results among URANS models, whatever the formulation used (SA, standard k-ε, SST). Indeed comparisons of long-time-averaged quantities are misleading, because of compensation of errors. When phase-averaged quantities are compared in details, the deficiencies of linear eddy-viscosity models clearly appear. In order to evaluate the ability of the models to reproduce the unsteady dynamics of the flow, it is much more relevant to focus on the behaviour of the vortex dipole, in terms of location, intensity, convection velocity and penetration into the ambient fluid: the Reynolds Stress Model is able to reproduce the correct dynamics, contrary to the k ′′ -ε ′′ model. However, since the turbulence is far from equilibrium in the region close to the slot exit, it is also difficult for a Reynolds-stress model to reproduce accurately the interaction between the turbulence and the vortex dipole. [START_REF] Bremhorst | Two-time-scale turbulence modelling of a fully-pulsed axisymmetric air jet[END_REF] showed the beneficial effect of multiscale k ′′ -ε ′′ modelling in an axisymmetric pulsed jet. Their study, together with the present study, suggests to investigate the performance of multiscale Reynolds-stress models (e.g., [START_REF] Schiestel | Multiple-time-scale modeling of turbulent flows in one-point closures[END_REF][START_REF] Cadiou | A two-scale second-moment turbulence closure based on weighted spectrum integration[END_REF].

Fig. 2 .Fig. 4 .Fig. 5 .Fig. 6 .

 2456 Fig. 1. Geometry, domain size and coordinate system.

  Fig. 7. Contours of phase-averaged velocity V at phase 90 (blowing). Dashed lines = negative contours.

Fig. 10 .Fig. 11 .

 1011 Fig. 8. Contours of phase-averaged velocity V at phase 225 (suction). Dashed lines = negative contours.

Fig. 16 .Fig. 17 .Fig. 18 .

 161718 Fig. 16. Definition of the angle θ between the eigenvectors of a ′′ ( → e 1 , → e 2 ) and those of S

  ′′ -ε ′′ Fig.14. Contours of P k ′′ and velocity vectors at phase=60.

	0.008					0.008
	006 0.007 1 8 y (mm) 2 3 4 5 6 7	m2/s2 1.5E+01 1.3E+01 1.0E+01 9.2E+00 7.3E+00 6.3E+00 6.0E+00 5.3E+00 5.0E+00 4.4E+00 4.0E+00 3.4E+00 3.0E+00 2.4E+00 2.0E+00 1.5E+00 1.0E+00 5.0E-01	y (m)	0.001 0.008 0.007 0.006 0.005 0.004 0.003 0.002	m2/s3 9.3E+04 8.6E+04 8.0E+04 7.3E+04 6.6E+04 5.9E+04 5.3E+04 4.6E+04 3.9E+04 3.3E+04 2.6E+04 1.9E+04 1.2E+04 5.7E+03 -5.0E+02 -5.0E+03 -1.0E+03	0.007 0.002 y(m) 0.003 0.004 0.005 0.006	m2/s2 1.5E+01 1.3E+01 1.0E+01 9.2E+00 7.3E+00 6.3E+00 6.0E+00 5.3E+00 5.0E+00 4.4E+00 4.0E+00 3.4E+00 3.0E+00 2.4E+00 2.0E+00 1.5E+00 1.0E+00 5.0E-01	y (m)	0.001 0.008 0.007 0.006 0.005 0.004 0.003 0.002	m2/s3 9.3E+04 8.6E+04 -5.0E+03 -1.0E+03 8.0E+04 7.3E+04 6.6E+04 5.9E+04 5.3E+04 4.6E+04 3.9E+04 3.3E+04 2.6E+04 1.9E+04 1.2E+04 5.7E+03 -5.0E+02 m2/s2 1.5E+01 1.3E+01 1.0E+01 9.2E+00 7.3E+00 6.3E+00 6.0E+00 5.3E+00 5.0E+00 4.4E+00 4.0E+00 3.4E+00 3.0E+00 2.4E+00 2.0E+00 1.5E+00 1.0E+00 5.0E-01
	-2	(a) RSM (a) PIV 0 2 (b) k x (mm) x(m) -0.003 -0.002 -0.001 0 0.001 0.002 0.003 0.001 x (m) -0.002 -0.001 0 -0.002 -0.001 0 0.001 0.002 0.001 0.002 (b) k ′′ -ε ′′ x (m) (c) RSM Fig. 15. Contours of k ′′ at phase=60.

Note that the symmetry of the flow is not imposed in this study, in order to enable the appearance of non-symmetrical solutions. However, only symmetrical solutions have been obtained.

Acknowledgements

The authors gratefully acknowledge Jacques Borée for many useful discussions during the course of the present research project.


Jean-Paul Chemla

Bernard Riera

HOME AUTOMATION WITH A PLC: FROM RE-QUIREMENT TO COMMISSIONNING

Keywords: PLC, Home Automation, Virtual home 1

This paper deals with the programming part of home automation. We will present a solution based on central PLC which controls all the basic devices in a home. The PLC is connected to buttons and sensors and can activate lights and shutters using relays. We will focus on the PLC program. We propose a methodology based on the definition and the use of a library of function blocks for each common device like buttons, lights and rollers. Once the library is defined, most of the programming task will consist in defining instances of the blocks, connecting them to define rules. The library should also take into account that there will be some human machine interface (HMI) to interact with the controlled home. Some solutions to build this kind of HMI are presented in the paper. In order to validate the controller, we propose to use virtual commissioning with a virtual house. In our case, we have used the software HOME I/O from Real Games, to perform this stage. At least, a real home with these solutions has been realized and comments about the use of this controlled home will be provided.

INTRODUCTION

Home automation

Home automation is a growing trend that is becoming more and more popular and affordable every day [START_REF] Dc Khedekar | Home Automation-A Fast -Expanding Market[END_REF]. Today this term is rather replaced by smart house. It is used in electronics technology, automation, IT and telecommunications to ensure safety functions (alarms, cameras and remote monitoring), living comfort (automation and job scheduling daily, remote control, home care ...) and energy management (optimization of heating, lighting ...). The operating principle of a home automation system is to centralize the control and monitoring. Two different trends are on this business: wired solutions often expensive but offer a complete solution for a house, and home automation boxes which uses the power of Internet and wireless. With or without subscription, they allow open use and can be controlled from a computer, a smartphone or a touch pad. Installation is very simple, and takes just a few minutes by a nonexpert user but it manages only a limited part of the house.

We consider here a global solution with a central PLC connected (often wired) to sensors, buttons and devices like lights, shutters, buzzers (for alarm) [START_REF] Figueiredo | Predictive control applied to builing automation : an experimental setup[END_REF].

PLC and Home automation

Unlike smart devices, this solution can provide the possibility to interact with both classical buttons in the home and smartphones. Moreover, the PLC with his relays will be able to switch all kind of plain devices (lights, roller shutters) without buying specific ones. Even if the behavior is classical (just ON or OFF for a light), there will be no compatibility issues. The proposed solution is intended to be installed at the same time that the electrical wiring of the home because it is intended to control all home devices. As it is a non-computer based, industrial device, the PLC is more secured [START_REF] Jacobson | A risk analysis of a smart home automation system[END_REF] and reliable than home automation boxes.

Context

The project will focus on basic devices of a home: lights, shutters and alarm. The main role of the PLC is to replace the classical wiring of the home to make it "smart". Examples of behaviour will be provided in the programming part of this paper. The PLC has on his inputs:

• Buttons or switches (usually used for lights and shutters)

• Motion sensors in corridors • Proximity sensor to detect if windows and doors are closed

The outputs of the PLC are connected to relays to activate:

• Lights • Shutters • Buzzers
It is possible to make a parallel with modern cars: the buttons and sensors are connected to a central controller which can switch on or off the headlights, lock or unlock the doors, go up or down the windows etc.

THE LIBRARY

As specified in the IEC 61131.3 norm, PLC can be programmed using programs, functions and function blocks. The environment provides built-in functions (as AND, XOR) and function blocks (as RTRIG, SR).

The user can define his own programs, functions and function blocks in multiple languages (IL, LD, ST, SFC, FBD).

A large home can have lots of each device (buttons, lights…) and programming can be a repetitive task and difficult to maintain. Moreover, as we want a smart behaviour, each device can lead to a complex program. That's why we will define our own function blocks for buttons, shutters, switches, lights, buzzers. These blocks will have input, output, input and output and private variables. They also will have a program part. Most of the programming software for PLC have the possibility to define new blocks. Let's define some of these function blocks.

• Button:

-Input: the button itself (an input of the PLC) -Outputs: short press, long press, double press

The button is simple: the function block detects if the button has been pressed shortly, pressed a long time (more than 1/2s) or doubled pressed. Most of our smartphones have a button with different functions according to these types of press.

• Light:

-Inputs: switch (bool), off (bool), timer (bool), time for timer (time) -Outputs: the light itself (an output of the PLC) -Input-Outputs: HMI_switch, HMI_timer

We thought that a light can be switched, like traditional ones: from ON to OFF or vice versa. Sometimes, it could be convenient to turn the light ON for a limited time. Then it will turn OFF after that the time is up.

In that case, we prevent 1 minute before time's up by turn the light OFF for 5 seconds. If the timer input becomes true during the timer is running, the light will shut off except if it is in the last minute (in that case, the timer will be relaunched). For the HMI, it could be nice to have only 1 variable which tells if the light is ON or OFF and to control the light. That's the role of the input-output variables: they are both information and control. HMI_switch is to switch the light, HMI_timer is for timer. For example, if the light is ON for a limited time, HMI_switch is OFF and HMI_timer is ON. We define that if Up or Down inputs are ON shortly, the shutter moves until his end (up or down). If these inputs are ON for a long time (long press on the buttons) the shutter will stop when these inputs return to OFF (button is released). As there is no sensor of position (neither if they are opened or closed) in most of shutters, we have to define the time taken to completely open or close the shutter. Inside the block, the position of the shutter is calculated. This position is an output of the block as it can be useful for the HMI. Another value is useful for HMI: the set point. This variable is a percentage of closure. It is of course an input as the user can change it through his HMI. It is also an output of the block because we want that this value is set to the actual position when the shutter is not moving (this value is both an information and a control).

Other blocks have been defined for buzzers and switches for example.

Programming the block is the most technical part of the PLC programming because we have to:

• Anticipate a lot of situations that can happen. For example: the shutter is moving down, what happens if the down button is pressed? (in this case, we decided that the shutter has to stop). Other example: to stop the shutter, the user can press both buttons Up and Down. But the two inputs won't be synchronous. • Be careful with all input and output variables because they can create instabilities. For example, if the user press a button to switch ON a light, the program sets ON the light and before that the HMI_switch is set to ON, the program can detect that the HMI_switch is in a different state from the light and light can return to OFF.

PROGRAMMING AND HMI

Once the library is defined, it remains:

• To declare all the inputs and output variables • To create instances of the blocks for each button, light, shutter, etc… • To create the variables for the communication with the HMI • To create and manage variables for global actions • To create an SFC to manage the alarm mode • Then to program the behaviour of the house mainly by connecting the blocks between them and to other variables (global actions, HMI, alarm mode). This last part is not difficult and similar to "If that then that" rules programming but using blocks.

Programming

We give some details about some of the above points.

The alarm mode

We defined 5 different states of the alarm: "Presence mode", "Verifications before leaving", "Absent mode", "Waiting for digital code", "Intrusion". The SFC below gives the conditions of the transitions between states. The leaving condition can be for example a double press on a button near the door. When in the state "Verification before leaving", buzzer can beep quietly if every window and door are closed but long beeps can occur if it is not the case.

Figure 1.

SFC of the alarm mode

Examples of program using blocks

In the figure 2 below, there is the program equivalent to "If the button near the door is pressed for a short time OR if the door is opened AND there is no daylight AND the outdoor light is OFF THEN the outdoor light is ON for 15 minutes". The global variable "global_outdoor_light_off" is calculated in an another part of the program. It can be set ON by multiple ways : the global light off of all the home or by a button on a HMI.

Human Machine Interfaces

As mentioned above, the human machine interface has been prepared by adding input/output variables in the definition of some blocks. In the program, variables have been created and localized to be accessible via the network with the Modbus TCP protocol. We have chosen this protocol because the PLC we used can manage it and as it uses TCP/IP, we can reach the PLC through different networks and 3G/4G. There is no need of an extra computer to have the interface (unlike [START_REF] Bingol | Kubilay Tasdelen Web-based Smart Home Automation: PLC-controlled Implementation[END_REF]). The interface itself has to be easy to use and has to give a very quick access to controls. A smartphone has these qualities.

In the Apple App store or in the Google Play store, we can find applications which can connect to PLC using industrial network protocols and we can create interfaces to interact with the PLC. For example: mySCADA (both stores) can use multiple protocols. Other apps are specific to a brand of PLC. We have chosen HMI Draw from SweetWilliams [START_REF]HMI Draw home[END_REF] in the Apple App store because it is ergonomic and has a large range of possibilities like using sensors of the device. In the figure 3 below, there is a view we can obtain with HMI Draw. In this view, we can see a switch for a light.

VIRTUAL COMMISSIONNING

In order to test the program of the PLC, we looked for a solution with a virtual home which can be connected to a real PLC. This solution is the software called HOME I/O [START_REF]Home I[END_REF]. This software is much more a serious game than a simulation. It is the result of a 3-year R&D project "DOMUS" (2011-2014) partially founded by the French Ministry of National Education, between CReSTIC lab from the University of Reims Champagne-Ardenne and Real Games, a Portuguese company [START_REF] Riera | a virtual house for control and STEM education from middle schools to Universities[END_REF]. The idea has been to bring a virtual house into the class room, usable from middle schools to universities, adapted to learners (Y generation) and teachers and suitable for control and STEM (Control, Science, Technology, Engineering and Math) education. HOME I/O has been designed by applying a systemic approach of education. For that, HOME I/O enables to study the house from different points of view (automation, energy efficiency, smart home …), as a whole or as a set of subsystems, but taking into account the environment. The possibilities of pedagogical scenarios (i.e. situations, questions, problems) offered by HOME I/O depend only on teachers' imagination. HOME I/O seems to be a very innovative pedagogical tool and it was awarded with the prize "Le Coup de Coeur" in the 6th Edition of "TROPHÉES DES TECHNOLOGIES ÉDUCATIVES 2014" at EDUCATEC-EDUCATICE in Paris.

Hence, HOME I/O is real time simulation software (Fig. 5) of a smart house and its surrounding environment, designed to cover a wide range of curriculum targets within Control, Science, Technology, Engineering and Math, from middle schools to universities. This virtual house, like a video game in first person (FPS, first person shooter), becomes a place of discovery and experimentation. HOME I/O is a realistic modern virtual house. It is necessary to make choices in equipment of the house and the possibilities of automation. It is possible to control lighting, rollers, electric heating, an alarm system for intrusion safety and home safety (fire). Altogether HOME I/O offers 174 interactive objects (lights, switches, inverters, roller blinds, gate, garage door, sirens, motion detectors, lighting, fire detectors, door opening sensor, remote control, radiators, thermostats) that can be used to automate the home or as a source of logical or numerical data (Fig. 6).

Figure 6

Interactive objects in HOME I/O activities, extremely immersive and addictive are a natural and understandable worry for parents and educators, computer games can teach young children to read and count, as well as help middle school students to learn about science and technology. It is why serious games are now considered as a possible educational modality. The vocation of a serious game is to make learning attractive using the techniques of the game (interactivity, challenge ...). The rise of serious games is the consequence of the development of the video game industry, whose sales worldwide now exceeds film and music industry.

Fig. 1. HOME I/O head-up display From the beginning of the project, the idea was to design a serious game where teachers and learners can play a role. Serious games for schools ask the essential questions of the role assigned to the teacher and how to effectively integrate them into an educational classroom sequence. That is why the real issue is not the serious game and its definition but the development of digital tools for students, effective in terms of learning, and integrating into innovative teaching practices where teacher has its place. Design a digital tool promoting investigative approach and allowing the teacher to imagine and implement multiple problem situations is the vision that led to the development of the software HOME I/O. The use of technology and video games technologies (rendered graphics and sound, interactivity, attractiveness) is seen as a way to promote "situation awareness," student motivation and thus foster their learning. In our case, serious games are considered primarily as a digital tool that the teacher will be and entertaining experience where the teacher remains in control of the lesson plan (i.e. the game design).

The design stage

The design of HOME I/O was performed by placing teacher and student in the heart of the approach. The initial specifications stated that the software should be developed to be used in France as part of the middle school technology curriculum, especially around the scope "comfort and home automation" and be complementary to Web resources and material available today. The virtual house should also be adapted to high school and universities. In middle school, it is important for the software to be autonomous and does not have to require other applications to be used. However, for high schools (but also higher education), the house could rather be seen as a source of data that can interface with third party devices and software (LabView, Matlab ...) or hardware (microcontrollers, programmable logic controllers, Arduino, Raspberry Pi ...). The idea was therefore to provide an educational tool contributing to the continuity of STEM education, from middle schools to universities. We applied a systemic approach (Riera and Debernard, 2003) of control and STEM education, it is to say: concentrate on interactions between elements, consider the effects of interactions, consider the global perception, deal to a multidisciplinary education, deal to an action by objectives, and knowledge of goals. For that, HOME I/O has been designed in order to be able to study the house from different points of view (automation, energy efficiency, smart home), as a whole or as a set of subsystems. Complexity can thus be adapted to the different levels and learning objectives. The design approach was also based on a thorough study of itions and day). e and lation mena, t the ). The es are ances lighting, fire detectors, door opening sensor, remote control, radiators , thermostats) that can be used to automate the home or as a source of logical or numerical data (Fig. 6). For example, it is possible to know and monitor power consumption.

Fig. 3. Earth map

Similarly, it is possible to change the weather conditions (wind speed and direction, cloud cover, humidity and minimum and maximum temperatures for one day). Environmental parameters play on the temperature and brightness inside the house (Fig. 4). A simplified model has been created for real-time simulation of heat transfer. This model includes radiation phenomena, convection and conduction, and takes into account the physical properties of building materials (no modifiable). The exchanges between air masses of different temperatures are also simulated. The model is influenced by the disturbances created by the opening of doors and windows. The temperature of the outside air acts on the temperature of the house (heat transfer and airflow). The wind facilitates the transfer of heat between the house and the outside air. The walls of the house located upwind are more reactive to heat transfer with the outside air. Cloud cover decreases the effects of radiation from the sun and sky. Finally, the moisture changes the dew point, which acts on how the outside air affects the temperature of the house. The goal set with HOME I/O is educational. This means that the house is a "black box", a data source and a place for experiments. Dynamic models in HOME I/O are purely descriptive. In a house, the physical phenomena have very different dynamics. For example, a garage door or a roller blind can be closed in seconds. However, to achieve the steady state lighting, fire detectors, door opening sensor, remote control, radiators , thermostats) that can be used to automate the home or as a source of logical or numerical data (Fig. 6). For example, it is possible to know and monitor power consumption.

Fig. 6. Interactive objects in HOME I/O

To be adapted to Y generation, interaction of the user and the virtual world obeys the uses of FPS games, and thus guarantees a rapid appropriation of this digital tool by students. Moving into the house is done using either the keyboard or mouse, or a joystick. HOME I/O is also compatible with the Oculus Rift virtual reality helmet. The functional realism is just as important as graphics rendering. We therefore based on current automation tools to avoid creating false problems. HOME I/O has been designed to highlight the evolution of home automation. Software solutions adopted by the home automation box manufacturers have therefore been studied from a functional point of view With another software, CONNECT IO (from the same company) this virtual home can be controlled by an external device like a PLC, using communication protocol (like Modbus TCP) or using a DAQ card with wired connections. This virtual commissioning gave us the opportunity to test function blocks, program and HMI interfaces.

REAL HOME

One of the author has a PLC in his home and applied the methodology given in this paper. This experience leads to these remarks:

-There is two kind of wires in the house: for control, with thin wires and for power with classical copper wires. Both of them are connected to the PLC. We could have chosen a network for the control wiring. For example, KNX is often used in that cases and lots of PLC can manage this network.

-If someone has a classical behaviour in the house, the actuators will also have a classical behaviour. Except the proximity sensors in the corridors and the timers for some lights, it seems to be a classical house. -Global control for shutters (with a long press on a shutter button) is very convenient -It is possible to switch on multiple lights with 1 button.

-We can add easily new features when we think about it by modifying the program section without creating new variables and without changing the function blocks. For example, turning on a light behind a door when opened like in figure 2 was done by adding the OR and AND blocks. -The alarms part of the program is mainly used to switch off all lights and to verify that all windows are closed. The PLC also sends email when someone remove the alarm which is convenient to know when children come back from school. -The PLC is running since more than 5 years without major trouble. Minor ones have been experienced at the beginning. We made some adjustments in function blocks and since years, we didn't change our library even if the program itself has changed many times. -We changed the app for the HMI independently to the program.

-Even if it is not very difficult to change rules, only people who are able to program a PLC are able to modify them.

CONCLUSION

We presented a methodology to program a central PLC which controls all standard devices in a home. The library we have developed is based on standard languages and can be implemented in quite every brand of PLC. We use a Scada application to provide a convenient interface. The original point of our paper is the commissioning on both a simulated and a real home. The behaviour of the programmed home offers functionalities like global controls, automated lights ON or OFF, interface with a smartphone, emails on events etc... To improve our work, we will focus on the possibility to add connected objects to our central PLC. One way to make this connection should be through an home automation boxes which can communicate with the PLC.

Figure 2

 2 Figure 2 Example of Function block diagram program for a light

Figure 3

 3 Figure 3Example of a view with the App HMI Draw

Figure 4 . 1 Figure 4 . 2 Figure 4

 41424 Figure 4.1 Figure 4.2 Figure 4.3 PLC configuration Tag configuration Switch configuration

Figure 5 HOME

 5 Figure 5 HOME I/O head-up display (left) and view of the kid room (right)

Fig. 2 .

 2 Fig. 2. HOME I/O kid room


Prof. V Yue Tang

Idir Benaissa

Dr Mathieu Huynh

Dr Laure Vendier

Dr Noël Lugan

Dr Stéphanie Bastin
email: stephanie.bastin@lcc-toulouse.fr

Philippe Belmont

Vincent César
email: vincent.cesar@lcc-toulouse.fr

Véronique Michelet
email: veronique.michelet@unice.fr

Guy Lavigne

Michelet

Prof P Belmont

An original L-shape, tunable N-Heterocyclic Carbene platform for efficient gold(I) catalysis

Keywords: N-heterocyclic carbene, gold, organometallic complexes, catalysis, atom-economical reactions

The synthesis and characterization of original NHC ligands based on an imidazo[1,5-a]pyridin-3-ylidene (IPy) scaffold functionalized with a flanking barbituric heterocycle is described as well as their use as tunable ligands for efficient gold-catalyzed C-N, C-O and C-C bonds formations. High activity, regio-, chemoand stereoselectivities are obtained for hydroelementation and domino processes, underlining the excellent performance (TONs and TOFs) of these IPy-based ligands in gold catalysis. The goldcatalyzed domino reactions of 1,6-enynes give rise to functionalized heterocycles in excellent isolated yields under mild conditions. The efficiency of the NHC gold 5 Me complex is remarkable and mostly arises from a combination of steric protection and stabilization of the cationic Au(I) active species by ligand 1 Me .

Over the last two decades, homogeneous gold catalysis has been recognized as a game changer in modern organic synthesis, due to the opportunity to access high molecular complexity from relatively simple substrates in an atom-economical step. [1] Monodentate ligands L play a major role in tuning the activity, stability, and selectivity of the gold catalysts. [1,2] Although several families of privileged ligands have emerged such as Buchwaldtype dialkyl(o-diaryl)phosphines [3] or N-Heterocyclic Carbenes (NHCs), [4] the quest for new and broad spectrum ligand systems in gold catalysis is still actively pursued. A rational and efficient strategy was recently shown to be the functionalization of Buchwald-type phosphines, either by grafting a basic group to direct the anti-nucleophilic attack of alkynes, [5] or by introducing a cationic charge to enhance the electrophilicity of the gold center. [6] Transposing this approach to NHC chemistry by derivatizing the standard 1,3-bis(2,6-diisopropylphenyl)imidazol-2-ylidene (IPr) ligand into bulkier and/or functionalized NHCs, met only limited success in gold catalysis, since nitrogen substituents remain too far from the gold center to exert an efficient steric environment. [START_REF] Hueber | For successful implementations of this approach[END_REF] The imidazo[1,5-a]pyridin-3-ylidene platform (IPy), first disclosed independently by the groups of Lassaletta and Glorius in 2005, [8] places the R group at the C5 position in close proximity to the metal center as a consequence of the annelation of pyridinyl and imidazolyl rings (Figure 1). It may thus be regarded as the geometrically carbene analogue of Buchwald-type phosphines. [9] Interestingly, IPy-based ligands and related N-fused heterobicyclic carbenes were already coordinated to gold(I) and led to efficient catalytic species in some instances (complexes A-C, Figure 1). [10] In line with our dual interests in the design of functionalized NHCs [11,12] and the development of new selective gold catalytic systems, [13] we report herein a novel class of functionalized IPy ligands for gold(I) catalysis, whose C5 position is substituted by an anionic (ligand 1 -) or neutral barbituric heterocycle (ligand 1 E ), as well as promising preliminary results in gold-catalyzed hydroelementation and domino processes. The barbituric heterocycle was chosen as it combines urea and malonate units into its structure, which actually may serve for further derivatization or during catalysis. [14] The zwitterionic airand water-stable precursor 1•H was prepared on gram-scale by reacting [2•H]Br with deprotonated hexahydropyrimidinetrione 3 through a facile aromatic nucleophilic substitution (SNAr) (Scheme 1). [12] Compound 1•H was fully characterized by spectroscopic and analytical techniques and its molecular structure was firmly established by single crystal X-Ray diffraction (XRD) (Figure 2). [15] Scheme 1. Synthesis of anionic gold(I) complex [4](PPN) and neutral gold(I) complexes 5 Me and 5 Tf . Dipp: 2,6-diisopropylphenyl.

The anionic gold(I) complex [AuCl(1)]Li ([4]Li) was then cleanly obtained by reacting the stable free NHC [1]Li with AuCl(tht) at low temperature. [16] For solubility reasons, a salt metathesis was carried out by adding bis(triphenylphosphine)iminium chloride (PPN)Cl to the crude solution of [4]Li in CH2Cl2 to afford the complex [4](PPN) in 66% yield after purification. The molecular structure of [4](PPN) was confirmed by an XRD experiment (Figure 2). In order to tune the overall electronic impact of the NHC ligand onto gold center, [10f,17] complex [4]Li was then reacted with methyl triflate or triflic anhydride to furnish the stable, neutral Au(I) complexes 5 Me , and 5 Tf , respectively, in which one oxygen atom of the malonate group is substituted. This Ofunctionalization of the malonate unit is accompanied by the loss of the symmetry plane present in the anionic complexes [4](Q), as illustrated by the splitting of the signals in the 1 H and 13 C NMR spectra of 5 Me and 5 Tf . Additionally, complex 5 Me was characterized by XRD (Figure 2). An in-depth analysis of the topographical steric maps of the ligands in the crystal structures using Cavallo's SambVCa software (Figure 2), [18,16] reveals that both ligands display high encumbrance over the gold center and even in proximity to the potential binding site in trans position to the carbene center (red and orange zones). Such a confinement of gold catalyst was previously shown beneficial for high activity and selectivity. [START_REF] Hueber | For successful implementations of this approach[END_REF] Additionally, the O-methylation is accompanied by a significant increase of the overall steric pressure of the NHC ligand from %Vbur = 47.8 for 1 -in [4](PPN) to %Vbur = 53.4 for 1 Me in 5 Me . Having the three complexes [4](PPN), 5 Me and 5 Tf in hand, we then evaluated their respective efficiencies in gold(I)-catalyzed hydroelementation and cycloisomerization/domino processes. As a first benchmark reaction, we chose the intermolecular hydroamination of terminal alkynes with aniline derivatives, [19] for which several NHC-Au complexes have shown to be active, [20] albeit with much less efficiency compared to phosphine-based gold(I) complexes. [21] Gratifyingly, under neat conditions, the reaction between phenylacetylene and aniline smoothly proceeded to completion at low temperature (30°C) using 0.5 mol% of catalysts [4](PPN), 5 Me or 5 Tf and 1 mol% of NaBAr F 4 (Ar F = 3,5-(CF3)2Ph) as chloride scavenger (Scheme 2, Eq. 1). At a much lower catalyst loading (0.01 mol%), pre-catalyst 5 Me appeared to be more active than the anionic [4](PPN) and 5 Tf , which arises from a fine-tuning of the electronics of the lateral malonate unit for efficient stabilization. [START_REF]The stoichiometric reactions of 5 Tf and 5 Me with an excess of aniline led to no reaction in the case of 5 Me and a displacement of the Tf group in 5 Tf to give PhNHTf and the zwitterionic complex[END_REF] Interestingly, under these conditions, complex 5 Me afforded the imine 6a in a 76% yield, which corresponds to a catalyst Turn Over Number (TON) of 7600, the highest TON reported so far for a NHC-Au catalyst in hydroamination reaction. Complex 5 Me was also highly active in the hydroamination reaction of a more reluctant substrate such as 1-dodecyne with aniline at 0.01% cat. loading and 80°C showing a TON of 3500 after 23h. This result compares well with L. Zhang's report using an amide-decorated Buchwald-type phosphine system (TON: 3900) under similar conditions. [5c] These promising results prompted us to investigate the challenging intermolecular addition of secondary anilines onto phenylacetylene, which had been reported only once using 5 mol% of a CAAC-Au(I) system. [23] Using 0.5 mol% of 5 Me at 80°C, diphenylamine reacted with phenylacetylene to give exclusively the Markovnikov adduct 7 in 74% yield (Scheme 2, Eq. 2).

Scheme 2 Gold-catalyzed intermolecular addition of amine and carboxylic acids to phenylacetylene and 1-dodecyne.

The IPy-Au complexes were then evaluated in the challenging intermolecular addition of carboxylic acids to alkynes (Scheme 2, Eq. 3). [5c,24] Gratifyingly, while complexes [4](PPN) and 5 Tf were inactive, complex 5 Me efficiently catalyzed the addition of benzoic acid and 4-methoxybenzoic acid to 1-dodecyne to produce the Markovnikov products 8a-b in good yields, 88%, and 74%, respectively, with a catalyst loading of only 1 mol%. This is remarkable considering that IPrAuNTf2 was previously reported inactive under these conditions. [5c] The third gold-catalyzed reaction was the atom-economical domino process implying 1,6-enynes and a nucleophile, for which carbon-and oxygen-functionalized nucleophiles are compatible. [25] These transformations have been described in the presence of phosphine-based ligands as well as NHC-based ligands. [25,7b] The three complexes [4](PPN), 5 Me and 5 Tf were firstly evaluated in the case of a carbon-electron rich aromatic indole derivative under standard conditions in diethyl ether at room temperature (Scheme 4, product 9a). Whereas moderate conversions were observed in the case of complexes [4](PPN) and 5 Tf , an excellent isolated yield of 95% was obtained in the case of gold complex 5 Me in only 45 minutes with 3 mol% catalyst loading. The lower reactivity of [4](PPN) may be explained by the poor electrophilicity and the lower steric protection of the postulated zwitterionic active species [Au + (1 -)] (see Figure 2), which experimentally evolves towards black nanoparticles, most presumably Au(0) nanoparticles. The lower stability and Lewis acidity, according to Gutmann-Beckett method, [16] of the gold complexes 5 Tf and 5 Me may account for the significant difference of activity between both derivatives (50 % vs. 95% yield). The high activity of complex 5 Me prompted us to further study its properties in other domino processes and also with a comparison with other NHC-based gold complexes such as IMesAuCl and IPrAuCl (Scheme 3, products 9a-9d). Scheme 3. Au-catalyzed domino cyclization/nucleophilic addition of 1,6-enynes with various nucleophiles. [a] in CH2Cl2; n.i. not isolated.

In all cases, the efficiency of complex 5 Me was outstanding compared to the other NHC-gold complexes in terms of yields, selectivity, or kinetic. The addition of electron-rich indole or 1,3,5trimethoxybenzene to N-tethered 1,6-enyne led efficiently to heterocycles 9b, and 9c in 71%, and 84% isolated yields, respectively. Low or moderate conversions were observed in the case of IMesAuCl or IPrAuCl. Moreover, the activity of complex 5 Me compared favorably with the activity of pseudo-half-sandwich NHCs bearing N-alkylfluorenyl arms. [7b] The addition of MeOH led to excellent yield for 5 Me and IMesAuCl, but a high difference of kinetic for the domino process (Scheme 3, product 9d) was observed in favor of 5 Me complex. The reactivity of carbontethered derivatives such as bis-sulfonated 1,6-enynes and diester was also studied and led to very good results as the functionalized alkenes 9e-9g were isolated in 83-89% yields. The large overall scope is noteworthy since the reaction conditions were amenable with carbon-as well as oxygen-based nucleophiles such as alcohols (allylic alcohol, MeOH, benzylic alcohol). The functionalized heterocycle 9c was obtained in high yield with TOF = 400 h -1 , 8 orders of magnitude higher than with IPrAuCl. Interestingly, in the case of MeOH as nucleophile, the higher efficiency of the gold 5 Me complex was again demonstrated as the observed TOF for the synthesis of 9d was 212 h -1 , compared to 35 h -1 in the presence of the standard IPr-based complex (see Supporting Information). With an efficient set of reaction conditions in hand, and considering the much faster conversions with gold complex 5 Me , we further challenged the efficiency of the catalyst by conducting the hydroarylation/ cyclization reaction on a 1-gram scale (3 mmol) and in the presence of 0.25 mol% of gold catalyst (see Supporting Information). Remarkably, the heterocycle 9c was isolated in 87% yield, which corresponds to a decrease in catalyst loading of more than 3 orders of magnitude from this work, and 10 orders of magnitude from previous studies. [25] In conclusion, we reported the straightforward and efficient synthesis of new NHC ligands based on an imidazo[1,5-a]pyridin-3-ylidene (IPy) scaffold functionalized by a flanking barbituric heterocycle and their use as tunable and efficient NHC ligands for original gold-catalyzed reactions. The efficiency of the related Au(I) complexes was evaluated in C-N, C-O and C-C bonds formation, according to hydroelementation reactions as well as domino processes, with high to excellent TONs and TOFs under optimized conditions. Further studies aiming at a deeper understanding of the key stereoelectronic features of this ligand class, [26] as well as on the possible transposition of these benefits to other types of transition-metal catalyzed reactions are underway in our laboratories.

Figure 1 .

 1 Figure 1. General depiction of the imidazo[1,5-a]pyridin-3-ylidene (IPy) platform, N-fused heterobicyclic carbene gold(I) complexes A-C , and ligand system (1 -, 1 E ). Dipp = 2,6-diisopropylphenyl.

Figure 2 .

 2 Figure 2. Top: Molecular structures of 1•H (left),[4](PPN] (center), and 5 Me (right) (ellipsoids drawn at 30% probability level). PPN cation and hydrogen atoms have been omitted for clarity. Bottom: Topographical steric maps of[4](PPN] (center) and 5 Me (right). Values in the four corners of the maps are the %Vbur of the NHC ligand in the corresponding quadrant.

Acknowledgements

This work was supported by the Agence Nationale de Recherche (ANR-16-CE07-0006), the Centre National de la Recherche Scientifique (CNRS) and the University Côte d'Azur. Y.T. and I.B. are grateful to ANR for PhD grants.


Reply to Comment on "Pressure enhancement in carbon nanopores: a major confinement effect", Phys. Chem. Chem. Phys., 2011, 13, 17163 Y. Long, aa J.C. Palmer, b B. Coasne, c K. Shi, a M. Śliwinska-Bartkowiak d and K.E. Gubbins a Van Dijk's Comment draws attention to the fact that the local pressure tensor, P(r), in an equilibrium system that is inhomogeneous is not uniquely defined at the nano-scale, as seen in Eq. (4) of our paper 1 , where the configurational term on the right involves an integral along a contour Cij connecting the centres of two molecules, i and j. The choice of this contour is arbitrary, resulting in the non-uniqueness of P. Physically, the nonuniqueness arises because the intermolecular forces are between two molecules located at ri and rj, i.e. they are non-local, whereas we wish to determine the pressure at a local point r, and there is no unique way to assign part of the non-local pair force to the local point r. To pursue the statistical mechanical approach to P one must first make an operational definition of the contour. All of this is well known in the statistical mechanical community. It is implicit in Kirkwood and Buff's 1949 paper 2 , and is discussed at some length in an Appendix of the seminal paper of Irving and Kirkwood (IK) 3 . Irving and Kirkwood 3 , writing 70 years ago, stated the situation clearly: "all definitions must have this in commonthat the stress 4 definition between a pair of molecules be concentrated near the line of centers. When averaging over a domain large compared with the range of the intermolecular forces, these differences are washed out, and the ambiguity remaining in the macroscopic stress tensor is of negligible order." Further discussion of these points were given later by Harasima (H) [START_REF] Harasima | [END_REF] , Schofield and Henderson 6 , Rowlinson 7 , Hafskjold and Ikeshoji 8 , Binder and coworkers 9,10 , Hatch and Debenedetti 11 and many others. One choice of contour that has been widely used is a straight line joining the centers of the two molecules. This was the contour used by Irving and Kirkwood (IK) in their 1950 paper 3 , and follows the line of force between the two molecules. Another quite widely used contour is that of Harasima 5 , in which the path runs from the position of particle i,   ,,

i i i
x y z , to position   ,,

j j i
x y z , and then to the position of particle j, x j , y j , z j ( ) . For a planar interface lying in the xy plane, there are two independent components to the pressure tensor, PN = Pzz normal to the interface, and PT = Pxx = Pyy acting parallel to the interface. It is the tangential pressure, PT, that is not uniquely defined.

In our paper we were at pains to point out this non-uniqueness of the pressure tensor, as we have in our previous papers on inhomogeneous systems 12,13 , and in subsequent papers quoted in the Comment. We say on page 2 "Eqn (4) is formally rigorous, but while the kinetic term is uniquely defined the configurational contribution is not, due to the arbitrary nature of the contour integral…." In a 2013 paper 14 (cited as ref. 4 of the Comment) we compare results for the tangential pressure, PT, for slit pores using both the IK and H definitions. As expected, the two definitions give different numerical results. However, they have some qualitative similarities. Next to the pore walls the positions of the peaks in PT coincide with the peaks in the density profile, but the heights and widths of these peaks differ. The integral of PT over z gives the same average of PT for both definitions, as it should. 6 In studies of adsorbates adsorbed on planar carbon surfaces or in slit-shaped carbon pores, provided the adsorbate wets the walls, as is the case in the paper 1 being discussed here, adsorption occurs in well-defined layers parallel to the wall (see Fig. 2

of ref. 1).

There is also evidence (described briefly below) of strong compression of the adsorbate molecules close to the walls. The area where there is disagreement, is how best to quantify this compression. One approach might be to try to define a microscopic pressure in the pore phase that results from the intermolecular forces at play, and that would give at least an approximate measure of the compression of the molecules. An alternate approach would be to discuss the compression in terms of the local density of the molecules near the wall.

Van Dijk proposes to use the pressure of the bulk phase in equilibrium with the pore phase, Pbulk. This, together with the temperature, determines the chemical potential in both the bulk and pore phases. The configurational part of the bulk pressure provides information about the average intermolecular forces acting on the molecules in the bulk phase. However, it provides no information about the local properties of the pore phase, or the intermolecular forces at play there, or how they change with the nature of the adsorbate, adsorbent, pore width and shape, etc., and so does not address the question of interest here. It is these latter intermolecular forces that determine the properties and behaviour of the pore phase. definitions, and two alternative contours), for LJ argon in a slit carbon pore of width H = 3σ at 87.3K, 1 bar bulk pressure. 15 Because the local tangential pressure is not uniquely defined, we must exercise caution in attaching physical meaning to the resulting values. However, as shown by Harasima 5 , and later by Schofield and Henderson 6 , the integral of PT(z) over z is unique, so that (for example)

P T ,IK (z) dz = P T ,H (z) ò ò dz (1)
where IK and H refer to the Irving-Kirkwood and Harasima definitions, respectively. For wetting systems, such as those in ref. 1, we find that it may be sufficient to integrate over a given peak in the local tangential pressure to obtain a well-defined effective tangential pressure for a given adsorbed layer. 16 An example is shown in Figure 1 for the same intermolecular forces as used in ref.

1. The effective tangential pressure shown by the horizontal line results from integration of PT(z) over the width of the adsorbed contact layer, t c , using four different definitions of the contour Cij; two of these definitions are the IK and H contours, and the other two were chosen randomly. When spatially averaged over the width of the adsorbed layer, these 4 definitions of the contour give pressures that agree within the accuracy of the simulation results. Thus, the spatial average tangential pressure shown in Fig. 1 is given by:

P T ,avg = 1 t c P T z min H /2 ò z ( ) dz (2)
where H = 3σ and

t c = H / 2 -z min .
The size of the dead space near the wall, z min , is the minimum distance of a molecular centre from the surface, z = 0, taken to be the value of z where the density profile departs from zero. The tangential pressure averaged in this way is approximately 0.36 GPa. If, instead of averaging over the width of a single peak in PT, we average PT over the entire pore width from z =0 to z = H, without accounting for the dead space, we again find the 4 definitions give the same averaged tangential pressure, with a value of approximately 0.145 GPa. While these averaged tangential pressures are substantially lower than the peak values of the local pressures, they still show an enhancement over the bulk phase pressure (10 -4 GPa) of more than 3 orders of magnitude. An alternate approach to quantify the compression is to measure a local density or average area occupied per molecule near the surface of the substrate. The surface area per molecule in a monolayer can be obtained from adsorption experiments, for example, as shown in the early pioneering papers of Brunauer and coworkers. 17,18 Using a local density or area approach does not escape the non-uniqueness difficulty at the nano-scale we must first define a system volume or area by defining a boundary, as is routinely done in defining a pore width, for example. Thus the partition function, and hence all thermodynamic properties, are not uniquely defined at the nanoscale. [START_REF] Gray | Theory of Molecular Fluids. 2. Applications[END_REF] A third approach to measure the compression effect would be to forgo macroscopic measures such as volume and pressure, and use some nano-scale measure such as the average distance separating neighboring molecules in a particular adsorbed layer. This is easy to measure in molecular simulations; an example from the current work is shown in Figure 2 for the contact layer next to the graphite wall. When the intermolecular pair force is zero (no compression) the distance separating the centres of neighbouring molecules is 0.382 nm. However, in the actual contact layer the intermolecular separation between neighbours is smaller; the snapshot of the simulation shown in Figure 2 Each of the above approaches has limitations. Defining an average tangential pressure for an adsorbed layer has the merit that the configurational contribution is a direct reflection of the intermolecular forces that lead to the compression. However, although both the IK and H definitions seem to work well for systems that show strong wetting and for which the substrate is planar, they may not work for systems that have different constraints, for example substrates with curved surfaces or weakly wetting systems that do not show strong layering. For example, the definition is known to give spurious results near the centre for spherical geometry, 8 and we observe the same limitation when used for cylindrical geometry. Defining a local density or area per molecule is also not entirely satisfactory, as such definitions are relatively insensitive to the intermolecular forces. The adsorbed films are at high density usually, so that a small decrease in intermolecular separation corresponds to a large increase in the force. Finally, we note that there is strong evidence of adsorption compression for wetting systems on carbons, and this effect can be large in the direction parallel to the surface of the carbon substrate. Starting with the early work of Brunauer and coworkers 17,18 , experimental studies of monolayers on carbons show that the surface area per molecule is consistently lower than that estimated from bulk liquid densities for simple adsorbates that wet the surface, and that this lowering is significant, ranging from 7-19% for the adsorbates argon, nitrogen, oxygen, carbon monoxide and carbon dioxide. Similar results are obtained from scattering experiments (e.g. ref. 20). The molecular simulation results (e.g. Figure 2) also show reduced surface area per molecule, suggesting that the experimental observations may well result from adsorption-induced compression in the tangential (xy) direction. Donahue and coworkers 21 , 22 have discussed evidence for adsorption compression from both experiment and molecular simulation. In particular, they note that compression can occur in higher adsorbed layers in addition to the contact layer, and that this can lead to oscillations in the isosteric heat as the bulk phase pressure is increased. Gor and coworkers 23,24,25 have studied the elastic properties of confined fluids using both molecular simulation and ultrasonic experiments, and find that the isothermal bulk modulus is enhanced by confinement. Recently, Vasu et al. [START_REF] Vasu | Van Der Waals Pressure and Its Effect on Trapped Interlayer Molecules[END_REF] studied the effect of confinement between two graphene layers of molecules that are susceptible to deformation under pressure in the bulk phase. By comparing Raman spectra for the confined molecules with those for the same molecules in the bulk phase they estimated the effective pressure in the confined phase to be about 12,000 bar, of the same order as the effective tangential pressures obtained in molecular simulations.

The compression effect is most readily observed in molecular simulations, as shown in the above figures and in ref. 1. We note that some authors, e.g. Nijmeijer et al. 27 , have reported calculated interfacial tensions, g , that are positive for weakly wetting walls, but become increasingly negative as adsorbate-wall interaction becomes more attractive (strong wetting). Although Nijmeijer et al. do not report calculated tangential pressures, the interfacial tension from the mechanical route is given by g

= P N -P T z ( ) é ë ù û dz ò
, where P N is the component of P normal to the wall. PN is known to be only slightly affected by the strength of wetting, whereas PT grows rapidly larger as wetting increases (see Fig. 5 of ref. 1). Thus in the strong wetting regime an increasingly negative g is not unexpected.

The underlying reason for the compression effect in the simulation is the very strong attractive force exerted on the adsorbate molecules by the carbon surface. This is illustrated in Figure 3, which shows the Lennard-Jones force between a LJ argon molecule and a single carbon atom in the planar substrate, and compares that with the total force exerted on the Ar molecule by all of the substrate carbon atoms with which it interacts. Although the pair interaction force is modest, with a minimum of -5.6 piconewtons, the total force at the minimum is -81 piconewtons, about 14 times more attractive. The main reason for this strongly attractive force-field is the high density of carbon atoms in graphite, ~114 nm -3 , about 2.6 times higher than for silica (density ~ 44 nm -3 ). The strong attractive force acting in a direction normal to the surface causes the adsorbate molecules to pack tightly on the surface, resulting in a repulsive force between adsorbate molecules in the tangential direction. This compression mechanism is most readily seen in videos of molecular dynamics simulations for this system.

In summary, that the local pressure in inhomogeneous systems is not well-defined is not in contention; this has been known at least since the early papers of Kirkwood's of the 1940's and 50's. We also agree that the methods used so far to quantify the adsorption compression are not of general application to all inhomogeneous systems, and all have some unsatisfactory features. It is for this reason that we included Figure 4, showing that the compression effect is supported by the results for the radial distribution function and the pair force between neighboring adsorbate molecules. Despite this difficulty in quantification, the combined evidence from experiment, simulation and density functional theory for adsorption compression on systems where the walls are wet is convincing, and the effect can be large. It is worthy of further study and understanding.

Figure 1 .

 1 Figure 1. Local tangential pressures using the IK and H definitions, and the spatially averaged tangential pressures (eqn. 2) for four different definitions of Cij (the IK and H

Figure 2 .

 2 Figure 2. An in-plane view of the contact layer next to the carbon wall for LJ argon in a carbon slit pore of width H = 3σ at 87.3 bar, 1 bar bulk pressure, showing separation distances between neighboring molecules. The open circles show the nearest neighbor separation, 0.382 nm, when there is no compression, i.e the intermolecular force is zero. Overlap of open circles implies strong repulsion between neighbouring pairs (see Figure 4 of ref. 1)

  indicates compression in the xy plane for most of the neighbor pairs, as indicated by overlap of the open circles that represent no compression. The result of this compression effect is shown in Figure 4(b) of ref. 1, where the pair force is plotted along with the reduction in average nearest neighbor distance. The average separation distance between nearest neighbours in the contact layer is seen to correspond to a strong repulsive force.

Figure 3 .

 3 Figure 3. Comparison of the isolated pair force between a carbon atom and an argon molecule with the total force on the Ar from all carbon atoms. The dashed line marks the minimum of the pair force.