Tiberiu Catalina
email: tiberiu.catalina@insa-lyon.fr

Joseph Virgone

Frederic Kuznik

Evaluation of thermal comfort using combined CFD and experimentation study in a test room equipped with a cooling ceiling

Keywords: cooling ceiling, full-scale experimental study, CFD simulations, thermal comfort 1. Introduction

This paper reports a full-scale experimental campaign and a computational fluid dynamics (CFD) study of a radiant cooling ceiling installed in a test room, under controlled conditions. This research aims to use the results obtained from the two studies to analyze the indoor thermal comfort using the predicted mean vote (PMV).

During the whole experimental tests the indoor humidity was kept at a level where the condensation risk was minimized and no condensation was detected on the chilled surface of the ceiling. Detailed experimental measurements on the air temperature distribution, surface temperature and globe temperature were realized for different cases where the cooling ceiling temperature varied from 16.9° to 18.9°C. The boundary conditions necessary for the CFD study were obtained from the experimental data measurements. The results of the simulations were first validated with the data from the experiments and then the air velocity fields were investigated. It was found that in the ankle/feet zone the air velocity could pass 0.2 m/s but for the rest of the zones it took values less than 0.1 m/s. The obtained experimental results for different chilled ceiling temperature showed that with a cooling ceiling the vertical temperature gradient is less than 1°C/m, which corresponds to the standards recommendations. A comparison between globe temperature and the indoor air temperature showed a maximum difference of 0.8°C being noticed. This paper also presents the radiosity method that it was used to calculate the mean radiant temperature for different positions along different axes. The method was based on the calculation of the view factors and on the surface temperatures obtained from the experiments. PMV plots showed that the thermal comfort is achieved and is uniformly distributed within the test room.

The majority of air-conditioning devices function using the principle of pulsated air, where the hot air of the room is partly recycled, cooled and returned into the room. The increase of the thermal loads in the buildings, mainly due to the arrival of office computers and lighting requirements, causes the installation of airconditioning systems necessary to neutralize these loads and to create a good indoor thermal comfort. Heating, ventilating and air conditioning (HVAC) systems, which consume large quantities of energy, have become a necessity for almost all the buildings [START_REF] Ashrae | Thermal environmental conditions for human occupancy ANSI-ASHRAE[END_REF] to provide a comfortable indoor environment.

Currently, the evacuation of these quantities of latent and sensible heats is done mainly with air treated, introduced by air diffusers. To maintain comfort under these conditions, a greater volume of cooled air must be provided to the working area. Many complaints about air-conditioning systems have been claimed, especially in summer by female occupants [START_REF] Imanari | Thermal comfort and energy consumption of the radiant ceiling panel system. Comparison with the conventional all-air system[END_REF]. Disadvantages such as noise, cold-drafts, air temperature differences between the human head and foot or energy wasting in certain cases, demand the investigation of other cooling systems which can create better indoor conditions and in the same time to be energy efficiently.

The use of water to cool the surfaces of the buildings is consequently a tempting alternative solution. It is even more appealing as water cooling requires much lower flow rates and higher temperatures of the water, so an energy consumption reduction could be achieved [START_REF] Stetiu | Energy and peak power saving potential of radiant cooling system in US commercial buildings[END_REF][START_REF] Niu | Energy saving possibilities with cooled-ceiling systems[END_REF]. The chilled ceiling radiant panels are room cooling systems for placement in the ceiling zone, their cooling surfaces being connected with closed circuit heat conducting pipework containing flowing chilled water.

The main difference between cooling ceilings and air-conditioning systems is the mechanism of heat transfer.

Classic air-conditioning systems are based only on the convection, while cooled ceilings employ a combination of radiation and convection. With cold ceilings, the transfer of radiant heat occurs by a clear emission of electromagnetic waves from the hot occupants and their environments to the radiant ceiling. The hot air arriving in contact with the cooled surface is cooled below the average temperature of the room and therefore falls down the zone of occupation [START_REF] Gan | Towards a better indoor thermal environment "CFD analysis of the performance of chilled ceiling systems[END_REF]. With a cooling ceiling the temperature of water or its flow rate can be modified so that the surface temperature adapt to the desired conditions.

Ceiling radiant cooling systems can be more comfortable than conventional air cooling systems due to small vertical temperature gradient, few air movements and reduced local discomfort for occupants during long stays in cooled room environments. Catalina et al. [START_REF] Virgone | Evaluation of performances, thermal comfort and energy consumption of a reversible radiant ceiling by capillary mat: Application for the prefabricated buildings[END_REF] found after a series of simulations that the cooling ceilings offer good thermal comfort, the mean radiant temperature being an important parameter in the comfort evaluation.

Nagano and Mochida [START_REF] Nagano | Experiments on thermal environmental design of ceiling radiant cooling for supine human subjects[END_REF] analyzed the thermal comfort of five subjects in an experimental test room equipped with cooling ceilings. Their study aim was to clarify the control conditions of ceiling radiant cooling systems for human subjects in a supine position. They concluded that the mean radiant temperature for a supine human body should be used in the design of ceiling radiant cooling. It was also found that some of the subjects voted the strong radiant sensation even though the temperature difference between the ceiling and the room air was less than 5°C in these experiments.

Kulpmann [START_REF] Kulpmann | Thermal comfort and air quality in rooms with cooled ceilings-results of scientific investigations[END_REF] reported an investigation of thermal comfort in a test room equipped with a cooled ceiling surface and supplied with upward displacement ventilation air. The results showed that a good value of thermal comfort was obtained and that the temperature of the room surfaces was lower or at least equal to the air temperature in the room.

Hodder et al. [START_REF] Hodder | Thermal comfort in chilled ceiling and displacement ventilation environments: vertical radiant temperature asymmetry effects[END_REF] investigated and determined thermal comfort design conditions for combined chilled ceilings/displacement ventilation environments. A typical office room was reproduced in a laboratory test room in which the ceiling temperatures were varied. It was found that vertical radiant temperature asymmetry has an insignificant effect on the overall thermal comfort of the eight subjects that participated in this experimentation.

Hodder et al. [START_REF] Kulpmann | Thermal comfort and air quality in rooms with cooled ceilings-results of scientific investigations[END_REF] reported that the existing guidance [START_REF] Fanger | Comfort limits for asymmetric thermal radiation[END_REF] is valid with no modification, for thermal comfort design in such design environments.

Imanari et al. [START_REF] Imanari | Thermal comfort and energy consumption of the radiant ceiling panel system. Comparison with the conventional all-air system[END_REF] compared the radiant cooling systems with air-conditioning systems in terms of thermal comfort, energy consumption and costs. The thermal environment as well as human response was tested by using a small meeting room equipped with radiant ceiling panels. It was shown in their work that the radiant ceiling panel system is capable of creating smaller vertical variation of air temperature and a more comfortable environment than conventional systems.

Simmonds [START_REF] Simmonds | Practical applications of radiant heating and cooling to maintain comfort conditions[END_REF] found that the traditional design criteria such as dry-bulb temperature and operative temperature were not always sufficient, mean radiant temperature having an important influence on the comfort results. Kitagawa et al. [START_REF] Kitawaga | Effect of humidity and small air movement on thermal comfort under a radiant cooling ceiling by subjective experiments[END_REF] investigated the effect of humidity and small air movements on the thermal comfort of subjects in a climatic chamber equipped with radiant cooling panels. The results of their questionnaires showed that the most comfortable sensation vote was obtained in the conditions whose thermal sensation vote was not neutral but approximately -0.5.

The thermal comfort is, by definition, a subjective sensation; different people will express different preferences and it is preferable to use real persons to evaluate this notion. However, this estimation method requires the use of a representative large number of subjects and hence can become expensive and time consuming. The approach used in this article study is based on experimental data and computational fluid dynamics simulations in order to obtain results on the parameters that influence the thermal comfort.

