Dr Pierre-Sylvain Mirade
email: mirade@clermont.inra.fr

CFD Modelling of Indoor Atmosphere and Water Exchanges during the Cheese Ripening Process

I. INTRODUCTION

Approximately a third of the world's milk production is used in cheese manufacture, which is a way of preserving milk since it is highly perishable. Cheese is a highly nutritious food that offers a diversity of flavours and textures [START_REF] Ny Farkye | Cheese technology[END_REF]. In the 25 countries of the European Union, total cheese production in 2004 exceeded 8,200,000 tonnes, almost 60% of which was in Germany, France and Italy. Indeed, France is considered the quintessential "Cheese Country", with over 400 different varieties of ripened cheese and a consumption per capita reaching 25 kg a year, i.e. one of the most highest consumption levels of ripened cheese in the world.

In the successive operations involved in cheese manufacture, i.e. standardization of milk, enzymatic or acid coagulation, draining and ripening, the ripening process is crucial, since this is the step where texture, aroma and flavour (i.e. cheese quality) are developing. Cheese ripening is an outcome of several biochemical and metabolic processes including proteolysis, which is the most complex, glycolysis and lipolysis.

It is therefore crucial to control both airflow, i.e. air velocity, air change rate and renewal in new air, and climatic conditions (air temperature and relative humidity, gas concentration) inside the cheese ripening rooms, since this will determine both the efficiency and the homogeneity of cheese ripening and the water losses of the cheeses. However, it is difficult to achieve homogeneity in the distribution of climatic conditions at every single point of a ripening chamber. Consequently, industrial plants experience significant differences in the distribution of indoor atmospheric conditions which causes damage to the cheeses being ripened, as has often been underlined in the few studies published on this topic. For example, the problem of damage to crust formation observed inside an Emmenthal ripening room was attributed by the author to differences of over 10% in the relative humidity field [START_REF] Pajonk | La solution flux variables alternants pour homogénéiser l'ambiance[END_REF]. Air velocities ranging from less than 0.05 to 0.40 m.s -1 were measured inside the stacks of 10 cm-diameter cheese models filling an 84 m 3 ripening room, thus leading to a threefold increase in the heat and water transfer coefficients [START_REF] Ps Mirade | Caractérisation expérimentale de l'aéraulique d'un hâloir de fromagerie et des échanges air-produit[END_REF]. Significant heterogeneity in air velocity distribution, with marked differences in the cooling of the cheeses due to air velocities ranging from 0.1 to over 2.5 m.s -1 , was also highlighted in two cheese cooling stores [START_REF] Falconer | Cooling air distribution in a cheese store[END_REF]. The heterogeneity in climatic conditions that exists in ripening chambers means that cheese-makers have to regularly move the cheeses to achieve even water losses and a uniform appearance of the cheese surface. Nevertheless, only a few studies on this very real problem can be found in the literature because they are performed directly by industrial manufacturers or else by scientists or engineers in close collaboration with them, and therefore they remain confidential. Even books dealing with cheese ripening processes do not provide accurate information on the interrelationship between ventilation, indoor atmosphere and cheese quality; the authors only recommend a homogeneous atmosphere with low air circulation around the cheeses (air velocity not exceeding 0.1 m.s -1 ), while maintaining a "high enough" air change rate per hour in the plant [START_REF] Eck | [END_REF][START_REF] Mahaut | Initiation à la Technologie Fromagère[END_REF].

In fact, the role played by ventilation at whole-room level is complex and remains poorly quantified.

Air circulation on the one hand allows the evacuation of the heat and moisture produced by the cheeses inside the stacks and, on the other hand, determines both water losses and gas concentrations in the atmosphere closely surrounding the cheeses, which itself influences the cheese ripening process. For example, during the manufacturing of the Camembert cheese, the presence of ammonia and oxygen in the atmosphere of the ripening room makes it easier to reduce acidity on the cheese surface [START_REF] Vassal | L'amertume des fromages à pâte molle de type Camembert: rôle de la présure et de Penicillium caseicolum, moyens de la contrôler[END_REF] and for the Penicillia to grow [START_REF] Roger | Respiration of Penicillium camemberti during ripening and cold storage of semi-soft cheese[END_REF]. Furthermore, the presence of carbon dioxide in the indoor atmosphere is reported to increase the opening of the curd of hard cheeses by stimulating propionic fermentation [START_REF] Mahaut | Initiation à la Technologie Fromagère[END_REF].

Cheese ripening was originally performed in caves or in vaulted cellars often underground, where natural ventilation with low air velocity allowed a perfect control of the indoor atmospheric conditions throughout the ripening time, which can last from a few days to several months depending on the cheese-making technology. AOC Roquefort cheese is still produced under these conditions inside underground cellars fitted out under an immense mass of fallen rocks in which humid and fresh air naturally penetrates through long faults. As a result of the industrialization of the manufacture of cheese that has occurred over the last fifty years, modern ripening rooms are now being fitted with mechanical ventilation in combination with air temperature and relative humidity conditioning. Air distribution inside the ripening chambers is increasingly realized through blowing ducts made of textile materials in order to ensure low air velocities and therefore a correct treatment of the products. There are various types of textile blowing ducts: some are partially or fully porous, while others are fitted with holes or slits. Although their design obviously impacts on airflow patterns inside industrial plants, very few blowing duct design studies can be found in the literature, especially in relation to the cheese ripening process [START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF]. The only closely-related studies published deal with the effect of the location of blowing areas on the effectiveness of ventilation in relation to comfort in the buildings [START_REF] Jones | Computational Fluid Dynamics for building air flow prediction: current status and capabilities[END_REF] or on the efficiency of contaminant removal in clean rooms [START_REF] Rouaud | Numerical investigation on the efficiency of transient contaminant removal from a food-processing clean room using ventilation effectiveness concepts[END_REF]. Other interesting studies concern the assessment of airflow within large rooms [START_REF] Wang | Distributed dynamic modelling of a refrigerated room[END_REF][START_REF] Hoang | Analysis of air flow in a cold store by means of computational fluid dynamics[END_REF][START_REF] Lam | CFD analysis and energy simulation of a gymnasium[END_REF] or the internal climate in greenhouses [START_REF] Reichrath | Using CFD to model the internal climate of greenhouses: past, present and future[END_REF][START_REF] Shklyar | Numerical model of the three-dimensional isothermal flow patterns and mass fluxes in a pitched-roof greenhouse[END_REF]. In all these studies, the authors most of the time combined experimental investigation with CFD modelling to assess airflow and indoor atmosphere patterns. This chapter presents recent studies dealing with cheese ripening chambers ventilated through blowing ducts made of textile materials as a review of the modelling of airflow and indoor atmosphere (air temperature and relative humidity fields, gas concentrations) patterns, while evaluating the consequences on water exchanges of cheeses bathing in this atmosphere.

II. OPERATION OF CHEESE RIPENING ROOMS

The geometry and operation of modern cheese ripening rooms are relatively simple. The airflow conditioned in temperature and relative humidity is generally supplied continuously through one or several circular blower ducts made of textile materials, depending on the width of the plant, and located level with the ceiling. The type and design of the blower ducts change according to cheese variety and sometimes in relation to the ripening phase, such as at the end of ripening just before the packaging step. After being introduced through the blower ducts, the air circulates continuously around the cheeses placed into stacks whose height is generally about half the height of the ripening room, thereby exchanging heat and, above all, moisture. The air is then extracted in the lower part of the room level with the floor via the air conditioning system located either at one end of the plant or, more rarely, in the middle.