Myhren et al. [START_REF] Myhren | Flow patterns and thermal comfort in a room with panel, floor and wall heating[END_REF] used Computational Fluid Dynamics (CFD) simulations to investigate possible cold draught problems, the differences in vertical temperature gradients and air speed levels in two office rooms equipped with different space heating and ventilation systems. CFD technique was also used by Rohdin et al. [START_REF] Rohdin | Numerical predictions of indoor climate in large industrial premises. A comparison between different k-e models supported by field measurements[END_REF] research study on the prediction of indoor environment in large and complex industrial premises. Stamou et al. [START_REF] Stamou | Evaluation of thermal comfort in Galatsi Arena of the Olympics "Athens 2004" using a CFD model[END_REF] have studied the thermal comfort in an Olympic Arena by using a CFD model. Calculated mean velocities and temperatures were used to determine the thermal comfort indices, predicted mean vote (PMV) and predicted percentage of dissatisfied person (PPD) and to evaluate the thermal conditions in the various regions of the Arena. CFD simulations were also found to be a good solution to determine the indoor thermal comfort in different research studies [START_REF] Nilsson | Thermal comfort evaluation with virtual manikin methods[END_REF][START_REF] Miyanaga | Simplified human body model for evaluating thermal radiant environment in a radiant cooled space[END_REF] by using CFD simulations of virtual human manikins.

For the present research article, CFD simulations were realized with the aim to analyze the air flow pattern and the velocity fields inside the test room. The thermal comfort was evaluated by using the PMV index [START_REF] Iso | Moderate thermal environments: determination of the PMV and PPD indices and specification of the conditions for thermal comfort, International Organization for Standardization[END_REF] that was calculated based using the data obtained from the experimental and CFD analysis. For the thermal comfort study, experimental results on the air temperature, mean radiant temperature and relative humidity were used, while the CFD study gave valuable results on the air velocity and air pattern. The particularities of this research study are related to a real-size controlled experimental set-up of using radiant ceiling panels equipped with capillary tubes. Another particularity of this article was the use of CFD simulations to examine the air flow pattern and to obtain conclusions on the eventual local discomfort that could appear under a cooling ceiling.

Experimentations

Test room presentation

Studying the influence of certain indoor parameters on the thermal comfort when using a cooling ceiling could be a delicate issue, an experimental investigation being mandatory. The experimental cell MINIBAT (INSA de Lyon, France), presented in Figure 1, is composed of two identical enclosures (Cell 1 and Cell 2), which dimensions are 3.10, 3.10, 2.50 m according to the coordinate directions (x, y, z). A single glazed façade isolates the test cell 1 from a climatic chamber where the air temperature is controlled by the means of an air-treatment system. The climatic chamber temperature can vary between -10°C and 40°C as a function of time. Using a battery of 12 spotlights, of 1000 W each, it is possible to simulate an artificial sunning for the test Cell 1. The spectrum of the gas-discharge lamps with metal halide is similar to the sun one. For this study, only in the test cell 2 were installed the radiant ceiling panels and were taken measurements. The cooling ceiling surface temperature was varied during the experiments, but the rest of the five exterior walls of the MINIBAT were kept at a constant temperature of 26°C, by the mean of a thermal guard controlled with a precision of ±0.5°C. The air distribution in the thermal guard was made by a ventilation network equipped with lateral diffusers and placed in the upper part of the controlled volumes. This distribution was chosen to avoid dead zones and to create a uniform air flow in the technical spaces. In Table 1 are presented the main elements of the MINIBAT envelope and in Table 2 are summarized the physical characteristics of the materials. For the data acquisition it was used a multiplexer-multimeter with more than 100 connections. The data acquisition central was linked using the software Labview with a computer, located outside the test facility. The time step chose for the measurement was 45 seconds and the total duration of the experiments was extended on several days.

Radiant cooling system description

In the test cell 2 were installed 6 radiant cooling panels (1.5 x 0.9 m -dimension for one panel), covering the entire ceiling area. The particularities of the installed system are the capillary tubes which are in contact with the interior part of the panels. With an exterior diameter of 3 mm and spaced between them by 1.5 cm, the capillary tubes offer good temperature uniformity along the panels' surfaces by comparison to similar systems.

The very small thickness of the heat exchangers capillary pipes allows them to be embedded into the surface of the walls, ceilings or floors. Passive elements of construction are thus transformed into heating and/or cooling surfaces by the use of the energy transferred between two radiant bodies. Each mat element was connected in parallel to the overall cold water collector in order to induce a homogeneous distribution in the ceiling's surface temperature. The water flowed through the mat at a temperature between 17°C and 20°C, lower values would have increased the risk of condensation. Such temperatures would not lead to condensation of moisture from air, but could offer opportunity for substantial energy savings since water at such temperatures could be provided without using refrigeration or active cooling. For our experiment it was used an air-water heat pump to obtain the desired water temperature for the cooling ceiling, along with a pipework system and a water pump.

Measurements

The test cell indoor air temperature was measured using 6 K-type thermocouples (Nickel-Aluminum) at different heights starting from 0.4 to 2.3 m and placed in different positions of the room. The test facility walls and floor have been equipped with K-type thermocouples in order to measure the internal/external surface temperatures; nine measurements were taken for each face of the room, a total of 108 surface temperatures being scrutinized.

The radiant ceiling temperatures were measured using Pt100 probes in 9 different points. The measurements error for the thermocouples was of ±0.5°C and for the Pt100 probes of ±0.3°C. The water temperature in the system was also measured at inlet/outlet position using Pt100 probes. The indoor humidity was measured using a sensor with a precision of ±2.5% for measurements of relative humidity between 10 and 90%. The black globe temperature was analyzed using a thermocouple positioned in a black globe with a diameter of 15 cm at a height of 1.2 m. The test cell was divided in four symmetric spaces, the analyzed measurement area being one of these parts (see Figure 2). The thermal guard temperature was kept at 26°C during the whole experimental campaign, but the ceiling temperature varied from 16.9°C to 18.9°C. The air temperature, relative humidity and globe temperature were measured in different locations of the rectangular grid of 0.43 x 0.43 m presented in Figure 2.

The measurements were taken also on different heights, a total number of 96 measurements (on x, y, z-axis) being realized. The experiments were conducted under natural convection conditions and no humidity sources were added. All of the experiments were carried out under steady-state conditions.

Experimental results

Vertical air temperature profile

Three experimental cases were carried out to investigate the effect of radiant cooling ceiling temperature (T RCC ) on the indoor environment. Figure 3 shows the vertical distribution of air temperature in the center of the test room. If the temperature difference between the head and the feet is sufficiently large, a person can have local warm discomfort and/or cold discomfort even though the overall average is thermally neutral. Nominal vertical air temperature gradients between 0.4 and 2.1 m heights were found to be around 0.71-0.77°C/m for the cooling ceiling at temperatures of 16.9°C, 17.6° and 18.9°C. The obtained experimental data are more than satisfied, in terms when the comfort standards call for a maximum temperature difference between head and feet of 3°C, regardless of the operative temperature.

Surface temperature

The MINIBAT test cell was equipped with 108 thermocouples to measure the surface temperature on the interior and exterior part of the room. It wasn't notice any major wall temperature fluctuation due to air treatment system of the thermal guard. The impact of surface's temperature on the human thermal comfort is considerable and a particular attention have be paid to this aspect. A 0.5°C difference was observed between the lower and the upper part of the walls (see Figure 4).

Mean radiant temperature

Even though they are not in direct contact with the body, cold surfaces (like radiant cooling ceiling) still greatly affect human perception of temperature. With radiant cooling, the energy balance on the human body is different than without the cooled ceiling in several ways. First, the heat losses from the human body by radiation are increased from about 35% without radiant cooling to 50% with radiant cooling [START_REF] Kulpmann | Thermal comfort and air quality in rooms with cooled ceilings-results of scientific investigations[END_REF]. Likewise, the heat losses due to convection decrease from about 40% without the chilled ceiling to 30% with radiant cooling.

The globe temperature, which integrates the effect of air temperature, radiation and air movement, was measured at z=1.2 m for different locations of the (x, y) plane. It was found that the globe temperature had lower values than the air temperature with 0.6°C to 0.9°C (see Figure 5).

Mean radiant temperature (MRT) is the most important parameter governing human energy balance, with the strongest influence on thermo-physiological comfort indexes such as PMV, its meaningful scope being the thermal comfort. To evaluate the MRT for all the axes and heights we used a radiosity method. Considering a spherical point located at a spatial position in a room, the radiative balance equation can be written for this point.

Taking into account the reciprocal relation of the view factor, it is obtained, by using a linearization of the radiative exchanges between the spherical point and the walls and neglecting the internal reflections, the thermal balance of the point [START_REF] Inard | [END_REF]:

4 1 4 0 ∑ = = N j Sj j T F MRT (1)
We suppose to have only rectangular surfaces for the different envelope components to consider. F oj , the diffuse view factor between a spherical point and a rectangular surface can be calculated with the Eq.4 according to the notations of Figure 6. In the Eq.3 the surface temperatures (T Sj ) temperature where obtained from the experimental campaigns.

                    + + ⋅ = 2 2 2 4 1 c b a c ab arctg F oJ π (2)
A comparison between the measurements taken on the mean radiant temperature and the results obtained using the radiosity method is presented in Figure 7. The results correspond to the position (x 1 ,y 1 ,z 3 ) to (x 4 ,y 1 ,z 3 ) as presented previously in Section 2.3. The results obtained using the radiosity model were in good agreement with the data obtained from measurements, a maximum difference of 0.4°C being noticed, value which can be explained by the measurements errors from the thermocouple sensors. The results of MRT calculations are presented as two temperature fields at position x 1 ,y 1 and x 4 ,y 1 on the y-axis. The analysed cases correspond to the position in the center of the test room and in the proximity of the wall. Figure 8 shows the mean radiant temperature distribution calculated using the Eq.3 on the two analysed cases. The effect of the radiant cooling ceiling on the radiant temperature is more noticeable for the x 1 ,y 1 position and is reduced when approaching the walls, due to radiant effect of the walls temperature.

Indoor humidity levels

One of main inconvenient of the radiant cooling ceiling is that it's presenting a relative high risk of condensation on the chilled area. Even if the indoor humidity has a small effect on the global thermal comfort, for the cooling ceiling systems it's a crucial design parameter. Usually the dew point temperature is calculated for each zone by monitoring all the time the ambient air temperature and the relative humidity. The inlet water temperature is controlled individually with values higher than the dew point of the ambient air. Therefore, if a risk of condensation is present, the water temperature is raised. A good solution to reduce the risk of condensation is to dehumidify the air but this solution will increase considerably the price of the system and it cannot be used alone, a control system being necessary for perfect safety against the condensation. For our experiment the relative humidity was kept at values between 42% and 43% (see Figure 9), so the risk of condensation on the ceiling was minimized.

Mathematical model of airflow

During the last years, Computational Fluid Dynamics (CFD) technology has progressed and made it possible to analyze complex situations where HVAC systems were simulated and compared with real data obtained from experiments [START_REF] Kuznik | Experimental and numerical study of a full scale ventilated enclosure : comparison of four two equations turbulence models[END_REF].

The mean air flow modelling is divided in two parts: a) mathematical model (governing equations, turbulence model, boundary conditions, et.) and b) numerical solution (discretisation, numerical scheme, solution convergence, etc.)

The CFD model was based on the description of the experimental test room. The airflow pattern and temperature distribution in the facility are governed by the conservation laws of mass, momentum and energy. The flow is assumed to be steady state, three-dimensional, incompressible and turbulent. The buoyancy effect is invoked in the momentum equation, k and epsilon. The Boussinesq approximation hypothesis is used for the buoyant force term. The radiation heat transfer is not integrated in the CFD model due to the constant temperature of surrounding surfaces. The numerical model is based on commercial CFD code called STAR-CCM+ [START_REF]Star CD 3.1B. User guide and methology[END_REF]. The Rayleigh-number has been calculated and its value of around 10 10 showed that the airflow is turbulent.

Turbulence modelling

For this study the turbulence is modelled with the AKN Low-Re K-ε model [START_REF] Abe | A new turbulence model for predicting fluid flow and heat transfer in separating and reattaching flows -1. Flow field calculations[END_REF]. This turbulence model was evaluated numerically stable with correct precision on the results [START_REF] Hashimoto | Numerical study on airflow in an office room with a displacement ventilation system[END_REF]. Hashimoto Y. [START_REF] Hashimoto | Numerical study on airflow in an office room with a displacement ventilation system[END_REF] indicated that the standard k-ε turbulence model could be applied for a fully turbulent flow but it is inappropriate for a buoyant flow. In his research article he studied the airflow in an office room with displacement ventilation system using a low-Re number k-ε turbulence model with a damping function to reduce the turbulent viscosity near a wall (AKN model).

The AKN model introduced Kolmogorov velocity scale, instead of the friction velocity, to account for the nearwall and low-Reynolds number effects in both attached and detached flows. The fluid is considered like an incompressible one with density computed via the ideal gas law considered as varying only with temperature.

The model used in this study can be written in a general form as:

φ φ φ φ φ ρ ρ   ∂ ∂ ∂ ∂ + - Γ =   ∂ ∂ ∂ ∂     j j j j u S t x x x ( 3
)
where φ is the variable, φ Γ the effective coefficient and φ

S

the source term of the equation. The expressions of the variables for the turbulence model used in this paper are summarized Table 3.

Boundary conditions

The CFD model boundary conditions were taken from the experimental data results of the walls', cooling ceiling and floor surface temperature. Table 4 resumes the surface boundary conditions used for the CFD simulations based on the three experimental series where the cooling ceiling temperature was varied from 16.9°C to 18.9°C.

The wall boundaries have been modelled using the no-slip condition with constant wall temperature. The low-y+ wall treatment was applied for our model. The wall treatment assumes that the viscous sublayer is well resolved and thus wall laws are not needed.

Discretization

The mesh is designed using the pre-processor STAR-DESIGN [START_REF]Star CD 3.1B. User guide and methology[END_REF] and the discretization of the computational domain is achieved by means of an unstructured mesh. In the boundary layer next to a non-slip wall, there are high gradients within a small region, so to capture these gradients accurately it was necessary to have fine mesh spacing normal to the wall. The volume mesh process consisted of several steps, including surface improvement, subsurface generation and then the actual interior volume meshing. The grid contains polyhedral elements, the final mesh being composed of 116,762 cells. The volume range was between 3.386E-6 to 1.5351E-3 and the minimum distance between centroids of neighbour cells was of 3.374e-3. A grid dependency analysis was conducted to ensure that the resolution of the mesh was not influencing the results.

Numerical scheme

The solution method is based on the following main hypothesis: the diffusion terms are second-order centraldifferenced, and the second-order upwind scheme for convective terms is used to reduce the numerical diffusion.

The velocity-pressure coupling method is the SIMPLE algorithm. The multigrid scheme allows to accelerate the convergence as our model contains a very large number of control volumes [START_REF] Kuznik | A Second Order Turbulence Model for the Prediction of Air Movement and Heat Transfer in a Ventilated Room[END_REF].

CFD simulations results

Grid independence analysis

Even if a convergent solution was obtained an additional step was necessary to achieve good accuracy. Chen et al. [START_REF] Chen | A procedure for verification, validation, and reporting of indoor environment CFD analyses[END_REF] recommends a verification of the model results by systematically refining the grid size realized generally by doubling the grid number and compare the two solutions. A grid independence test was performed and the data results showed that the mesh density chosen was sufficient for a grid independent solution of velocity fields, no noticeable error being observed. Figure 10 shows the velocity results difference between the two analysed grid mesh (1) 116,762 cells and (2) 233,524 cells.

Comparison of the experimental data to CFD simulations

The model validation is important to ensure adequate results and was a necessary step before analyzing the air flow pattern obtained from the CFD simulations. The validation was focused on the comparison of air temperature for the three experimental series. The simulated results were analyzed and compared with the experimental data on the position (x1, y1) along the z-axis according to the analyzed experimental configuration presented previously in the article.

The analyzed data points are the most relevant in terms of temperature distribution in the test room. In Figure 11 the air temperature distribution is confronted with the experiment data. The comparisons between the experimental and the numerical data confirm the validity of the numerical model used in the CFD simulation; the differences between the CFD and experimental results could be explained by the measurements errors of ±0.5°C.Furthermore the accuracy of the results was analyzed by a sensitivity study of the air velocities values and air pattern flow to a change in the cooling ceiling temperature. This temperature was modified accordingly to the measurements errors of the Pt probes used in the experimental campaign (±0.3°C). The conclusions of this sensitivity analysis showed that the air flow pattern is not modified and the differences between the air velocities were insignificant.

Air velocity fields