Figure 1 provides an illustration of the typical geometry of a modern cheese ripening room, representing a pilot ripening room built in our laboratory whose operation is truly representative of current standards in industry. All the results presented in this paper were obtained in this plant, and therefore a full description of its geometry and operation has been given. This pilot ripening room was 5.8 m long, 4.95 m wide and 2.95 m high, which gives an overall volume of about 84 m 3 . As this pilot room was to be used for the assessment of airflow patterns and not the heat and water transport stemming from the interaction between the cheeses and their surroundings, real cheeses were replaced by empty cans 10 cm in diameter and 4.4 cm high, i.e. inert objects presenting the same resistance against air circulation. Six rows (visually three, since they were placed two by two) of 7 stacks of 16 racks of 21 cans (i.e. a total of 14,112 cans) were then installed inside this room to obtain a filling pattern that was representative of current industrial practice. Free space, i.e. the height between two consecutive racks, was 10 cm. An air conditioning system composed of two fans and two batteries and installed in a space located above the ceiling of the pilot ripening room controlled the temperature and flow rate of the air blown into the room. Inside the pilot ripening room, the ventilation system was designed to mimic a common industrial configuration, i.e. an air conditioning system placed at the end of the room, extracting air in its lower part and blowing the conditioned air through a duct in its upper part.

This system was composed of a 340 mm blowing duct made of textile material and a 315 mm suction duct. In its standard configuration, the blowing duct running along the ceiling at half-width in the room was fitted on each side with three rows of several hundred holes 6 mm in diameter. Blowing ducts fitted with eight rows of more than one thousand holes 3 mm in diameter or with one row of several tens of holes 20 mm in diameter were also used to ventilate the pilot room. After being blown into the room, the air was extracted at 35 cm from the ground by means of a suction duct placed against a vertical wall at half-width in the room. The suction duct was connected to the space located above the ceiling of the pilot ripening room where the fans were installed. The full airflow rate blown into the room was 1600 m 3 .h -1 , i.e. an air change rate of 19 volumes per hour, which corresponds to normal industrial practice. The air velocity magnitudes were 11 m.s -1 at the output of the holes, regardless of the blowing duct considered.

III AIRFLOW PATTERNS IN CHEESE RIPENING ROOMS

Few experimental and numerical studies of the overall operation of modern cheese ripening rooms are available in the literature because most of them are confidential as a result of being performed directly by industrial manufacturers or jointly with scientists. Nevertheless, this section presents CFD results validated by experimental measurements and performed at the laboratory, highlighting heterogeneity in air velocity distribution inside the pilot ripening room whose geometry is depicted in Figure 1 [START_REF] Ps Mirade | Caractérisation expérimentale de l'aéraulique d'un hâloir de fromagerie et des échanges air-produit[END_REF][START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF][START_REF] Ps Mirade | Computational fluid dynamics prediction and validation of gas circulation in a cheese-ripening room[END_REF].

A. Description of the CFD Models Constructed

Based on the geometrical configuration of the pilot ripening room, the Fluent 6.0.20 code [START_REF]Anonymous. Fluent 6: User's Guide[END_REF] was used to construct numerical models based on an unstructured 3D mesh of about 1.2 million hexahedral and tetrahedral cells. The inside of the blowing and extraction ducts was considered as being outside the computational domain, thus air inflow corresponded to the output from the holes of the blowing duct and air outflow to the bottom of the extraction duct. About five thousand hexahedral cells were used for meshing each of the rows of the stacked racks of cheese models. Each hole of the blowing duct was meshed with at least 8 triangular cells, while the outlet area was meshed with 328 triangular cells. To make it possible to link the fine mesh of the holes to the coarser mesh built into the other parts of the blowing duct, hexahedral 'boundary layer'-type cells were placed in the vicinity of the holes.

In the inlet area corresponding to the holes, an air velocity of 11 m.s -1 determined from in situ measurements performed with a hot-film anemometer (model 8465, TSI, Minnesota, USA) and a turbulence rate arbitrarily valued at 10% were specified. A classical outflow-type boundary condition was During the numerical calculations, airflow was considered as steady, incompressible, isothermal and turbulent. To assess the sensitivity of the calculations to turbulence modelling, main flow turbulence was taken into account using the very popular standard k-ε model [START_REF] Be Launder | Mathematical Models of Turbulence[END_REF], the Renormalization group (RNG) k-ε model [START_REF] Yakhot | Renormalization group analysis of turbulence: I. Basic theory[END_REF] or the k-ω model [START_REF] Wilcox | Turbulence Modelling for CFD[END_REF] far from the walls, which were assumed to be smooth and where the standard wall function was applied. The Simple algorithm [START_REF] Sv Patankar | A calculation procedure for heat, mass and momentum transfer in three-dimensional parabolic flows[END_REF] was chosen for coupling pressure and velocity and introducing pressure into the continuity equation. The first-order or second-order upwind differencing scheme was also chosen in the computational models as a discretization scheme for the convection terms of each governing equation. Although first-order schemes are known to increase numerical discretization error and therefore to give less accurate results, especially when the flow is not aligned with the mesh, better convergence of calculation is obtained when using first-order versus second-order schemes [START_REF] Hirsch | Numerical Computation of Internal and External Flows: Computational Methods for Inviscid and Viscous Flows[END_REF].

Calculations were performed on an Athlon 1900XP+ PC with 1.5 Go of RAM and computation time was about 25 to 28 h depending on the model being solved.

In all the models, the material filling the pilot ripening room was taken into account as a porous medium.

Porous media are modelled by adding a momentum source term to the standard fluid flow equations. The source term is composed of two parts: a viscous loss term (the first term on the right-hand side of Eq. ( 1)) and an inertial loss term (the second term on the right-hand side of Eq. ( 1)). In the case of a simple homogeneous porous medium, source term S i is formulated in the CFD code following the equation:

z y x i where v v C v S i i i i , , 2
1 2 ≡ ⎟ ⎠ ⎞ ⎜ ⎝ ⎛ + - = ρ α µ (1)
To model turbulent flows through packed beds, Ergun [START_REF] Ergun | Fluid flow through packed columns[END_REF] established a semi-empirical correlation applicable over a wide range of Reynolds numbers. Comparing the Ergun equation with Eq. ( 1), the permeability α and inertial loss coefficient C 2 can be identified as:

2 3 2 ) 1 ( 150 ε ε α - = p D (2) 3 2 ) 1 ( 5 . 3 ε ε - = P D C ( 3
)
where ε is defined as the ratio of the volume of voids divided by the volume of the packed bed region.

Although a stack of cheeses is not exactly like a packed bed owing to distances between cheeses varying greatly, unlike the elements filling a packed bed, and although the cheeses could not be considered as spherical particles, the Ergun formulation was nevertheless applied to the stack. However, the D p and ε values were assumed to be variable according to the three spatial directions due to the anisotropy of the porous medium. Thus, in the purely empirical approach developed, D p was redefined as being a characteristic length of the elements filling the volume considered as a porous medium, and ε remained a void fraction but defined as the surface of voids divided by the total surface of the face of the porous volume perpendicular to the flow direction under study. In the stack of cheese models, the parameter D p corresponded to the diameter of one can in the vertical direction and to its height in both transverse and longitudinal directions. The void fraction ε could be easily calculated in relation to the vertical direction from the geometry of the stack. Hence, a value of 40% was determined [START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF]. On the other hand, for the transverse and longitudinal directions, the ε values were identified by visually comparing the measured fields to the air velocity fields calculated from many CFD calculations taking into account the filling of the ripening room as a porous medium with void fractions equal to 40% in the vertical direction and to 60%, 70%, 80%, 85%, 90% or 95% in the two other directions. Indeed, the calculation of the void P.S. Mirade fraction from the actual geometry of one stack proved impossible in these directions due to serious difficulties in accurately determining the total surface of cans that was perpendicular to the air circulation and consequently blocked it. Hence, a value of 90% was determined further to the visual comparison of the calculated and measured air velocity fields, meaning a very poor resistance of the cheese models against the air circulation in the transverse and longitudinal directions [START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF]. This poor resistance resulted from the large free spaces between the cans, particularly between the top of the cans of one rack and the bottom of next rack located above.

From a scientific point of view, this method of adjusting the coefficients in the porous medium is far from satisfactory. A methodology based on an anisotropic mega-porous medium coupled with the Darcy-Forchheimer model is still in progress in the laboratory; in this approach, viscous resistance, inertial resistance and porosity factors had to be adjusted according to the three spatial directions. The first results obtained showed that in a food industry setting such as a stack of cheeses, the viscous resistance factor was equal to zero whatever the spatial direction considered [START_REF] Ps Mirade | Airflow modelling by computational fluid dynamics in an industrial plant filled with food products[END_REF]. This amounts to neglecting the Darcy term in comparison with the Forchheimer term in the expression of the momentum source added to the standard fluid flow equations. In other words, it means that the mean static pressure gradient is linear with the squared velocity within the porous medium, and not proportional to the velocity as expressed by the Darcy formulation. This assumption is often made in bioclimatology when assessing interactions between airflows and plants in greenhouses or wind circulation in forests [START_REF] Bartzanas | Numerical simulation of the airflow and temperature distribution in a tunnel greenhouse equipped with insect-proof screen in the openings[END_REF].

B. Heterogeneity in Airflow Patterns

The air velocity fields calculated and measured within the pilot ripening room were represented as coloured velocity intensity maps ranging from 'white' colour areas in which air velocities were lower than 0.1 m.s -1 to 'black' colour areas where air velocities exceeded 0.4 m.s -1 . To make the flow pattern easier to assess inside the stacks where air velocities were very low, all magnitudes higher than 0.4 m.s -1

were plotted in black, including all those that reached 1 m.s -1 in the vicinity of the outlet holes of the blowing duct. Given the symmetry in the distribution of air velocity according to a vertical plane crossing the blowing duct, only the half-width of the room is plotted in figures 2-5.

Air Velocity Fields Calculated

Figure 2 shows the air velocity patterns simulated around and into the stacks according to a vertical section located at half-length in the pilot ripening room and obtained from a CFD model accounting for turbulence using the standard k-ε model and the first-order differencing scheme, for the blowing duct with holes 6 mm in diameter [START_REF] Ps Mirade | Computational fluid dynamics prediction and validation of gas circulation in a cheese-ripening room[END_REF].

This figure shows that the air blown through the holes flows along the ceiling and the lateral wall at a velocity higher than 0.4 m.s -1 (the values even reached 1 m.s -1 in proximity to the blowing duct), before being separated into two bodies when reaching the top of the side stacks. From here, the first body of air continues to flow down along the wall before entering the stacks, while the second body of the airflow appears to travel towards the blowing duct, giving rise to a swirl above the side stacks with air velocities equal to 0.2 m.s -1 . The swirl leads to the formation of a poorly ventilated area above the stack located at half-width in the room and underneath the blowing duct, and in which air velocities do not exceed 0.1 m.s -1 . In the passages between the rows of stacks, the numerical model calculated air velocities ranging from 0.1 to 0.2 m.s -1 , with a slight increase in velocity when approaching the ground.

Inside the side stacks of the cheese models, a marked gradient in relation to height appears as the air velocities peak at 0.3 m.s -1 in the lower part of the stack while the upper part is clearly poorly ventilated, except in the last rack at the top of the stack owing to the swirl located above. Furthermore, ventilation within the stack located in the middle of the room was clearly poor whatever the height considered (Figure 2). Thorough analysis of the air velocities calculated at that site revealed that more than half these values were below 0.05 m.s -1 . The heterogeneity in air distribution highlighted by these calculations certainly has a strong impact on cheese ripening factors such as heat and water exchanges.

The other calculated values confirm the airflow pattern depicted in Figure 2, i.e. poor ventilation in the stacks located in the middle of the chamber and a marked gradient in the side stacks, but with a slight variation in magnitude as a result of a three-dimensional effect due to the presence of the suction duct at one end of the pilot ripening room.

Validation of the Numerical Results

To rapidly measure the air velocities above the stacks of cheese models and between the rows of the stacked racks, a fast method set up at the laboratory [START_REF] Ps Mirade | A new experimental method for measuring and visualising air flow in large food plants[END_REF] was applied using a specially built system (Figure 1) to support and automatically move a measurement system at a slow and constant velocity of 1.5 cm.s -1 . The measurement system comprised a multi-directional hot-film anemometer (Model 8465, TSI, St Paul, Minnesota, USA) connected to a data logger (Model DT600, DataTaker, Rowville, Australia). Inside the stacks where automatic and constant movement of the probes was impossible, recordings of air velocity had to be averaged over at least 40 s at each measurement point to obtain a constant value for mean velocity. The sensors (Model 8475, TSI) were moved into the free space of 50 mm between two stacked racks and positioned at the exact measurement point via a telescopic antenna. To easily connect the measurements performed with the two methods, a regular experimental mesh was set up in the three spatial directions, with one measurement point each 19 cm in width, 21 cm in length and 20 cm in height.

From a qualitative point of view, Figure 3 that shows the air velocities measured at half-length of the room confirms the results of the CFD model, i.e. higher air velocities along the ceiling and the side wall, the presence of a large swirl above the stacks, a marked gradient in air velocity distribution according to height in the side stacks, and poor ventilation in the stacks located in the middle.

Comparison between Figure 3 and Figure 2 shows that the CFD model quite correctly simulates the airflow patterns inside the room and the stacks of cheese models, but with some discrepancies in the prediction of airflow magnitudes. Indeed, the numerical model underestimates the air ventilation level in areas with strong gradients, such as in the swirl or at the side stacks. Three possible explanations for this limited accuracy may be put forward: (i) the integration into the model of the standard k-ε model and of the first-order upwind differencing scheme, (ii) the use of porous media which only roughly describe the preferential paths for airflow over the cheese models on the racks, or (iii) the potential impossibility of obtaining truly independent results due to an insufficiently fine mesh, despite the steadily increasing power of computers and despite the 1.2 million cells used [START_REF] Borth | Influence of mesh refinement on the numerical prediction of turbulent air flow in rooms[END_REF]. Dependence on mesh size was impossible to quantify because no additional cells could be added to the model due to the limited amount of memory size on the computer used.

From a quantitative point of view, when considering the 6 mm blowing duct, Table 1 indicates that the discrepancy between simulation and measurement calculated from 4,200 points of comparison came to about 0.12 m.s -1 in the half-volume of the room, 0.07 m.s -1 in the right-hand-side stacks, and 0.03 m.s -1 in the stacks located in the middle, when using the standard k-ε model for modelling turbulence and the first-order upwind differencing method as discretization scheme for the convection terms in the equations. As can be seen in Table 1, changing the turbulence model and working with second-order schemes does not improve the accuracy of the prediction of air velocity in the room or in the stacks, possibly due to the mesh created which did not have enough cells to give more accurate results.

P.S. Mirade

Regardless of the models chosen, the discrepancy between calculation and measurement still ranged from 0.12 to 0.13 m.s -1 in the half-volume of the room due to the poor performance of modelling in the swirl area, from 0.05 to 0.07 m.s -1 in the side stacks and around 0.03 m.s -1 in the other stacks. On the other hand, computation time increased dramatically and convergence/divergence of computation did not occur when using a second-order discretization scheme or the RNG k-ε turbulence model; the residuals remained high and stable around a value lower than the value at the first iteration, but slightly higher than the convergence criterion. Hence, the standard k-ε model and the first-order upwind differencing scheme proved to be robust, time-saving and, given the mesh built, as accurate as the other models and schemes tested. The mean discrepancy in the air velocity prediction was 0.12 m.s -1 in the half-volume of the room and 0.05 m.s -1 inside the stacks.

C. Effect of Design of Blowing Duct on Ventilation Homogeneity

In the same pilot ripening room, the influence of the design of the blowing ducts made from textile materials on both ventilation homogeneity and ventilation level was experimentally and numerically tested through three blowing ducts with holes of different diameters: 3, 6 and 20 mm, at constant total air inflow rate and velocity magnitudes at the output from the holes [START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF].

CFD Results

Figure 4 shows the air velocity field calculated around and into the stacks of cheese models according to a vertical section located at half-length in the pilot ripening room, for the blowing ducts with holes 3 mm in diameter (Figure 4a) and 20 mm in diameter (Figure 4b), respectively. CFD results concerning the 6 mm blowing duct are depicted in Figure 2. Although the airflow rate blown into the room and the velocity magnitudes at the output from the holes were the same in all three configurations, analysis of Figures 4 and 2 reveals different ventilation levels and ranging heterogeneity around the cheese models. However, regardless of the design of the blowing duct, qualitatively speaking, air circulates in the same way once blown through the holes, giving rise to velocities higher than 0.4 m.s -1 along the ceiling and the lateral wall, a poor ventilation underneath the blowing duct and in the stack located in the middle of the plant, a large swirl above the stacks, and a marked gradient in air velocity distribution with height in the side stacks. On the other hand, the gradient of air velocities varies with the diameter of the blowing duct holes: smaller holes gave lower air velocities and lower heterogeneity around the cheese models. At a duct hole diameter of 3 mm, the air velocities ranged from less than 0.1 m.s -1 to 0.3 m.s -1 (Figure 4a), whereas with a hole diameter of 6 mm they ranged from less than 0.1 m.s -1 to nearly 0.4 m.s -1 (Figure 2), and with a hole diameter of 20 mm, peak air velocity even reached 0.6 m.s -1 (Figure 4b). This variation probably results from the different number of holes between the three configurations. Indeed, reducing the hole diameter implies a strong increase in the number of holes in order to maintain the same airflow rate blown into the chamber, since there was no variation in air velocity magnitudes at the hole outlets. Consequently, the number of the air jets around the blowing duct increases, therefore increasing friction between them, causing a loss of energy in the airflow that reaches the top of the side stacks and penetrates inside them with a velocity that decreases further as the hole diameter tends towards 3 mm. Increasing the hole diameter contributes to an increase in the inlet jet momentum, thus allowing the jet to remain stable for longer and enabling high velocities to predominate until the airflow reaches the top of the side stacks.

Validation of the CFD Results

Figures 5a, 3 and 5b show the air velocity values measured at half-length of the pilot ripening room when ventilated with the 3 mm blowing duct, the 6 mm blowing duct and the 20 mm blowing duct, respectively. Comparison between numerical and experimental air velocity values reveals a quite correct qualitative prediction of the airflow patterns inside the room and the stacks of cheese models, but with an underestimation of the air ventilation level in areas with strong gradients, such as in the vicinity of the blowing duct where the friction phenomena between the different jets take place and in the swirl located above the stacks or in the lower part of the side stack. In addition, Figures 5a andb confirm the numerical results indicating that an increase in blowing duct hole diameter from 3 mm to 20 mm leads to an increase in ventilation level and air velocity gradient in the side stack.

Table 1, which also details absolute values for the mean discrepancy between calculation and measurement for the 3 mm-and 20 mm-diameter blowing ducts, again underlines that coupling the k-ε turbulence model and the first-order upwind differencing scheme constitutes a good compromise solution for saving time while yielding accurate prediction of air velocities. Although the parameters of the porous medium modelling the filling of the pilot ripening room were previously adjusted based on the results obtained with the 6 mm blowing duct, the accuracy of the numerical models for the two other blowing ducts was absolutely satisfactory. The greatest discrepancy observed between prediction and measurement, which was 0.15 m.s -1 , corresponded to the 20 mm blowing duct, i.e. the hole diameter leading to the strongest gradient in air velocities within the side stack and in the swirl area visible above the same stack. Using the k-ε or k-ω model and especially the first-order differencing scheme to predict velocity magnitudes in this configuration performed poorly compared to the two other cases in which velocity gradients were lower, thereby increasing mean discrepancy with the measured values. Better accuracy would almost certainly have been achieved by working with the RNG k-ε model or/and second-order discretization schemes, provided that convergence of calculations had been fully satisfied, but this condition was never fulfilled (Table 1).

Although no variation occurred in the total airflow rate blown in the ripening room nor in the velocity magnitudes at the output from the holes, both experimental and numerical results revealed that changing the hole diameter led to differences in ventilation levels in the chamber and around the cheeses, almost certainly due to a variation in the friction phenomena occurring between the different air jets in the vicinity of the blowing duct. Hence, careful attention should be paid to the choice and use of textile blowing ducts for ventilating cheese ripening rooms and, more generally, industrial food plants.

IV INDOOR ATMOSPHERE IN CHEESE RIPENING ROOMS

Besides air velocity patterns, indoor atmosphere variables including temperature, relative humidity and gas (CO 2 , O 2 , NH 3 ) concentrations in ripening rooms are also crucial parameters during the cheese ripening process. Therefore, CFD models initially constructed to calculate air velocities in the pilot ripening room were then adapted to predict air temperature and relative humidity fields as well as CO 2 concentrations [START_REF] Ps Mirade | Computational fluid dynamics prediction and validation of gas circulation in a cheese-ripening room[END_REF]29].

A. Air Temperature and Relative Humidity Fields

Once the air velocity field in the pilot ripening room was determined, the numerical model was modified to solve the energy conservation and species transport equations in order to calculate air temperature and relative humidity fields [29]. An air temperature of 286 K and a relative humidity of 86% were specified in the outlet from the blowing holes. Constant heat and water vapour source terms were directly introduced into the porous media to account for the interaction between air and cheeses.

Their values were determined from industrial data, i.e. a heat flux of 10 W.m -3 and a water vapour flux of 3 kg.s -1 .m -3 corresponding to a daily water loss of 1% for a cheese weighing 400 g.

The temperature field calculated at half-length of the pilot ripening room is shown in Figure 6 covering the full width of the plant. Figure 6 indicates that apart from one relatively large area located underneath the blowing duct in which air temperature exceeded 13.4°C, air temperature everywhere else underwent low variations of only 0.2°C, even including the inside of the sides stacks in which a strong gradient in air velocity distribution was highlighted (Figures 23). Figure 6 also illustrates that any effect of this gradient was noticeable in the air temperature patterns. On the other hand, the area with high air temperature, located approximately at half-width of the room and in the upper part of the middle stack, could be unambiguously correlated to the poorly ventilated area visible in Figures 2 and3 at that site. It therefore appears that in this part of the middle stack and just underneath the blowing duct, the ventilation was not efficient enough to evacuate the heat generated inside the stack by the cheeses, giving rise to an accumulation of heat and thus an increase in air temperature.

Figure 7 depicting the relative air humidity field calculated at half-length of the pilot ripening room reveals that air moisture equal to 86% immediately at the output of the blower duct holes rapidly reaches 90-92% near the walls of the plant following a mixing of different air streams. As in the case of air temperature, the poor ventilation in the upper part of the middle stack and underneath the blower duct almost certainly caused the increase in relative air humidity that locally peaked at 97%. This increase revealed a strong accumulation of water vapour in air that was further increased by the fact that air temperature also increased. Indeed, in the range of high humidity, a 1°C increase in air temperature reduces the relative air humidity by 5-6%, obviously provided that the water mass in the air does not vary in the meantime. Calculations indicated a combined increase in air temperature and moisture, which means that the water mass in air increased in reality by 1 g.kg -1 of dry air (i.e. by about 12%) as a result of low air velocities at that location. Figure 7 also highlights a slight dissymmetry in the relative humidity field according to a plane crossing the blowing duct, almost certainly due to a slight dissymmetry in the air velocity field, although this dissymmetry was far from obvious when examining the air velocity patterns in the full width of the room. In reality, the rows of stacked racks of cheese models were not accurately distributed in relation to the blowing duct; the median plane of the filling of the ripening chamber was out of line with the median plane of the blowing duct by about 15 cm.

B. Mean Age of Air

As previously described in Chapter 8 -Section IV, the notion of Mean Age of Air (MAA) can be used to assess ventilation effectiveness in industrial food plants where air is used to treat products. MAA distribution can easily be determined through the resolution of an additional user-defined scalar transport equation in a CFD model. Inside a room, local MAA values provide information on the average time it takes for air to travel from the inlet area to any point of the room, and thus on the "freshness" of the air [START_REF] Rouaud | Numerical investigation on the efficiency of transient contaminant removal from a food-processing clean room using ventilation effectiveness concepts[END_REF][START_REF] Fisk | Air change effectiveness and pollutant removal efficiency during adverse mixing conditions[END_REF][START_REF] Abanto | Airflow modelling in a computer room[END_REF].

By solving the specific scalar transport equation allowing MAA to be computed by means of a userdefined function incorporated into the Fluent code, Chanteloup [START_REF] Chanteloup | Implémentation de la notion d'âge moyen de l'air à l'aide d'une « User-Defined Function[END_REF] calculated MAA distribution in 3 dimensions in the pilot ripening room and compared the results obtained with the air temperature and relative humidity fields previously established [29]. Figure 8 shows a strong heterogeneity in MAA distributions calculated at half-length of the pilot ripening room, since MAA values ranged from less than 150 s near the ceiling and the lateral walls of the room to approximately 220 s at half-width and half-height of the room. The lower MAA values obviously corresponded to the most ventilated areas while the higher values were concentrated in a location where air velocities were clearly lower than 0.1 m.s -1 . Comparing Figure 8 with Figures 6 and7 again proved that MAA appeared to be a better and more sensitive parameter than air velocity in detecting insufficiently ventilated areas leading to a combined accumulation of heat and moisture that is potentially harmful to the cheese ripening process.

The MAA criterion was so pertinent that the dissymmetry in airflow patterns due to a dissymmetry in the filling of the room in relation to the blowing duct was clearly highlighted.

C. Gas Circulation

In order to model how a gas added to air (namely CO 2 ) circulated in the volume of the pilot ripening chamber, the CFD model was adapted with the aim of determining the optimal injection point for introducing an exogenous gas to accelerate the ripening process [START_REF] Ps Mirade | Computational fluid dynamics prediction and validation of gas circulation in a cheese-ripening room[END_REF].

Two possible configurations were investigated: first, using the blowing duct to inject the additional gas via the holes, which is in theory the simpler solution, and second, directly injecting the gas into the core of the stacks of cheese models by means of injection points represented in the numerical model by nine small cubes of 20 cm 3 . Each of the faces of these nine cubes injected CO 2 at a velocity of 11 m.s -1 to maintain the same quantity of gas introduced into the ripening chamber atmosphere when compared to the first method where injection was performed through the blowing duct holes. The nine injection volumes of cubic shape were located at half-width, half-height, and at quarter, half and three-quarters of the length of the three rows of stacks.

CFD Modelling of Gas Transport

In order to calculate CO 2 transport into the chamber, the following additional convection-diffusion equation needed to be solved:

2 2 2 2 2 . ) .( ) ( CO CO CO CO CO S R J Y v Y t + + ∇ - = ∇ + ∂ ∂ r r ρ ρ (4)
where Y CO2 is the local mass fraction of the CO 2 species predicted by the Fluent code [START_REF]Anonymous. Fluent 6: User's Guide[END_REF] through the solution of Eq.( 4), R CO2 is the net rate of production by chemical reaction, and was zero in this configuration since no chemical reaction took place, S CO2 is the rate of creation by addition from the disperse phase plus any user-defined sources, and

2 CO J r
is the diffusion flux of CO 2 species arising from the concentration gradients.

In turbulent flows, the Fluent code [START_REF]Anonymous. Fluent 6: User's Guide[END_REF] computes the diffusion flux in the following form:

2 2 2 . , CO t t m CO CO Y Sc D J ∇ ⎟ ⎟ ⎠ ⎞ ⎜ ⎜ ⎝ ⎛ + - = µ ρ r (5)
where D CO2,m is the diffusion coefficient for CO 2 species and Sc t is the turbulent Schmidt number with a default setting of 0.7.

In practice, the computation of Eqs.( 4) and ( 5) is performed after total convergence in computation of the continuity and momentum equations has occurred, i.e. once the velocity field has been determined. Many assumptions were made in the numerical model built, i.e. that no interaction occurred between the CO 2 species and the cheeses, natural convection was not taken into account given that the energy equation was not solved, and that the gas extracted by the suction duct was not recycled back into the chamber through the blowing duct.

Computation time reached 60 h on an Athlon 1900 XP+ PC with 1.5 Go of RAM, when the additional transport equation for modelling gas circulation that required an unsteady computing was solved.

Comparison of the Two Methods of Gas Injection

CFD modelling aimed to dynamically determine where and how an exogenous gas has to be introduced into the pilot ripening room so that a homogeneous distribution is reached as quickly as possible. effect due to the extraction of the air at just one end of the room, slight differences between the three sections appear on Figure 9b, especially for the section farthest from the extraction duct during the first nine minutes of injection. However, using the blowing duct to add an exogenous gas appeared to be an efficient solution for evening out gas distribution throughout the whole volume of the ripening room.

As regards the second method consisting in injecting the gas by means of nine cubic volumes directly placed inside the stacks, Figure 10, which is a top view giving the CO 2 concentration calculated in a horizontal plane crossing the nine injection cubic volumes, clearly shows that this method was inefficient in evening out gas distribution inside the room. After 15 min of injection, calculations indicated that CO 2 concentration was far from homogenous since it ranged from nearly 0% to 1% just around the injection volumes. Immediately outside these volumes, CO 2 concentrations inexorably dropped as the gas was carried away by the airflow patterns before being extracted by the exit duct.

Hence, mean gas concentration values were higher in the stacks located in the middle that were poorly ventilated than in the side stacks. The calculations also revealed that even after 1 h of injection (data not shown), CO 2 concentrations had not varied since the first 10 min, regardless of the location in the room, and thus remained heterogeneous. Unlike the first method, adding gas by means of injection volumes placed at different locations in the stacks led to a heterogeneous distribution in the unit, as determined by the numerical results obtained. More homogeneous values would almost certainly have been obtained if the number of injection volumes within the stacks had been increased, but this option remains inconceivable from an industrial point of view. Indeed, each injection point would require a duct connected to a gas generator, thereby hampering cheese-makers who regularly move the cheese stacks.

V CONSEQUENCES ON WATER EXCHANGES OF CHEESES

Both calculations and measurements showed differences in air velocities in the side stacks of cheese models inside the pilot ripening room, with magnitudes ranging from less than 0.05 m.s -1 to only 0.3 m.s -1 , 0.45 m.s -1 or more than 0.5 m.s -1 according to the hole diameter of the blowing ducts tested.

It therefore appeared relevant and useful to assess the effect of this heterogeneity on the water losses of real cheeses that would have been placed at a few specific points inside the pilot ripening room. The location of these points is indicated in Figures 6 and7 by the digits 1 to 6.

A. Heat and Water Transfer Coefficients

P.S. Mirade

In order to calculate heat and water exchanges between airflow and cheeses, the convective heat and water transfer coefficients, h and k, first have to be determined. These coefficients are related to both the geometrical characteristics (shape and dimensions) of the object considered and to the airflow characteristics (mean air velocity and turbulence intensity), but not to the composition of the object. A few years ago, the laboratory adapted a psychrometric method based on wind tunnel experiments in order to accurately measure the heat and water transfer coefficients for different shapes and dimensions of objects found in the food industry in relation to mean air velocity and turbulence intensity [START_REF] Kondjoyan | Determination of transfer coefficients by psychrometry[END_REF]. In an airflow where air velocity and turbulence intensity were controlled, the psychrometric method consisted in regularly recording water losses due to evaporation of a fully wetted plaster cast of the object, together with its surface temperature and the air temperature. The mean values of h and k are then calculated independently according to the two following equations:

s air s air s s air vap m T T T T T T L h - - - - Φ - = ) ( 4 4 σ ε (6) )) ( ) ( ( s vap ws dew vap m T P a T P k - Φ = (7)
In these equations, the mean water flux ( m Φ ) was determined from the measured water loss from the plaster cast. The air temperature (T air ) was measured by means of a K-type thermocouple placed at 5 mm from the surface of the object. The mean surface temperature of the plaster cast (T s ) was measured by means of three K-type thermocouples placed just under the surface. Their number and location were judiciously chosen based on the results obtained during a previous study [START_REF] Ghisalberti | A thermographic method to map the local transfer coefficient on the complete surface of a circular cylinder in an air flow[END_REF]. The water activity at the surface of the object (a ws ) was equal to 1 owing to the saturation of the plaster cast.

Figures 11 and12 show the experimental values of the mean convective heat transfer coefficient (Figure 11) and the mean water transfer coefficient (Figure 12), both determined by psychrometry in the laboratory wind tunnel on plaster casts of cylinders 100 mm in diameter and 40 mm high, in an air velocity range lower than 0.5 m.s -1 [START_REF] Ps Mirade | Caractérisation expérimentale de l'aéraulique d'un hâloir de fromagerie et des échanges air-produit[END_REF]. During the experiments, turbulence intensity was set at 12%, i.e. a value corresponding to the mean value measured in different types of industrial food plants [START_REF] Jd Daudin | Influence de l'indice de turbulence de l'écoulement sur les procédés de traitement thermique de solides par l'air[END_REF].

Although h and k were measured independently, their values obeyed the Lewis relation [START_REF] Lewis | The evaporation of a liquid into a gas[END_REF] and were proportional even for these low air velocities. Moreover, the experimental values of h and k were smoothed by an exponential law, which is the most appropriate mathematical function for representing the mixed convection that occurs at this low velocity range. Indeed, power laws that are widespread in heat transfer problems are normally valid only for forced convection, i.e. at higher air velocities [START_REF] Ozisik | Heat Transfer: a Basic Approach[END_REF]. In addition, an exponential law gives coefficient values that have a real physical meaning when the air velocity tends towards zero, i.e. when natural convection predominates. Figures 11 and12 indicate a threefold increase of the measured heat and water transfer coefficients when the air velocity around the plaster cast increases from less than 0.1 m.s -1 to 0.5 m.s -1 . This proves that the location of the cheeses in the ripening room is an important factor in explaining variation in the water losses.

Figure 13 confirms that location has an influence on heat and water exchanges as it shows the measured water losses of 18 plasters casts 100 mm in diameter and 40 mm high, fully saturated with water and placed at different locations within the stacks filling the ripening room. The location of the plaster casts was chosen following the airflow study [START_REF] Ps Mirade | Caractérisation expérimentale de l'aéraulique d'un hâloir de fromagerie et des échanges air-produit[END_REF].

B. Water Losses of Real Cheeses

Water losses of cheeses can also be calculated from the two following equations: P dm = k . S . (a ws . P vap (Teq) -P vap (Tdew) ) . 3600 . 1000 . 24 (8)

P.S. Mirade

h . (T eq -T air ) + σ . ε s . (T eq 4 -T air 4 ) + k . L vap . (a ws . P vap (Teq) -P vap (Tdew) ) = 0 (9)

A programme was created using Matlab software (The MathWorks, Natick, USA) to calculate water loss from cheese, i.e. P dm , from Eq. ( 8). In this equation, the value of the water transfer coefficient k is determined from the exponential law of Figure 12, given the mean air velocity around the cheese calculated by the CFD model. Other variables are imposed, including total surface S of the cheese (equal to 0.03 m 2 ), water activity at the cheese surface (0.98 or 1), air temperature T air (14°C) and relative humidity RH (95% or 98%), which go to define dew temperature (T dew ). The last variable of Eq. ( 8) to be evaluated, i.e. the mean temperature of the cheese surface (T eq ), results from equilibrium between heat transfers occurring by convection, radiation and evaporation; T eq is determined numerically by solving Eq. ( 9), in which the convective heat transfer coefficient h takes a value based on the experimental law of Figure 11, given the mean air velocity calculated by the CFD model. Concerning the radiation exchange formulation in Eq. ( 9), the temperature of all the surfaces surrounding the cheese placed in the stack was assumed to be equal to the air temperature, thus artificially increasing heat transfer by radiation and therefore causing an overestimation of the mean surface temperature. In fact, calculations showed that the radiation term, although overestimated, had little effect on the water losses from the cheese, since the change in the mean surface temperature ranged from 0.03°C to 0.29°C at the very most.

Based on a room air temperature of 14°C, Figure 14 gives values of water losses from cheeses initially weighing 400 g as a function of water activity at the cheese surface, relative air humidity and mean air velocity magnitudes, in a range corresponding to the heterogeneity highlighted in the three configurations of the pilot ripening room investigated [START_REF] Ps Mirade | Effect of design of blowing duct on ventilation homogeneity around cheeses in a ripening chamber[END_REF]. Regardless of the relative air humidity and water activity at the cheese surface, water losses logically increase as mean air velocity increases around the cheese, owing to the variation of the water and heat transfer coefficients with air velocity indicated in Figures 11 and12. Figure 14 illustrates how modifying the diameter of the holes of the blowing duct can strongly affect the disparity in water losses according to the location of the cheeses in the ripening chamber. For example, replacing the blowing duct with 20 mm-diameter holes by the one with 3 mm-diameter holes can double the water losses from cheeses placed in the lower part of the side stacks near to the ground, although no variation occurred in the full airflow rate.

Table 2 presents the daily water losses of six cheeses initially weighing 400 g located as indicated by digits 1 to 6 in Figures 6 and7, and calculated from Eq. ( 8) considering the air velocity, relative humidity and temperature fields simulated. The water activity at the surface of the cheeses was assumed to be either equal to 1 as in the beginning of the ripening process, or 0.97 as occurs during the ripening process.

Analysis of Table 2 reveals that cheese water loss could increase by 500 to 600% depending on whether it is located in the upper part of the middle stack (location 1) where there is very low air velocity and high relative air humidity and temperature, or in the lower part of the side stacks (locations 4 and 6) where air velocity reached 0.4 m.s -1 and where air temperature and relative humidity were lower.

C. Simulation of the Coupling between Airflow and Water Losses

The previous results allow the water losses of cheeses and air temperature and relative humidity fields to be determined without numerically solving the coupling that exists between the airflow and the food products bathing in this air, which can cause differences, for example, in the fluxes of water evaporation between different points of the same stack of products. In the broad approach previously presented, heat and water fluxes were assumed as constant in the calculations of relative air humidity and temperature patterns in the pilot ripening room. However, in reality, during air treatment operations heat and moisture fluxes vary with time, and their time-point values depend on air and product surface characteristics. Indeed, heat and moisture fluxes are proportional to the difference in temperature and in water vapour pressure, respectively, between the air and the food product surface. Moreover, both

P.S. Mirade

fluxes are affected by airflow characteristics and product shape as a result of the influence of these parameters on heat and water transfer coefficients.

As general purpose CFD codes such as Fluent, Star-CD or CFX were designed for solving turbulent fluid flow problems coupled with heat and mass transfers, it is theoretically possible to calculate the spatial distribution of the air characteristics (velocity, temperature, relative humidity) and time variations of product temperature and weight. Nevertheless, there are several major factors that limit the capacities of these codes and the accuracy of the results obtained, including: (i) the lack of efficient modelling procedures to describe water diffusion inside solids, and (ii) the use of the wall function approach in modelling turbulent boundary layers near product surfaces, which leads to large errors in transfer coefficient predictions or which leads to very fine meshes and therefore to a cell number requiring an immense memory size well in excess of the computers currently used.

In order to accurately predict the air velocity, air temperature and relative humidity as well as the water loss of stacked food products, a specific user-defined function which can be compiled in the CFD code 'Fluent 6.1.22' was recently developed in the laboratory [START_REF] Ps Mirade | Simulation numérique de l'écoulement de l'air et des échanges de chaleur et d'eau au sein d'un empilement de produits alimentaires[END_REF]. Validation of the function implemented was achieved by comparing the results obtained for a few hot and moist cylinders placed in a cold and dry airflow with those obtained from the solving of analytical solutions using a programme written in Matlab programming language, giving the kinetics for temperature and water concentration at any point of a cylinder radius.

Figures 15a andb show the relative air humidity calculated after 37 hours (Figure 15a) and 45 days (Figure 15b) using the specific function developed around six cheeses 100 mm in diameter and 45 mm high, arranged in a row and placed in an airflow whose characteristics at the inlet area were a velocity of 0.01 m.s -1 , a temperature of 12°C and a relative humidity of 95.4%. The 3D computational model developed used symmetry conditions for all the external faces except for the inlet and outlet areas, and thus simulated water exchanges as they occurred in the core of a stack of cheeses. Airflow was considered as turbulent -with modelling using the standard k-ε model -, non-isothermal and unsteady, since heat and water transfer phenomena vary naturally with time. Salt transfers occurring inside the cheeses and which interfere with water transfers and water activity were neglected. A total of 350,000 cells were used for meshing the computational model, and the full convergence of all the discretized equations was reached after about 3.2 h on a Pentium Xeon 3.2 GHz PC with 4 Go of RAM.

Figures 15a andb show relative air humidity fields calculated after about 37 process hours (Figure 15a)

and 45 process days (Figure 15b) which vary both spatially and over time. As the air circulates around the cheeses as it moves from the inlet to the outlet, its relative humidity increases due to the accumulation of the moisture evaporating from the wet surface of the cheeses, thus generating the higher values around the cheeses located farthest from the air inlet. Around an individual cheese, relative air humidity peaks locally in the wake generated by the circulation of air around the cylinder, which is known to be a low velocity region. Comparison between Figure 15a and Figure 15b reveals that the user-defined function implemented in the CFD code correctly simulated the variations in the water transfer phenomenon over time, i.e. a decrease in the water evaporation on account of a decrease of the water content inside the cheeses and in water activity at the cheese surface.

Figure 16 presents the water loss kinetics of three cheeses (cheeses 1, 3 and 6 of the previous row)

predicted by the user-defined function. The numerical model logically reported differences in water losses according to the location of the cheeses, with the higher water loss for the first cheese in contact with "dry" air and the lower water loss for the cheese located near the outlet of the computational domain where the air is moistest. These differences decrease over time and tend towards equilibrium, since water content and activity at the surface of the first cheese decreases more rapidly compared with the others due to its higher water loss.

P.S. Mirade

The user-defined function developed and incorporated in the CFD code Fluent then makes it possible to compute temperature and water concentration fields and determine the water loss kinetics of the cheeses together with the temperature and relative humidity fields of the air flowing through the stacked food products. The laboratory is currently extending the approach developed to accurately calculate airflow and the water exchanges within a 3D stack filled with one hundred cheeses. A full quantitative experimental validation will be performed afterwards.

VI. CONCLUSION

Based on various recent studies performed in a pilot ripening room constructed in the laboratory, this chapter has clearly demonstrated that CFD techniques can be very useful tools for assessing the operation and indoor atmosphere of industrial cheese ripening rooms, even though the representation of the filled plant using cheese models was significantly simplified in the numerical models. In general, the calculated predictions of airflow patterns within the pilot cheese ripening chamber were in fairly close agreement with actual measurements. However, in some regions, i.e. in the side stacks or in the swirl located above these stacks, the measurements indicated stronger velocity gradients than predicted by the model. Moreover, careful attention must be paid to the choice of the blowing ducts in cheese ripening chambers, since blowing duct design has a very real influence on airflow patterns, and therefore on water losses from cheeses. Indeed, the model clearly highlighted differences in ventilation levels around cheese models following a simple change in the diameter of the holes of the blowing duct. Further calculations were performed to identify a solution for introducing an exogenous gas into the ripening room in a homogenous way. To achieve this, the calculations indicated that it was far more efficient to use the blowing duct than introduce the gas directly into the stacks.

However, in light of the results presented in this chapter, it is very difficult to give quantitative recommendations on ventilation and indoor atmosphere in ripening chambers, and equally difficult to evaluate the impact of the calculated gas distributions on the cheese ripening process. Further data is required in order to fully understand and quantify the interaction between the indoor climate of ripening rooms -air velocity, temperature and relative humidity, gas concentrations, renewal in fresh air, etc. -and the cheeses being ripened (weight losses, gas consumption and production at the cheese surface by microbial flora, microbial growth, etc.). Studies are in progress in the laboratory, together with research aimed at determining appropriate coefficients to represent the rows of cheese stacks in the CFD model through porous media coupled with the general Darcy-Forchheimer formulation.

The development of a system based on high-performance sensors coupled with information on airflow patterns would give ripening room operators greater flexibility in looking to improve the consistency and quality of the cheese-making process. To achieve this, research will need to focus on controlling cheese ripening by monitoring indoor atmosphere and on obtaining evenly-distributed ripening conditions around the cheeses. Further progress can be expected in the years to come due to the increasing calculation power of computers and the increasingly flexible CFD codes available. Thus, it will likely be possible to fully simulate how a ripening room operates by implementing unsteady numerical models capable of accurately predicting heat, mass and gas exchanges between indoor atmosphere and several hundred cheeses being ripened. the location of the cheeses for which calculated water losses are reported in Table 2. the location of the cheeses for which calculated water losses are reported in Table 2. The discretization scheme used in the numerical calculations was either a first-order or-second order upwind differencing scheme.

b Since convergence was not reached and divergence did not occur (the residuals remained stable around a value significantly lower than the value at the first iteration but slightly higher than the convergence criterion), discrepancy between calculation and measurement was evaluated when the residuals were minimal during solving.

Centrale d'acquisition

  applied at the bottom of the extraction duct. This kind of boundary condition is conventionally used to model flow exits where details of the flow velocity and pressure are unknown prior to solution of the flow problem.

Figure 9a displays the

  Figure9adisplays the kinetics of the mean concentration in CO 2 calculated for each of the three rows of

Figure 9b ,

 9b Figure9b, which gives the kinetics of the mean concentrations in CO 2 calculated according to three

Figure 1 .

 1 Figure 1. Illustration of the geometry of the pilot ripening room constructed in the laboratory: (a) an

Figure 2 .

 2 Figure 2. Vertical section located at half-length of the pilot ripening room showing the calculated air

Figure 3 .

 3 Figure 3. Vertical section located at half-length of the pilot ripening room showing the measured air

Figure 4 .

 4 Figure 4. Vertical section located at half-length of the pilot ripening room showing the calculated air

Figure 5 .

 5 Figure 5. Vertical section located at half-length of the pilot ripening room showing the measured air

Figure 6 .

 6 Figure 6. Vertical section located at half-length of the pilot ripening room showing the calculated air

Figure 7 .

 7 Figure 7. Vertical section located at half-length of the pilot ripening room showing the calculated

Figure 8 .

 8 Figure 8. Vertical section located at half-length of the pilot ripening room showing the distribution of

Figure 9 .

 9 Figure 9. Variation over time of the calculated mean concentrations of CO 2 previously injected into the

Figure 10 .

 10 Figure 10. Distribution of the calculated mean concentrations of CO 2 15 min after an injection performed directly into the stacks of cheese models by means of nine injection volumes of 20 cm 3 ,

Figure 11 .Figure 12 .

 1112 Figure 11. Influence of air velocity on the mean convective heat transfer coefficient values measured

Figure 13 .

 13 Figure 13. Influence of air velocity on the measured weight loss of plaster casts of cylinders 100 mm

Figure 14 .

 14 Figure 14. Influence of air velocity, relative air humidity and water activity at the cheese surface on the

Figure 15 .

 15 Figure 15. Distribution of predicted relative air humidity for a row of six cheeses using a specific user-

Figure 16 .

 16 Figure 16. Evolution over time of the water losses from three cheeses (cheeses 1, 3 and 6) placed in a

  Figure 1

1 )Figure 16

 116 Figure 14

  Miradecalculations that formed much of the results on which this chapter is based. The author is also grateful to the French Ministry of Research and 'Arilait Recherches' for their financial support to a major part of the studies presented.

	α		Permeability	m 2
	ε		Void fraction	-
	ε s		Emissivity of the object (0.91 for plaster)	-
	Φ	m	Mean water flux	kg.m -2 .s -1
	NOMENCLATURE µ Dynamic viscosity of the fluid	kg.m -1 .s -1
	ρ		Fluid density	kg.m -3
	a ws C 2 σ	Water activity at the surface of the plaster cast or the cheese Inertial resistance factor Stefan-Boltzmann constant (5.67.10 -8 )	-m -1 W.m -2 .K -4
	D p	Mean particle diameter	m
	h , h	Mean convective heat transfer coefficient	W.m -2 .K -1
	k , k	Mean water transfer coefficient	kg.m -2 .Pa -1 .s -1
	L vap	Latent heat of water vaporization	J.Kg -1
	P dm	Water loss from cheese	g.24h -1
	P vap	Saturating vapour pressure at the considered temperature	Pa
	RH	Relative air humidity	%
	S		Surface of the cheese	m 2
	T air	Air temperature around the plaster cast or the cheese	K
	T dew	Dew temperature calculated from air temperature and relative humidity	K
	s T , T eq	Mean surface temperature of the plaster cast or the cheese	K
	v i		Air velocity according to the direction i	m.s -1
	x, y, z	Three spatial directions	-
	Greek symbols

P.S.

Table 1 .

 1 Mean discrepancies in air velocity calculated as the absolute value of the difference between the calculated and measured magnitudes divided by the number of points of comparison, according to model used for solving the problem and blowing duct tested.

	Blowing duct tested	Models used for solving	Convergence	Mean discrepancy between measured and calculated
			reached	air velocity (m.s -1 )
	Hole diameter	Turbulence modelling /	Yes/No	Right-hand-side stacks	Stacks located in	Half-volume of the
		Discretization scheme a			the middle	ripening room
	6 mm	k-ε / 1 st order	Yes	0.07	0.03	0.12
	6 mm	k-ω / 1 st order	Yes	0.06	0.03	0.12
	6 mm	RNG k-ε / 1 st order	No	0.05 b	0.02 b	0.13 b
	6 mm	k-ε / 2 nd order	No	0.05 b	0.03 b	0.12 b
	6 mm	k-ω / 2 nd order	No	0.06 b	0.03 b	0.13 b
	20 mm	k-ε / 1 st order	Yes	0.08	0.06	0.15
	20 mm	k-ω / 1 st order	Yes	0.07	0.05	0.15
	3 mm	k-ε / 1 st order	Yes	0.07	0.03	0.12
	3 mm	k-ω / 1 st order	Yes	0.06	0.03	0.11
	a

ACKNOWLEDGEMENTS

The author sincerely thanks A. Kondjoyan, J.D. Daudin, A. Lasteyras, J.M. Auberger, S. Portanguen, L. Lasteyras, T. Rougier, D. Peyne, A. Carles, T. Laboue, R. Lardy, V. Chanteloup and C. Chevarin for their active participation and efficient contribution to the discussions, experiments and numerical

P.S. Mirade

Table 2. Calculated water losses of six cheeses initially weighing 400 g versus airflow characteristics (air velocity, temperature and relative humidity) determined by the CFD model, and versus two water activities (1 and 0.97) at the cheese surface. The location of the six cheeses inside the pilot ripening room is indicated on Figures 6 and7.

Air velocity (m.s -