The air speed could be an important parameter for the human thermal comfort; an increased air speed will aid the evaporation of sweat thus leading to a cooling effect, particularly if loose clothing is worn. However, if the air velocity is too high it may cause discomfort and a sensation of draughtiness. The thermal comfort cannot be ameliorated by an increase of air velocity if the temperature or humidity is too high. Mean air speed should be less than 0.15 m/s in the occupied zone to prevent a feeling of discomfort for the most sensitive. The CFD technique was found a valuable technique to complement the experimental results, allowing a detail access to velocity fields on different planes. In Figure 12, the flow fields at the two analyzed position (x1, y1) and (x4, y1) on the yz-axis are plotted for the experimental study where the radiant cooling temperature was set at 17.6°C.

The hot air heated by the floor rises due to Archimedes effect; once the ceiling is reached, the air is cooled and is then goes down. It can be noticed that the air speed is higher than 0.2 m/s in the feet and ankles zone and may cause problems for the occupants. Because the vertical walls temperatures are not strictly equal, the air flow is not symmetric in the test room. The results showed that the air velocity levels were still below the standards limits in most of the occupied zones.

Thermal comfort analysis

The use of standardized scales to gather ratings of thermal sensation from groups of individuals has permitted the researchers to identify the vital elements of the thermal environment and to evaluate the average response of the subjects. Thermal sensation ratings given by real subjects are recommended, but these studies are expensive and demand a large number of subjects in order to get correct results. For our study the predicted mean vote PMV and the predicted percentage of dissatisfied PPD [START_REF] Myhren | Flow patterns and thermal comfort in a room with panel, floor and wall heating[END_REF] were used to evaluate the thermal comfort. Both indices have been used with success during the past decades and heir valid range extends to our case of air conditioning with cooling ceiling panels.

The PMV model assumes a relation between optimal thermal conditions, using the steady state heat balance equation for the human body and thermal comfort ratings from panels of subjects. The PMV is related by a 7point thermal sensation scale of a group of people exposed to a certain environment to the calculated result of the basic heat balance equation. The PMV range values are within -3 and +3, where the positive values stand for hot feeling and the negative for cold feeling, 0 being a neutral value.

The elements necessary to calculate the PMV are the air temperature, the mean radiant temperature, the air velocity, the relative humidity and two subjective parameters which are the human metabolism rate and the clothing insulation of the person. For this study these parameters were obtained by the mean of an experimental study, by calculations and by a CFD analysis. Figure 13 shows how the data obtained from the different studies are connected with the thermal comfort analysis. The main purpose of the CFD simulation was to give important information on the air flow pattern and on the air velocity fields. The experimental study was necessary to obtain data on relative humidity and air temperature, but also to provide the boundary conditions mandatory for the CFD analysis. The mean radiant temperature was calculated by using a radiosity method and had its support on the results of surface temperatures provided by the experimental tests.

In Figure 14 the PMV distribution on the planes (x 1 ,y 1 and x 4 ,y1 on the yz-axis) are investigated. Figure 15 presents the PMV index when the radiant cooling ceiling is set to 17.6°C, for different metabolism rates and clothing insulations. The hypothesis for the analysis were that the metabolism (M) corresponds to an office work and its value was 1.2 met (70 W/m²) and the clothing insulation (I cl ) of 0.58 clo (0.089 m²K/W). The effect of the radiant ceiling is more intense in the center part of the room where the mean radiant temperature is lower with an overall thermal comfort which is acceptable and with values between -0.2 and 0.2. The two plots of PMV demonstrate that it can be obtain a uniform thermal comfort distribution in the room when using the chilled ceiling to cool the indoor air.

Based on the data obtained from the studies, the PMV/PPD were calculated for different heights and for different levels of metabolism rate and clothing insulation. Even for higher levels of I cl the thermal comfort is still good, exceeding the comfort range described in ASHRAE 55 (-0.5÷0.5) [START_REF]Thermal environmental conditions for human occupancy[END_REF] only in the case where the metabolism rate takes higher values than 1.3 met.

Conclusions

The air cooling using chilled ceiling panels was found to be an interesting alternative to traditional airconditioning system. The research work was based on a fully experimental study which had the support of an extended analysis on the air flow pattern by using the CFD technique. The boundary conditions necessary for the CFD simulation were taken from the experimental data and then a validation on the air temperature was realized.

The air velocity fields obtained from the CFD study showed that a local discomfort at feet/ankle zone was observed but with good values of air velocity for the rest of the test room. The mean radiant temperature, a necessary element for the PMV calculations, was evaluated using a radiosity method by calculating the view factors in each point of the analyzed planes and with the data from the experiment. The experimental results obtained for different chilled ceiling temperature showed that with a cooling ceiling the vertical temperature asymmetry is less than 1°C, value which corresponds to the standards recommendations. The globe temperature was measured and compared with the air temperature, a maximum difference of 0.8°C being observed.

The relative humidity was kept at normal values of 43% in order to avoid the condensation that could have been produced on the cold surface of the ceiling. The thermal comfort was studied using the PMV index scale. The elements necessary to calculate the PMV were taken from the two studies, CFD and experiment. The PMV distribution plots showed that thermal comfort is achieved and is uniformly distributed no matter the position in the room for the different planes in the test room. The results reported in this paper confirm that the cooling ceiling has advantages like low vertical air gradient and that the thermal comfort is obtained even for higher metabolism rates or clothing insulation.

Nomenclature

  k [m²s -3 ] ρ mass density [kgm -3 ] τ turbulent time scale [s]


Jona-Lasinio Mattia

Ludovic Pricoupenko

Three Resonant Ultra-Cold Bosons: Off-Resonance Effects

Keywords: numbers: 34, 50, Cx 03, 65, Nk 03, 75, Kk 05, 30, Jp

We solve a finite range two-channel model for three resonant identical bosons. The model provides a minimal description of the various magnetic Feshbach resonances in single species ultra-cold bosonic systems, including off-resonant scattering. We obtain important insights into the interpretation of seminal experiments: the three-body recombination rate measured in Sodium and the Efimov resonances observed in Caesium. This approach quantifies non universal effects appearing for a finite magnetic field detuning.

One of the main issues in current ultracold physics is to achieve highly correlated quantum gases. In these studies the magnetic Feshbach resonance is a crucial ingredient: it permits to tune the strength of the two-body interaction measured by the s-wave scattering length to an arbitrarily large value, while the atomic density remains constant. However for dilute bosonic gases, such highly correlated states are very unstable as a consequence of three-body recombinations into deep molecular bound states [1,2]. Therefore, in the last decade an impressive effort has been made both theoretically and experimentally, for a deep understanding of three-body properties in these systems. Universality concepts borrowed from Nuclear physics together with the specificity of ultracold atoms where an explicit energy scale separation occurs in scattering processes have led to unsuspected physical insights in this domain [3,4]. For example, the so-called Efimov states have been observed for the first time in a resonant ultracold Caesium gas [5]. For scattering lengths much larger than the range of inter-atomic forces, Universal Theory [4] permits a powerful analysis of three-body properties without any knowledge of the short range details of the real interatomic forces. For narrow resonances like the one observed in Ref. [2], it is possible to fully determine the Efimov spectrum at resonance using the effective range approach [6,7]. Despite their success, these approaches are not designed for a description of three-body properties at finite detuning, where off-resonant scattering effects come into play. Experiments on Efimov states in Ref. [9] clearly exhibit some non universal behavior which cannot be taken into account by the Universal Theory. Moreover, the peak in three-body losses in Ref. [2] occurs for relatively large magnetic detuning where the effective range approach is inoperant. In this letter, we solve the three-boson problem by using a two-channel model including the short range character of interatomic forces. The model describes all the various type of magnetic Feshbach resonances: broad, narrow or in the neighborhood of a shape resonance. Results of the model compare quantitatively with experiments of Ref. [2], highlighting the importance of off-resonant effects in narrow resonances. Concerning the experiments in Refs. [5,9], our results show that the observed violation of universality follows from the fact that the resonance in the atom-dimer scattering appears in a domain where the range of interatomic forces is not negligible as compared to the scattering length.

We first introduce the finite range two-channel model. The two channels refer to the 'open' channel populated by atoms which are identical ultra-cold bosons of mass m, and to the 'closed' channel where fundamentals entities are couples of two tightly bound atoms that we here call 'molecules'. The model takes a simple form in the second quantized form, where the operator a k annihilates an atom in the open channel with momentum k, while b k annihilates a molecule of wavevector k in the closed channel. Both a k and b k obey standard bosonic commutation rules

[a k , a † k ′ ] = [b k , b † k ′ ] = (2π) 3 δ(k -k ′ )
, corresponding to the choice r|k = exp(ik • r) for the plane wave state. Any other commutator vanishes. The Hamiltonian is similar to the ones introduced for two-component or fully polarized fermions in Refs. [10,11]:

H = dk (2π) 3 ǫ k a † k a k + ǫ k 2 + E mol b † k b k + g 0 2 dkdKdk ′ (2π) 9 χ * k χ k ′ a † K 2 -k ′ a † K 2 +k ′ a K 2 +k a K 2 -k + Λ dkdK (2π) 6 χ * k b † K a K 2 -k a K 2 +k + h.c. . (1)
The first two terms in Eq. ( 1) are the kinetic operators in the open and closed channel respectively:

ǫ k = 2 k 2 /(2m)
, and E mol is the internal energy of the molecular state in the closed channel, defined with respect to the zero energy in the open channel. The magnetic tunability of the interaction strength is due to the fact that E mol is an affine function of the magnetic field B with the slope δµ, where δµ is the difference between the magnetic moments for an atomic pair in the open and closed channel. The second line in Eq. (1) mimics the interatomic force in the open channel which is responsible for the background scattering. In a real system it is characterized by an attractive van der Waals tail, and the van der Waals constant C 6 sets up the short range scale of this pairwise interaction: R vdW = 1 2 (mC 6 / 2 ) 1/4 . In the present model, for simplicity we choose a separable interaction, with a Gaussian weight χ k = exp(-k 2 b 2 /2) imposing a cut-off at short distances: b ≡ O(R vdW ). The last term in Eq. ( 1) describes the coupling between the two channels and models the Feshbach resonance mechanism (we choose Λ ∈ R). The inter-channel coupling is also taken into account via the same Gaussian weight χ k as in the open channel interacting term.

We now determine the different parameters of the model from measured two-body properties. For this purpose we solve the two-body scattering problem in the center of mass frame. For an incident plane wave of energy E = 2 k 2 0 /m and momentum k 0 , the two-body state is a coherent superposition of one molecule plus two atoms:

|Ψ = (βb † 0 + dk (2π) 3 A k a † k a † -k )|0
, where A k is the atomic wave function in the form:

A k = (2π) 3 δ(k -k 0 ) + 4π 2 m f (E) 2ǫ k -E -i0 + , (2)
and f (E) is the scattering amplitude. Remarkably, the s-wave scattering length (denoted by a) obtained in this two-channel model can be exactly identified with the following expression which is known to be very accurate in the vicinity of a Feshbach resonance [12]:

a = a bg 1 - ∆B B -B 0 . (3)
In Eq. ( 3), a bg is the background scattering length, and ∆B is the width of the resonance located at the magnetic field B 0 . One finds a bg = bg 0 √ π/(g 0 -g c 0 ), with g c 0 = -4π 3/2 2 b/m, an energy detuning from resonance ν = E mol -2Λ 2 /g 0 + δµ∆B where ν = δµ(B -B 0 ) and a resonance width ∆B = 8π 2 Λ 2 a bg /(mg 2 0 δµ). This last relation implies that the energy width δµ∆B always has the same sign as the background scattering length a bg . We checked that this property is indeed verified for the various resonances reported in Ref. [12]. Using these parameters, the scattering amplitude admits a simple expression for E = -2 q 2 /m < 0 (k 0 = iq and q > 0):

1 f (E) = qerfc(qb) - e -q 2 b 2 a bg 1 - δµ∆B E -ν + δµ∆B . (4)
We use experimental or theoretical spectroscopic data on the two-body bound states as a way to choose a precise value for b (in absence of such data, we arbitrarily set b = R vdW ). In order to avoid any confusion the two-body bound states are denoted by 'dimers' and are distinct from the molecular state. Their binding energies

E = -E dim = -2 q 2 dim /m < 0 are poles of f (E).
In what follows, we briefly sketch their spectrum as a function of the energy detuning ν. For ν < 0, there is a branch terminating at zero energy for ν = 0 -: this branch results from the interchannel coupling and is denoted below as the Feshbach dimer's branch. Furthermore, in the case a bg > b √ π, another branch exists for all possible values of the detuning. Away from the Feshbach resonance and for ν > 0, it results from the direct coupling in the open channel and we denote it as the 'background dimer' branch. For decreasing values of ν there is an avoided crossing between the two branches. We checked for several resonances that a choice of b of the order of R vdW permits to describe the lowest dimer's branch over a wide range of magnetic field detuning. In Fig. (1), we compare the results of the model with the experimental data of the resonance located at B 0 ∼ -11.7 G for Caesium [8,9] (∆B ≃ 28.7 G, a bg ≃ 1720a 0 , R vdW ≃ 101a 0 and δµ ≃ 2.3|µ B | [12]). The system is in the vicinity of a shape resonance (a bg ≫ R vdW ) and the spectrum displayed corresponds to the branch of the background dimer. We found very good agreement for b = 0.7R vdW . We now turn to the central part of the present work by investigating the problem of three interacting bosons. The general three-body eigenstate of Eq. (1) of energy E in the total center of mass frame, is a coherent superposition of three atoms plus one atom and one molecule:

|Ψ = dK (2π) 3 dk (2π) 3 A K,k a † K 2 +k a † K 2 -k a † -K |0 + dK (2π) 3 β K b † K a † -K |0 . ( 5
)
In Eq.( 5) A K,k and β K are respectively the atomic and the atom-molecule wavefunctions. For E > 0, the atomic wavefunction contains an eigenstate of the atomic kinetic operator A

K,k , and for all E ∈ R, A K,k can written as:

A K,k = A (0) K,k + Λχ k β eff K /(E rel -2ǫ k + i0 + ), (6)
where

E rel = 2 k 2 rel /m = E -3 2 K 2 /(4m)
is the relative energy between a pair of atoms and the third atom. For convenience, we introduced in Eq. ( 6) the effective atommolecule wavefunction β eff K :

β eff K = β K × (E rel -ν + δµ∆B)/(E mol -ν + δµ∆B). (7)
The three-boson problem for this model is solved by finding the solutions of the integral equation which is deduced from the stationary Schrödinger equation:

m|χ k rel | 2 β eff K 4π 2 f (E rel ) -2 dk (2π) 3 β eff k χ * K 2 +k χ K+ k 2 ǫ k + ǫ K + ǫ K+k -E -i0 + = - dk (2π) 3 χ * k Λ A (0) K,k + 2A (0) -K 2 +k,-3K 4 -k 2 . (8)
Remarkably all the two-body physics, contained in the scattering amplitude f (E), appears in the diagonal part of Eq. ( 8). Nearby a resonance, in the limit where |a| → ∞, Eq. ( 8) converges asymptotically in the lowenergy regime (i.e. for kb ≪ 1 and Kb ≪ 1) toward the so-called Skorniakov Ter-Martirosian equation of Ref. [13], while the high energy limit of the integral kernel acts as an ultra-violet cut-off and the Thomas collapse is avoided [14]. As the magnetic detuning |B -B 0 | is increased, off-resonant effects come into play and the model can be used to quantify deviation from the universal theory [3,4]. To this end, we numerically solved Eq. ( 8) using two independent codes for several Feshbach resonances. We report in the following our results concerning the resonances studied experimentally in Refs. [1,2,5,9]. First, we investigate the spectrum of trimers made of Caesium atoms for the Feshbach resonance at -11.7 G. Trimers are obtained by searching negative energy solutions (or E < -E dim if a shallow dimer exists) of Eq. (8) where A (0) = 0. At resonance (ν = 0), we recover the Efimov spectrum: E n = E 0 e -2πn/s0 , (n ∈ N) where E 0 is the energy of the lowest Efimov state [of the order of -2 /(mR 2 vdW )] and s 0 = 1.00624 [3,4]. In Fig.

(2) we plot the spectrum as a function of the magnetic field. Our model predicts four important features: (i) The existence of two Efimov branches extending in the regions of positive and negative B. Each branch is continuous through the formal limit |B| = ∞. Other Efimov branches are located in a small interval near B 0 ; (ii) The threshold of the first Efimov branch at 10.3 G observed in Ref. [5] at 7.8 G; (iii) The two trimer branches hit the background dimer branch in a non universal region where short range details of the interatomic forces are not negligible -see inset of Fig. (1)-and where our model gives qualitative informations only; (iv) In experiments reported in Ref. [9], the magnetic field was decreased from a large detuning and an atom-dimer resonance loss was found at ∼ 25 G. The present model shows that the observed threshold corresponds to the second trimer branch. The disagreement between theoretical and experimental thresholds follows from item (iii). Moreover the asymptotic behavior of the wave function at distances r ≪ a is a crucial ingredient in Universal Theory and follows from taking the 'unitary approximation' of the scattering amplitude

f ∼ -1/ik at intermediate momentum 1/a ≪ k ≪ min(1/R vdW , 1/r e ),
where r e is the effective range. However, at the observed threshold ∼ 25 G, a/R vdW ∼ 4 and from Eq. ( 4), a/r e ∼ 2. Hence, in this detuning region there is no clear separation of momentum scale and the 'unitary approximation' for intermediate momentum is not correct. This fact gives an important insight in the deviation from prediction of Ref. [4] concerning the ratio between the scattering lengths at the two observed thresholds.

In this last part, we study the recombination of three incoming atoms of vanishing total energy into one twobody bound state and one atom. This is the main process responsible for the short life time of a resonant BEC in a dipolar trap. Experimentally, it is measured via the atomic loss rate which is defined for N atoms trapped in a cubic box of size L by:

Ṅ = -α rec N (N -1)(N -2)/L 6 , (9)
where α rec is the three-body recombination constant.

The atomic loss rate plotted as a function of the magnetic field detuning exhibits a large peak centered at the resonance. We limit our analysis to the standard regime where for a given detuning there exists at most one shallow dimer only (a bg is of the order of R vdW or nega-There are two distinct regimes in the recombination process: (a) In the regime of negative energy detuning (ν < 0) the formation rate of Feshbach dimers is the dominant loss mechanism. This dimer is an eigenstate of the model Hamiltonian in Eq. ( 1) and the recombination constant can be computed exactly; in what follows, it is denoted by α Fesh rec ; (b) In the regime ν > 0 there is no shallow dimer and the deep bound states populated in inelastic scattering processes are not described by our model FIG. 3: Three-body recombination constant αrec as a function of the magnetic field B for the narrow Feshbach resonance in 23 Na at 907 G. The three different symbols correspond to different ramp speeds of the magnetic field across the resonance in Ref. [2]. Solid line: prediction of the two-channel model, dashed line: prediction of the effective range model in Ref. [6]. The vertical line at resonance is due to the existence of narrow oscillations of αrec also predicted in Refs. [6,7].

Hamiltonian. In this case, we denote the recombination constant by α deep rec . In both regimes the source term in Eq.( 6) is A (0) K,k = (2π) 6 δ(K)δ(k) and corresponds to the incoming atomic wave function in the three-body scattering process. In the first regime [case (a)], the outgoing dimer and atom created by inelastic scattering have a relative momentum 2 √ 3 q dim . The dimer formation manifests itself as an outgoing wave in the atom-molecule wave function, and the wave function β K has a pole at K = 2 √ 3 q dim + i0 + ; we denote the residue of this pole by γ. In order to evaluate α Fesh rec , we enclose the three incoming atoms in a fictitious box of arbitrary large size L and impose periodic boundary conditions. The wave function is then deduced from Eq. ( 5) with |Ψ box ≃ |Ψ /L 9/2 . The rate of molecule formation is obtained in configuration space by computing the total flux of the probability current associated with the atom-molecule wave function (relative particle of reduced mass 2m/3) through the box surface. Since the flying dimer has a non vannishing probability (p closed ) to be in the closed channel, we di-vide this total flux by p closed thus obtaining the atomic loss rate Ṅ . By comparing this expression with the case N = 3 in Eq. ( 9), one obtains:

α Fesh rec = 2 √ 3 q 3 dim |γ| 2 /(9πmp closed ). (10)
In the regime ν > 0 [case (b)], the recombination rate can be evaluated on a qualitative basis only. The idea is to compute in the fictitious box of size L the probability P < of finding the three atoms or the atom and the molecule in a volume of the order of R 3 vdW . This calculation is performed in configuration space by using Eq. ( 5). Since R vdW gives the typical size of deep bound states, the loss rate is obtained by dimensional analysis with Ṅ ∝ -P< To conclude, we presented a rather simple formalism to capture the main physical features in resonant threebosons systems. The short range details of interatomic forces are described by one parameter only, so that we avoided the complexity of more detailed models [15,16].

FIG. 1 :

 1 FIG. 1: Comparison between experimental dimer's energies in133 Cs obtained at Innsbruck (black losanges)[8] and results of the two-channel model for b = 0.7R vdW (black line). Inset: detail of the deviation of the spectrum from the universal law E dim = 2 /(ma 2 ) (dashed line).

FIG. 2 :

 2 FIG. 2: Spectrum of bound states in 133 Cs (B0 ∼ -11.7 G). Solid lines: trimers; dashed line: dimers; dotted horizontal line: "High energy" limit -2 /(mR 2 vdW ); doted vertical lines: position of the Efimov resonances observed in the experiment (at 7.5 G and 25 G resp.); Crosses: thresholds computed with the two-channel model (at 10.3 G and 23.2 G resp.).

mR 2 vdW, 2 vdW.

 22 and finally the recombination constant is estimated by: α deep rec = L 6 P< mR We applied this formalism for the narrow resonance [i.e. δµ∆B ≪ Na at B 0 ≃ 907 G [1] (∆B ≃ 1 G, δµ ≃ 3.8|µ B |, a bg ≃ 63a 0 and R vdW ≃ 44.5a 0 [12]). Fig. (3) shows a dramatic agreement of our results (for b = R vdW ) with the experiments in Ref. [2]. Three-body properties in narrow resonances are usually described within the effective range approximation by using the two parameters a and R ⋆ = 2 /(ma bg δµ∆B) [6]. As shown in Fig. (3), this latter approach gives reasonable results very close to the resonance only. In the zero range limit (b → 0), when the parameters a and R ⋆ are held fixed, and a bg ≡ 0(b) vanishes (but a bg = b √ π), the present model coincides with the effective range approach exactly. Therefore, in realistic situations where b and a bg are both finite, this two-channel model quantifies consistently off-resonant effects. We verified that α rec in Fig. (3) is not very sensible to the precise choice for b ≡ 0(R vdW ), showing that for a narrow resonance the present model gives quantitative results for the deviation from universality.

We thank Y. Castin, F. Ferlaino, R. Grimm, F. Werner for fruitful discussions and W. Ketterle for providing us with experimental data. LPTMC is UMR 7600 of CNRS and its Cold Atoms group is associated with IFRAF.


Laetitia Chapel
email: laetitia.chapel@cemagref.fr

Guillaume Deffuant

Sophie Martin

Christian Mullon
email: christian.mullon@wanadoo.fr

Defining Yield Policies in a Viability Approach

Keywords: Viability theory, marine ecosystem, fisheries management, Support Vector Machines

proposed a dynamical model of biomass evolution in the Southern Benguela ecosystem, including five different groups (detritus, phytoplankton, zooplankton, pelagic fish and demersal fish). They studied this model in a viability perspective, trying to assess, for a given constant yield, whether each species biomass remains inside a given interval, taking into account the uncertainty on the interaction coefficients. Instead of studying the healthy states of this marine ecosystem with a constant yield, we focus here on the yield policies which keep the system viable. Using the mathematical concept of viability kernel, we examine how yield management might guarantee viable fisheries. One of the main practical difficulties up to now with the viability theory was the lack of methods to solve the problem in large dimensions. In this paper, we use a new method based on SVMs, which gives this theory a larger practical potential. Solving the viability problem provides all yield policies (if any) which guarantee a perennial system. We illustrate our main findings with numerical simulations.

Introduction

The viability theory [START_REF] Aubin | Viability theory[END_REF] aims at controlling dynamical systems with the goal to maintain them inside a given set of admissible states, called the viability constraint set. Such problems are frequent in ecology or economics, when systems die or badly deteriorate if they leave some regions of the state space. For instance Béné et al. (2001) studied the management of a renewable resource as a viability problem. They pointed out irreversible overexploitation related to the resource extinction. [START_REF] Bonneuil | Making ecosystem models viable[END_REF] studied the conditions the prey-predator dynamics must satisfy to avoid extinction of one or the other species as a viability problem. [START_REF] Cury | Viability theory for an ecosystem approach to fisheries[END_REF] consider viability theory to advise fisheries. [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF] proposed a dynamical model of biomass evolution of the Southern Benguela ecosystem, involving five different groups (detritus, phytoplankton, zooplankton, pelagic fish and demersal fish). They studied this model in a viability perspective [START_REF] Aubin | Viability theory[END_REF], trying to assess, for a given constant yield, whether each species biomass remains inside a given interval, taking into account the uncertainty on the interaction coefficients. The aim was to identify constant yield values that allow persistence of the ecosystem. We extend the problem and we focus here on the yield policies which keep the system viable, instead of considering a constant yield.

Using the mathematical concept of viability kernel, we examine how yield management might guarantee viable fisheries. The viability kernel designates the set of all viable states, i.e. for which there exists a control policy maintaining them within the set of constraints. Outside the viability kernel, there is no evolution which prevents the system from collapsing. [START_REF] Aubin | Viability theory[END_REF] proved the viability theorems which enable to determine the viability kernel, without considering the combinatorial exploration of control actions series. These theorems also provide the control functions that maintain viability. This general approach shows several interesting specific aspects:

-It can take into account the uncertainties on the parameters which are generally high in ecosystem modelling. Here, we manage the uncertainties like in [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF]).

-The viability kernel can define a variety of different policies, which respect the viability constraints. Therefore, it offers more possibilities for negotiations and discussions among the concerned stakeholders than techniques which propose a single optimal policy.

The main limitation of the viability approach is its computational complexity. The existing algorithm for viability kernel approximation [START_REF] Saint-Pierre | Approximation of viability kernel[END_REF] supposes an exhaustive search in the control space at each time step. This makes the method impossible to use when the control space is of a 51 dimensions like in our problem. [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF] solved this problem with a method which is only adapted to linear equations of evolution. Here, we use a new method, based on support vector machines, which can be applied to non-linear models as well [START_REF] Deffuant | Approximating viability kernels with support vector machines[END_REF].

We present the viability model of the Southern Benguela ecosystem and we recall the main concepts of the viability theory. Then, we describe our main numerical results. We show the shape of the found viability kernel, and the corresponding possible yield policies. Finally, we discuss the results and draw some perspectives.

The viability model of the Southern Benguela ecosystem

Following a classical approach [START_REF] Walters | Structuring dynamic models of exploited ecosystems from trophic mass-balance assessments[END_REF], we suppose that the variation of the biomass of species i due to its predation by other species j depends linearly on the recipient and donor biomasses (B j and B i ), with respective coefficients r ji and d ji . The biomass lost by species i due to the predation by the other species is expressed by equation (1):

∑ + - = → j i ji j ji i B d B r dt i dB ) ( ) ( . (1)
The variation of the donor biomass B i due to this interaction takes into account the assimilation of the biomass of other species j, multiplied by a growth efficiency coefficient (denoted below by g i ). Therefore, the biomass gained by species i, because of its consumption of other species, is expressed by:

∑ + = ← j j ij i ij i i B d B r g dt i dB ) ( ) ( . ( 2
)
For the detritus, the variation of the biomass follows the same principle, but it also integrates the non-assimilated biomass of the other species, except phytoplankton, which is added to the detritus biomass B 1 (multiplied by its growth efficiency g 1 ):

∑ ∑ > + - = - ٢ ١ ١ ) )( ١ ( ) ( j k k kl l lk l B d B r g g dt d assimilate non dB . ( 3
)
The model of the Southern Benguela ecosystem considers trophic interactions (predation, consumption and catch) among 5 components: detritus (i = 1), phytoplankton (i = 2), zooplankton (i = 3), pelagic fish (i = 4), demersal fish (i = 5). In total, the biomass evolution can be written as follows:

١ ١ ١ ١ ١ ) ( ) ١ ( ) ١ ( Y dt d assimilate non dB dt dB dt dB dt dB - - + → - ← = i i i i Y dt i dB dt i dB dt dB - → - ← = ) ( ) ( . (4)
where g i is the growth efficiency of species i, Y i is the yield of species i. Figure 1 shows the structure of the ecosystem. [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF] take into account the uncertainty on parameters r ij and d ij , which is

expressed by: [ ] ij ij ij ij ij r r r r r δ δ + - ∈ , , [ ] ij ij ij ij ij d d d d d δ δ + - ∈ , (5)
They consider this model in a viability perspective, in order to study the persistence of the ecosystem and to define the impact of the fisheries. Given a constant yield, they define scenarios which result is a "healthy" system.

Extending the work of [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF], we incorporate the fisheries in this study as a control variable of the system, in order to find the yield policies which allow keeping the system viable. To guarantee a perennial system, the viability constraints are defined by:

٠ i i i m B M ≤ ≤ ≤ , max min ٠ y Y y i ≤ ≤ ≤ , [ ] y y Y i δ δ + - ∈ ′ , , i = 4,5, (6)
where m i is the minimum level for the resource, M i the maximal biomass that can be contained in the ecosystem, y min is the minimum level for yield for demersal and pelagic fish, and y max the maximum level. The parameter y δ limits the evolution of the fisheries between two time steps. We suppose that the levels of yields of pelagic fish and demersal fish are the same.

These constraints, which attain critical values of a "healthy" system allow one to link yield objectives with the principle of ecosystem persistence.

The viability analysis control problem and viability kernel approximation

In the viability problem, the controls are the yields on pelagic fish (Y 4 ), demersal fish (Y 5 ), and the uncertainty on coefficients r ij and d ij . This means that for any state of the system located in the viability kernel, there exist values of these parameters for which the system remains in the viability kernel at the next time step. Adding the constraints on the derivatives of Y 4 and Y 5

implies to add two dimensions to the state space, which would then be 7. This reaches the current computational limits, therefore, we supposed that Y 4 = Y 5 = Y. This hypothesis is of course not realistic, but we thought it would nevertheless be an interesting first step.

The viability control problem is to determine a control function: Saint [START_REF] Saint-Pierre | Approximation of viability kernel[END_REF] proposed an algorithm to approximate the viability kernel from the problem defined on a grid but the result is a set of points that is viable and it requires an exhaustive search in the control space, which is not possible in our case because the control is in the dimension 51.

) ( ' ), ( ' ), ( ), ( ٥ ٤ t Y t Y t d t r t ij ij → with i,j = 1,2,3,4,5 (7)
To approximate the viability kernel of the Southern Benguela ecosystem, we use a new algorithm [START_REF] Deffuant | Approximating viability kernels with support vector machines[END_REF] (see Appendix 1) which is built on previous work from Saint-Pierre (1994), using a discrete approximation of the viability constraint set K by a grid.

Its main characteristic is to use an explicit analytical expression of the viability kernel approximation, in order to make it possible to use standard optimization methods to compute the control. This analytical expression is provided by a classification procedure, the support vector machines (SVMs) (Vapnik, 1998, Cristianini andShow-Taylor, 2000). This algorithm is interesting in the case we study, because the analytical expression of the viability kernel allows to use optimization techniques in order to find the best evolution in high dimensional control spaces.

Numerical simulations

The donor and recipient control coefficients are derived from a mass-balanced Ecopath model for the ecosystem [START_REF] Shannon | Trophic flows in the southern Benguela during the 1980s and 1990s[END_REF]. We use the evaluation of the parameters provided in [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF]. Table 1 gives the values of the viability constraint set and we put The following figures present some results for given values of biomasses of each species. The boundaries of the axes are the constraints defined on the species represented. The approximation of the viability kernel is represented in grey. Inside the viability kernel, there is at least one viable path which allows keeping a healthy system and outside, there is no evolution which prevents the system from collapsing. We focus here on values of detritus biomass = 2000 tons/km² because this ensures the existence of a viability kernel for almost all the values of the others compartments. For a level of detritus biomass = 1600 tons/km², given values of zooplankton and phytoplankton are necessary to guarantee a viable path. For lower detritus biomass, there is no viable path: a threshold of detritus biomass is necessary for ensuring a perennial system.

In the algorithm used to approximate the viability kernel [START_REF] Deffuant | Approximating viability kernels with support vector machines[END_REF], we used a grid with 6 points per dimension (46000 points in total) and 1642 support vectors are necessary to define the boundary of the kernel.

We focus on the effects of fisheries on demersal and pelagic fish. The levels of pelagic fish, demersal fish and yield have an influence on the boundary of the viability kernel:

Effects of fisheries on demersal fish

-For low values of pelagic fish biomass, the demersal fish biomass must not be too high and consequently intensive fishery must be avoid (see the circle at the top left of Figure 2, Pelagic 5); -In the same way, when the biomass of pelagic fish is high, the value of demersal fish biomass must not be too low to guarantee a perennial system and some low levels of catch must be avoided (see the circle at the bottom of Figure 2, Pelagic 60);

-For mean values of pelagic fish biomass, there is no restriction about the fisheries.

Figure 3 presents a 2D slice of the viability kernel, when detritus biomass = 2000 tons/km², phytoplankton = 400 tons/km² and zooplankton = 130 tons/km². We note that the viability kernel is smaller: a high level of pelagic fish represents a non-viable situation. Again, some high and low levels of fisheries must be avoided. In general, when the value of zooplankton is higher, the viability kernel is smaller and there is no viable path starting from pelagic fish biomass = 60 tons/km².

Effects of fisheries on pelagic fish

We explore now the impact of fisheries on pelagic fish, keeping the same values for others species.

Figure 4 presents the viability kernel where detritus biomass = 2000 tons/km², phytoplankton = 90 tons/km², zooplankton = 100 tons/km² and for demersal fish = 5, 15, 30 tons/km².

We notice that fisheries affect the boundary of the viability kernel only when the demersal fish biomass is too low: the more the pelagic biomass is, the more the catch can be important.

However, whatever the level of demersal fish, the level of fisheries must be controlled to guarantee a healthy system. For mean values of demersal fish, the system is not viable for low and high values of pelagic fish. For high values of demersal fish, the pelagic biomass must not be too low to guarantee the persistence of the ecosystem.

8

For high values of zooplankton (see Figure 5), the viability kernel is smaller: some values of demersal fish and fisheries are necessary to ensure a viable path:

-For low values of demersal fish, the level of fisheries must be carefully set; lower and higher values of catch represent non-viable situation;

-For mean value of demersal fish, the system is not viable for high values of pelagic fish;

-Fisheries have an influence for high biomass of demersal and pelagic fish: a minimum level of yields is necessary to ensure ecosystem persistence (area surrounded in Figure 5).

Main results

Our study illustrates the potential utility of the viability kernel to help the definition of viable fishery policies: given values of the biomass of the five species, the viability kernel provides the levels of catch to avoid. In addition, the viability kernel defines some conditions in which the fisheries can be increased without compromising the viability. We notice that the maximum thresholds for fisheries used by [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF] can also be increased.

Discussion and conclusion

Solving the viability problem provides all yield policies (if any) which guarantee a perennial system. This study shows that it is possible and interesting to integrate fisheries as a control parameter of a viability problem. We made strong simplifications: we supposed the same yield for the two species, and we should obviously take other parameters into account, like social and economics issues [START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF]. Nevertheless, we think that this work illustrates the potential of the viability approach to help the definition of fishery policies.

One of the main practical difficulties up to now with the viability theory was the lack of methods to solve the problem in a large number of dimensions. The use of learning procedures such as SVMs gives this theory a larger practical potential. However, to deal with a problem of six dimensions with the current algorithm can only be done with a very rough precision and several improvements are necessary to get more reliable and accurate results.

Moreover, it will be interesting to define yield strategies which allow the system to come from a non-viable state back to a viable state in minimum time, or minimizing some cost. This relates to the definition of the resilience proposed in [START_REF] Martin | The cost of restoration as a way of defining resilience: a viability approach applied to a model of lake eutrophication[END_REF]). ( 10)

h h K K ∈ ∃ ∈ ∀ x x , such as ) (h h β < -x x
At each step n , we define a discrete set where i x is a real vector and

{ } ١ , ١ - ∈ i y
, SVM define a function f

which separates examples of each labels:

(11) b k y f i i n i i + = ∑ = ) , ( ) ( ١ x x x α with ٠ ≥ i α and         - - = ٢ ٢ ² exp ) , ( σ x x x x i i k .
In [START_REF] Deffuant | Approximating viability kernels with support vector machines[END_REF], we show that it is possible to find an optimal control vector u*, which defines the position the most inside the current approximation of the kernel among all possibilities in G(x) (we use a gradient algorithm).

The steps of the algorithm are the following:

• Initialize the sets

h h K K = ٠ and K K L h = ) ( ٠ .
• Iterate:

o Define the discrete set ١ + n h K from n h K and n f as follows: { n h h n h K K ∈ = + x ١ such that ( ) ١ *) , ( - ≥ + dt f h h n u x x ϕ and ( ) } K dt h h ∈ + *) , ( u x x ϕ o If n h n h K K ≠ + ١
then run the SVM on the learning sample obtained with the points h x of the grid h K , associated with the labels ١

+ if ١ + ∈ n h h K x
, and with labels

١ -otherwise. Let ١ + n f be the obtained classification function. ) ( ١ + n h K L is defined as follows: { K K L n h ∈ = + x ) ( ١ such that } ١ ) ( ١ + = + x n f o Else stop and return ) ( n h K L .
13

List of tables

Tab 1 -Estimation of the minimal and maximal biomasses (Bi) for the five species..............15 Pelagic: 5 Pelagic: 30 Pelagic: 60

List of figures

  keep the viability constraints (6) satisfied indefinitely. Solving this problem requires to determine the viability kernel, which is the set of states for which such a control function exists.

  (no catches at all), max y = 5 tons/km² (the minimal level of the biomass of pelagic and demersal fish, corresponding on the maximum constant value tested by[START_REF] Mullon | Viability model of trophic interactions in marine ecosystems[END_REF]), y δ = 0.5 (which represents a variation of 10% of the maximal yield). The yield for others species has been set to 0, except for detritus (

Figure 2

 2 Figure 2 presents a 2D slices of the viability kernel where detritus biomass = 2000 tons/km², phytoplankton = 100 tons/km², zooplankton = 90 tons/km² and for different values of pelagic fish biomass. Horizontal axis represents demersal fish and vertical one the fisheries.

  of the discrete set and which constitutes the current approximation of the viability kernel. The boundary of this set is define thanks to a particularly procedure, the support vector machines (SVM), which is a method for data classification. Given a

Figure 1 -FiguresFigure 1 -

 11 Figure 1 -Components and structure of the Southern Benguela ecosystem. Arrows represent the flux between compartments (from Mullon et al., 2004).....................................................16 Figure 2 -Approximation of viability kernel. The horizontal axis represents demersal fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 90 tons/km² and phytoplankton = 100 tons/km²..................................................................................................16 Figure 3 -Approximation of viability kernel. The horizontal axis represents demersal fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 130 tons/km² and phytoplankton = 400 tons/km²..................................................................................................16 Figure 4 -Approximation of viability kernel. The horizontal axis represents pelagic fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 90 tons/km² and phytoplankton = 100 tons/km²..................................................................................................17 Figure 5 -Approximation of viability kernel. The horizontal axis represents pelagic fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 130 tons/km² and phytoplankton = 400 tons/km²..................................................................................................17

Figure 2 -

 2 Figure 2 -Approximation of viability kernel. The horizontal axis represents demersal fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 90 tons/km² and phytoplankton = 100 tons/km².

Figure 3 -Figure 4 -Figure 5 -

 345 Figure 3 -Approximation of viability kernel. The horizontal axis represents demersal fish, vertical axis fisheries, detritus = 2000 tons/km², zooplankton = 130 tons/km² and phytoplankton = 400 tons/km².

Acknowledgements

The authors are grateful to I. Alvarez, J.P. Aubin, N. Bonneuil and P. Saint-Pierre for useful discussions.

Appendix 1: Algorithm of SVM viability kernel approximation

We consider a given time interval dt and we define the set-value map

Considering the compact viability constraint set K, the viability kernel of K under G is the largest set included in K such that, for any x in Viab(K):

We define a grid h K as a finite set of K such that: