V Sambou

B Lartigue
email: lartigue@cict.fr

F Monchoux

M H A S E Adj

Theoretical and experimental study of heat transfer through a vertical partitioned enclosure: application to the optimization of the thermal resistance

Keywords: partitioned enclosure, heat transfer, thermal resistance, optimization

Enclosures divided by multiple vertical diffusive partitions have high insulating qualities and may provide tangible benefits as a construction material. In this study, we have developed onedimensional analytical model of coupled heat transfer (conduction, convection, radiation) in such enclosures. This model is numerically and experimentally validated for an important number of configurations and size of alveolus. The variation of the thermal resistance versus the number of partitions has been experimentally validated. The model is then used to find the partitions number that maximizes the thermal resistance of a partitioned enclosure. Effects of thermal and geometrical parameters on the maximal thermal resistance of the partitioned enclosure have been also investigated. We have shown that the thermal resistance could be improved by decreasing the thermal conductivity of walls; decreasing the emissivity of the partitions faces; or using very thin partitions. The use of relatively thick exterior vertical walls marginally degrades the optimal thermal resistance while at the same time increases the thermal inertia. The combination of the different parameters gives a set of solutions when one desires to obtain precise characteristics for an insulating envelope. The presented model can help to determine the most suitable combination of parameters allowing to get the desired maximal thermal resistance.

Introduction

Residential and commercial buildings represent 46% of the total energy consumption in France [1], becoming the first final energy consuming. Corresponding greenhouse gas emission represents a quarter of the total emission. The Kyoto protocol aims reduce energy consumption and greenhouse gas emissions in residential as well as commercial buildings. As a consequence, in its thermal standards, the building code recommends industrial development of products with strengthened thermal insulation for opaque building walls [START_REF]Réglementation thermique 2000[END_REF]. Yet, high insulation which decreases heating energy during winter period leads to overheating during hot days. To overcome this phenomenon using thermal inertia is an efficient way [START_REF] Antonopoulos | On the dynamic thermal behaviour of indoor spaces[END_REF]. Construction materials with non-localized insulation aim to have this double property: they both are insulating materials and have a high heat capacity that increases thermal inertia of the building. Among these materials, alveolar materials have an important place.

Enclosures with vertical partitions may describe such construction material [START_REF] Lacarrière | Numerical study of heat transfer in a wall of vertically perforated bricks: influence of assembly method[END_REF]. They have been efficiently used to reduce conduction, convection and radiation heat transfer. In addition to alveolar materials, other engineering applications recover from such a system, like alveolar solar collectors, cryogenic storage, etc.

Laminar natural convection in differentially heated rectangular cavities has been studied extensively [START_REF] Catton | Natural convection in enclosures[END_REF][START_REF] Bejan | Convection Heat Transfer[END_REF]. More recently, several works outline the decrease of heat transfer using partitioned enclosures. Andersen and Bejan reported in [START_REF] Anderson | Heat transfer through single and double vertical walls in natural convection: theory and experiment[END_REF] the decrease of the overall heat transfer in single and double partitioned enclosure, both theoretically and experimentally. Other studies focus on the only natural convection heat transfer in enclosures divided by one [START_REF] Tong | Natural convection in partitioned air-filled rectangular enclosures[END_REF][START_REF] Nishimura | Natural convection heat transfer in enclosures with an off-center partition[END_REF] or multiple [START_REF] Nishimura | Natural convection heat transfer in enclosures with multiple vertical partitions[END_REF] nondiffusive partitions, showing that partitioned enclosures compare favorably with fully insulating the enclosure with a porous material [START_REF] Tong | Natural convection in partitioned air-filled rectangular enclosures[END_REF]. The coupled conduction-convection heat transfer is studied when considering thick partitions by Ho and Yih [START_REF] Ho | Conjugate natural convection heat transfer in an air-filled rectangular cavity[END_REF] in a single partitioned enclosure as well as by Kangni et al. [START_REF] Kangni | Natural convection and conduction in enclosures with multiple vertical partitions[END_REF] and Turkoglu and Yücel [START_REF] Turkoglu | Natural convection heat transfer in enclosures with conducting multiple partitions and side walls[END_REF] in multiple partitioned enclosures. These studies mainly emphasize the decrease of heat transfer by reducing the size of the alveolus. Radiation heat transfer is rarely taken into account although it is as important as the convection heat transfer. The complexity of three-dimensional coupled heat transfer as well as the complexity of the geometry of alveolar materials make difficult to optimize these products.

In this study, we have established an analytical model of conduction-convection-radiation heat transfer for an enclosure with multiple vertical diffusive partitions, representing more closely real applications. The performances of the system calculated with the analytical model have been compared with both the results of a numerical model and experimental results. The validated analytical model is then used to evaluate the sensitivity of geometrical parameters and thermal properties to the thermal resistance.

Heat transfer through partitioned enclosures

Geometry

The studied system is a two-dimensional partitioned enclosure of internal height H, internal width L and a supposed infinite depth (Fig. 1). The four external walls of the enclosure have a thickness e w .

N partitions of thickness e p are equally placed in the enclosure. The enclosure has therefore (N+1) alveolus of thickness e a calculated by: p a L N e e N 1

- = + (1)

One-dimensional analytical heat transfer model

Constant uniform temperatures T H and T C are imposed at external surfaces of vertical walls. Heat transfer is supposed one-dimensional in the enclosure. The enclosure is segmented into three visible parts on Fig. 1. The heat transfer is evaluated for a unit length in the Oy direction.

Heat transfer through horizontal walls (parts 1 and 3)

In parts 1 and 3 the thermal conduction resistance is calculated using:

h c w 1 3 w w w T T L 2 e R R Q k e - + = = = (2)
where k w is the thermal conductivity of the walls.

Heat transfer through the alveolus (part 2)

In the part 2, heat flow calculation should be computed taking into account either conduction or convection occurring in the alveolus, as well as radiation between the alveolus faces. If heat flux through the alveolus occurs by conduction, calculation is performed with:

i i i h c cd a a T T Q k H e - = (3)
If the heat flux is convective, it is then given by:

i i i h c cv a i a T T Q k Nu H e - = (4)
The Nusselt number in the alveolus i (Nu i ) is given by Bejan correlation [START_REF] Bejan | Convection Heat Transfer[END_REF]:

1 1 4 i i Nu 0.364 Ra A - = (5)
where Ra i is the Rayleigh number and A i the aspect ratio of the alveolus i.

( )

i i h c 3 i i i i T T Ra g H - = β α ν and a H A e = ( 6
)
where β i is the thermal expansion coefficient; α i the air thermal diffusivity and ν i the air kinematic viscosity calculated at

i i i h c m T T T 2 + = .
Calculation of the radiation heat transfer in the cavities supposes that air is a non-participant media, and that the radiative exchange of horizontal faces is negligible. Vertical faces are assumed to be perfectly diffuse with a constant emissivity coefficient. Then, radiation flux between cold and hot faces of the alveolus i is given by: ( )

i i i i rad rad h c Q h H T T = - (7)
The radiative coefficient i rad h is given by [START_REF] Saccadura | Initiation aux transferts thermiques[END_REF]: The global heat flux through the alveolus i is therefore given by: ( )

i i i i i i 3 rad m h c h c hc 1 h 4 T 1 1 1 F = σ -ε -ε + + ε ε (8)
i i i i g g h c Q h H T T = - (9)
with

( ) i i a a i g r a max k ; k Nu h h e = + ( 10
)
Heat flux through the left and right external vertical walls of the envelope is given respectively by:

( )

1 1 h h w w w T T Q k H e - = and
( ) ( ) N 1 2 c c w w w T T Q k H e + - = (11)
Similarly, the flux through a partition j is given by: ( ) ( )

j j 1 j c h p p p T T Q k H e + - = (12)
Finally, using the conservation of the flux in steady state, the total flux can be calculated from the temperature difference (T h -T c ) by:

( ) i h c tot N 1 p w w p g i 1 T T Q H e e 1 2 N k k h + = - = + + (13) i h T , i c
T and tot Q computation is made by solving iteratively equations [START_REF] Nishimura | Natural convection heat transfer in enclosures with multiple vertical partitions[END_REF] to [START_REF] Turkoglu | Natural convection heat transfer in enclosures with conducting multiple partitions and side walls[END_REF]. Thermal resistance of the part 2 is then given by:

( ) i N 1 p w w p g i 1 h c 2 tot e e 1 2 N k k h T T R Q H + = + + - = = (14)

Total resistance of the whole partitioned enclosure

The total resistance of the whole enclosure is determined with the thermal resistances R 1 , R 2 and R 3 in parallel. The total resistance for a unit of surface (r tot ) is therefore given by: ( )

1 2 3 tot w 1 2 1 3 2 3 R R R r H 2 e R R R R R R = + + + (15)

Two-dimensional numerical heat transfer model in partitioned enclosure

Actually, as soon as the convection takes place in alveolus, heat transfer is two-dimensional. The proposed numerical model consists in computing Navier-Stockes and energy equations. Using the Boussinesq approximation to account for density variation (in vertical direction), the steady state laminar non-dimensional continuity and momentum equations can be written as:

U W 0 X Z ∂ ∂ + = ∂ ∂ ( 16
)
1 2 2 2 2 2 1 2 2 2 2 2 U U P Pr U U U W X Z X Ra X Z W W P Pr W W U W X Z Z Ra X Z ∂ ∂ ∂ ∂ ∂ + = - + + ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + = - + + + θ ∂ ∂ ∂ ∂ ∂ (17)
( )

2 2 1 2 2 2 U W Ra Pr X Z X Z - ∂θ ∂θ ∂ θ ∂ θ + = + ∂ ∂ ∂ ∂ ( 18
)
where U and W are the dimensionless velocity components in the x and z direction, respectively. P and θ are the dimensionless pressure and temperature. Dimensionless parameters are defined as:

x X H = ; z Z H = ; 0 u U W = ; 0 w W W = with ( ) 0 0 h c W g T T H = β - 2 0 p P W = ρ ; 0 h c T T T T - θ = - with h c 0 T T T 2 + =
Pr is the Prandtl number defined as

a a Pr ν = α
Radiation between alveolus faces is calculated using the discrete ordinates method. Air in the alveolus is supposed semi-transparent non-diffusing medium with very small absorption coefficient.

Alveolus faces are assumed to be opaque and perfectly diffuse for the infrared radiation. The radiative flux leaving the surface of an alveolus face is given by: ( )

4 rad in S q = 1-q + T (19)
where the incident radiative heat fluxψ ψq in can be given by solving the radiative transfer equation (RTE) [START_REF] Colomer | Three-dimensional numerical simulation of convection and radiation in a differentially heated cavity using the discrete ordinates method[END_REF].

In solid parts the two-dimensional steady state heat transfer is governed by:

2 2 2 2 0 X Z ∂ θ ∂ θ + = ∂ ∂ ( 20
)
with the following boundary conditions:

( )

1 0, Z 2 θ = ; L 1 , Z H 2 θ = -; Z 0 0 Z = ∂θ = ∂ ; w 2e Z 1 H 0 Z = + ∂θ = ∂ (21)
A uniform, tight (401 x 603) (X,Z) grid -ensuring independence of the solution on the grid size -is used. The set of equations was computed using the computational fluid dynamics software Fluent [START_REF]Fluent, Fluent 6.1, User's guide[END_REF], based on the finite volumes method. In order to compare the results with those of the analytical model, the mean temperature is calculated at alveolus faces by: ( )

w w e 1 H m e H (X) X, Z dZ + θ = θ (22)
Average Nusselt number in an alveolus is calculated from:

w w e H cv e cv ref 1 q (z) dz H Nu q + = (23)
where ref q is the average conduction heat transfer in the case where convection does not develop. cv q (y) is the local parietal flux, calculated as:

( ) i i ' cv a T (z) T (z) q (z) k (z) x - = ∆ ( 24
)
where T i (z) is the temperature of the face i and T i' (z) is the temperature of the very close node in the viscous air layer.

Similarly, radiative Nusselt number is defined by:

w w e H rad e rad ref 1 q (z) dz H Nu q + = (25)
where rad q (y) is the local radiation flux at the alveolus face.

In order to validate the numerical code, one reproduces in Table 1 results of Wang et al. [START_REF] Wang | Étude numérique du couplage de la convection naturelle avec le rayonnement de surfaces en cavité carrée remplie d'air[END_REF] for an air-filled square cavity differentially heated in which natural convection is coupled with radiation.

Results of other authors [START_REF] Abdelbaki | Simulation numérique des transferts thermiques couplés à travers les parois alvéolaires des bâtiments[END_REF] on a same square cavity in which only convection is studied are also reported. The agreement is good: the maximal deviation with [START_REF] Wang | Étude numérique du couplage de la convection naturelle avec le rayonnement de surfaces en cavité carrée remplie d'air[END_REF] in Nusselt number is 1%.

Experimental study

A PVC partitioned enclosure has been built (Fig. 2a). Internal dimensions of the enclosure are

Validation of the analytical model

Although all the 7 configurations of the enclosure have been numerically and experimentally studied, only results of configurations G (N=3, T h -T c =23.2 K) and D (N=11, T h -T c =29.7 K) are presented in detail, since they are characteristic of two different flow regimes (Fig. 4) in air-filled cavities.

Numerical validation

Tables 3 and4 report convective and radiative analytical Nusselt numbers. In the enclosure G, Nu cv >>1 (Table 3) in all alveolus indicates that convection develops. Radiative flux represents 2/3 of the total heat transfer. In the enclosure D, convection is slight with cv Nu 1 ≅ (Table 4). Radiative flux is still higher than the convective one. These results show the predominance of radiative flux over convection in alveolus. This is due to the high emissivity of the material (PVC) used here and to the low local Rayleigh number in cavities (about 1.9 10 5 for case G and 2.3 10 3 for case D). Therefore one can not study this type of enclosure without taking into account the radiation [START_REF] Nicolas | Etude numérique des transferts de chaleur couplés dans une paroi alvéolaire[END_REF]. When convection decreases, the temperature stratification decreases on the faces of cavities and the difference between Nu cv given by the two models decreases and tends to zero when convection stops.

Analytical, numerical and experimental dimensionless mean temperatures at alveolus faces are plotted in Fig. 6. These graphs show a great agreement between analytical and numerical temperatures. Indeed, the maximal difference between analytical and numerical temperatures is about 0.1K for both enclosures G and D.

Except the values of Nu cv presenting a non negligible difference when convection develops, the numerical results match the analytical ones; this is a first validation of the analytical model.

Experimental validation

Fig. 6 compares also analytical and experimental dimensionless temperatures at alveolus faces.

These graphs show a good agreement between analytical and experimental results. The maximal difference between analytical and mean experimental temperatures is 0.6K for enclosure G and 0.9K for enclosure D. Although these maximal values are higher to the uncertainty of the thermocouples estimated to 0.3K, the results shown here can be considered as another validation of the analytical model. that the third dimension according y does not intervene in the heat transfer in the enclosure even with a small partitions number. This remark justifies limiting the numerical study to two-dimensional model.

Comments

The analysis of numerical and experimental results allows us to use with confidence the analytical model presented in 2.2, to describe heat transfer through the partitioned enclosure.

Indeed, this analytical model is a simplified expression of the two-dimensional modeling of the thermal transfer in enclosures. Radiation heat transfer calculation takes into account the geometry of the alveolus through the view factor F ij . The aspect ratio is sufficiently high to neglect the horizontal side effects. The two-dimensional aspect of the convection heat transfer is taken into account with the aspect ratio in Bejan relation [START_REF] Bejan | Convection Heat Transfer[END_REF]. The radiation has the effect to attenuate temperature stratification on faces and then justifies the validity of the use of Bejan relation. In the case of conduction heat transfer -preponderant when convection does not occur -only the presence of horizontal walls causes the two-dimensionality of heat transfer. The thickness of these walls is very low in relation to the height of the cavity; therefore, the heat flux through these horizontal walls is very small in comparison with the one through the partitioned enclosure.

The model is therefore validated, regardless of the number of partitions placed in the enclosure. A later validation will be performed with the determination of the number of partitions that maximizes the thermal resistance of the material.

Application: optimization of the thermal resistance

Maximal thermal resistance

Thermal resistance through alveolar material as described above varies with the number of partitions. Indeed, the increase of the number of partitions reduces radiative and convective flux in alveolus. On the other hand, this increase furthers conduction in the solid parts of the enclosure.

Therefore limiting heat transfer through the material implies to find the optimal partitions number leading to the maximal thermal resistance. This calculation is made theoretically using the analytical model of the partitioned enclosure that has been validated. Experimental determination is made by measuring the thermal resistance of the enclosure in configuration A to G.

Fig. 9 shows the variation of the thermal resistance versus the number of partitions. Bars of experimental resistance uncertainty are added in the graph. This curve shows a good agreement between theory and experiments. The existence of a maximum indicates that there is an optimal number of partitions that induces a maximal thermal resistance.

The maximal thermal resistance achieves a compromise between two phenomena: the conduction and the radiation heat transfer. Indeed, the more the number of partitions, the more conduction heat flux since the partitions are more conductive than air. Thus, the conduction thermal resistance decreases when N increases. On the other hand, the addition of a partition causes a reduction of the temperature difference in cavities and therefore the radiation heat transfer decreases. Then the radiation thermal resistance increases with N. The thermal resistance being the sum of two values varying in an opposite way with N, it exists an optimal number of partitions for which this resistance is maximal.

This number of partitions is equal to 45 in our case (PVC enclosure). The optimal number of partitions obtained with the analytical model is constant if the temperature difference (T h -T c ) varies over a constant mean temperature value; it increases slowly with the mean temperature. The optimal number varies with different parameters as shown in following.

Case of a construction material: influence of parameters

The validated analytical model has been then used to determine the thermal resistance of an enclosure made with a construction material. A parametric study for the partitioned enclosure has been made according to the thermal conductivity of the material of the enclosure (Fig. 10a), the emissivity of alveolus faces (Fig. 10b), the thickness of partitions (Fig. 10c) and the vertical wall thickness of the enclosure (Fig. 10d). The geometry of the enclosure is the same described in Fig. 1. The considered parameter only varies.

Thermal conductivity of the partitioned enclosure

Fig. 10a illustrates the influence of thermal conductivity of the material of the partitioned enclosure on the optimal number of partitions. The graphs outline two phenomena. First, the optimal number of partitions slightly decreases to a limiting value when the thermal conductivity increases. The limiting value is 40 partitions corresponding to an alveolus thickness of about 0.004 m for our case.

Material thermal conductivity has a small influence on the optimal number of partitions (from 40 to 43 in Fig. 10a), but of course a great one on the thermal resistance which increases (from 1.2 to 3.6 ) when thermal conductivity decreases. Secondly, when k is high, the variation of the thermal resistance presents a tray around the optimal value. Indeed, as soon as the convection stops in cavities, the reduction of radiation by addition of a partition has the same order of magnitude than the increase of conduction flux. Therefore, the total heat flux varies slowly. For instance, for the case of

- - = 1 1 k 2.
2 W.m .K , the optimal number could be reduced of half without the thermal resistance is modified.

Walls and partitions emissivity

In order to model a classic brick, material is considered to be terra cotta (k=1.04 W.m -1

.K -1

). Fig. 10b shows that the diminution of the alveolus faces emissivity improves the optimal thermal resistance of the enclosure. This result is predicable due to the fact that radiation represents the highest part of heat transfer in alveolus. It is worth noting that the optimal number of partitions decreases quasi linearly with the emissivity. The emissivity variation from 0.9 to 0.1 decreases the optimal number of partitions from 40 to 16 corresponding to an increase of the alveolus thickness from 0.004 m to 0.014 m. The aspect ratio A decreases from a value of 46 to 13.2. The diminution of the partitions number can constitute an economy on the material used for the manufacture.

Thickness of partitions

Considered material is still terra cotta. The partitions thickness has a strong influence on the thermal resistance of the enclosure when the partitions number is high (Fig. 10c). This influence is weak when the partitions number is low because the quantity of material is small in this case. The optimal resistance decreases when the partitions thickness increases. This result is predicable because by increasing the partitions thickness one increases the contribution of conduction in the solid parts into the total heat transfer through the whole enclosure. The use of thin partition improves the thermal resistance. For example, using partitions of 0.001 m of thickness instead of 0.003 m partitions increases the optimal resistance of about 20%.

Thickness of vertical external walls

Fig. 10d shows that the increase of the thickness of external vertical walls of the partitioned enclosure has a modest influence on the resistance. Indeed, the use of vertical walls of 0.020 m of thickness instead of 0.008 m reduces the maximal thermal resistance only about 4%. This result is important because it shows that one can improve the thermal inertia of the partitioned enclosure by increasing only the vertical wall thickness without degrading very significantly the thermal resistance.

Conclusion

  cold faces temperatures of the alveolus i, respectively.

εT

  are the emissivity coefficients of hot and cold faces of alveolus i, respectively; F hc is the view factor; i m is the mean temperature in the alveolus i.

  H=0.184 m, L=0.285 m and 0.275 m of depth. External walls have a thickness e w =0.008 m. 1 to 47 vertical partitions of thickness e p =0.003 m are placed between the vertical walls of the enclosure, creating therefore different configurations from 2 to 48 equal width alveolus. Thermal properties of the solid parts are the thermal conductivity of the PVC [2] ( faces emissivity ( =0.9). Physical air properties are β i = 3.49 10 -possible, thermocouples (type K) are placed on each partition at points A, B, C and D as indicated on Fig. 2b. Points of measures A to C permit to highlight the temperature stratification on Oz direction. Measures on points B and D are used to check the two-dimension of the problem. The thermocouples are standardized and the temperature uncertainty is estimated to 0.3K. Heat transfer through the partitioned enclosure is experienced in a calibrated hot-box (Fig. 3). Two highly insulated chambers are clamped tightly together to surround the test enclosure. The cold and hot chambers are metering and thermally controlled. Two fluxmeters (100 x 100 mm 2 ) are placed on the hot and cold sides of the enclosure. The fluxmeter precision given by the constructor is 5%. The climatic chamber is the one used for U-value measurements of a whole wall according to the standard ISO 8990 [19]. The apparatus has been modified for the characterization of a single brick [20]. Temperatures and parietal flux are regularly measured with a datalogger driven by a computer until the steady state is reached. The experimental thermal resistance of the partitioned enclosure is determined by: hot face temperature, exp c T is the average cold face temperature and q exp the average parietal flux measured on the hot face and on the cold face. Table 2 reports the problem parameters of the seven configurations of the partitioned enclosure experimentally studied, i.e. the partitions number N, the aspect ratio A of a cavity and the temperature difference (T h -T c ) of the whole enclosure. The temperature difference values are imposed by the experiments. These configurations have been chosen to characterize different modes of heat transfer, as can been shown on Fig. 4. Indeed, according to Yin et al. [21] and Batchelor [22], pure conduction occurs in the cavities of configurations A to D while established convection regime is expected for configurations F and G. Configuration E is in the transition zone.

  The analytical model allows to compare contributions of radiation and convection. If temperature difference in cavities is low, i.e. T h close to T c or N high (configurations A to D), convection does not occur. Only conduction heat transfer occurs in air cavities. Radiation and the conduction heat transfer vary simultaneously with the temperature difference in cavities and the ratio between them remains more or less constant. Moreover, if the temperature difference in cavities increases sufficiently -it occurs theoretically when (T h -T c ) reaches a limit value and/or when N is low (configurations F and G) -convection takes place in cavities. The convection heat flux increases more rapidly than the radiative one until the two modes balance. In Fig. 5 numerical isotherms contours are plotted for configurations G and D, respectively. In addition to the expected two-dimensional flow, these patterns also confirm the complete establishment of the convection in the enclosure G contrarily to the enclosure D, in accordance with the prediction of the diagram on Fig. 4. Tables 5 and 6 present radiative and convective numerical Nusselt numbers in the alveolus of the enclosures G and D, respectively. The numerical values of Nu rad are very closed to the analytical ones. The maximal difference is 1.3%. On the other hand, the numerical values of convection Nusselt numbers are lower than those obtained analytically. The maximal difference is 25% for enclosure G and 18% for enclosure D. This maximal difference takes place in the first and last cavities for the two cases. The analytical model overestimates the convection heat flux in cavities. It is due to the fact that Bejan correlation used here is established for cavities with isotherm faces.

Fig. 7

 7 Fig. 7 presents dimensionless temperatures at points A, B and C placed on a vertical line on partitions for the partitioned enclosure G and D. These graphs show measured temperature stratification on alveolus faces, confirming the two-dimension of the temperature field. The maximal temperature difference between points A and C is 2.9 K for enclosure G and 2.5 K for enclosure D. The mean value of this difference is 2 K for enclosure G and 1.5 K for enclosure D. These values confirm that convection is more established in enclosure G than in enclosure D. In order to determine the importance of the three-dimension of the problem, Fig. 8 presents dimensionless experimental temperatures at points B and D placed on a horizontal line on partitions for the partitioned enclosures G and D . These temperatures are practically the same. Indeed, the maximal temperature difference between points B and D is 0.3 K for enclosure G and 0.9 K for enclosure D. For enclosure D, only two values of the temperature difference between points B and D are superior to the thermocouple uncertainty. These values are erroneous. These results show

AFig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :Fig. 10 :Fig. 1 :Fig. 2 :Fig. 3 :Fig. 4 :Fig. 5 :Fig. 6 :Fig. 7 :Fig. 8 :Fig. 9 :Fig. 10 :

 234567891012345678910 Figures

Table 1 :

 1 Comparison of Nusselt numbers for a differentially heated square cavity (Ra=10

	Tables
							6 )
			De Vahl	Le Breton et	A. Abdelbaki and	H. Wang et	This study
			Davis [21]	al. [21]	Z. Zrikem [21]	al. [20]
	ε=0	Nu cv	9.270	8.794	9.446	8.852	8.947
	ε=0.2	Nu cv	N/A	N/A	N/A	2.337	2.338
		Nu rad	N/A	N/A	N/A	8.399	8.458
	ε=0.8	Nu cv	N/A	N/A	N/A	7.873	7.941
		Nu rad	N/A	N/A	N/A	11.208	11.229

Table 2 :

 2 Configurations of the partitioned enclosure -T c 30.2 29.4 23.2 29.7 21.4 28.8 23.2

	Name	A	B	C	D	E	F	G
	N	47	23	15	11	7	5	3
	A	61.3 20.4 12.3	8.8	5.6	4.1	2.7
	T

h

Table 3 :

 3 Analytical Nusselt numbers for the enclosure G (N=3,T h -T c =23.2 K)

	Alveolus	1	2	3	4
	Nu rad	11.5 10.9 10.3	9.7
	Nu cv	5.8	5.9	6.0	6.0

Table 4 :

 4 Analytical Nusselt number for the enclosure D (N=11, T h -T c =29.7 K)

	Alveolus	1	2	3	4	5	6	7	8	9	10	11	12
	Nu rad	4.2	4.1	4.0	3.9	3.8	3.7	3.7	3.6	3.5	3.4	3.3	3.2
	Nu cv	1.4	1.4	1.4	1.5	1.5	1.5	1.5	1.5	1.5	1.5	1.5	1.5

Table 5 :

 5 Numerical Nusselt numbers for the enclosure G (N=3, T h -T c =23.2 K)

	Alveolus	1	2	3	4
	Nu rad	11.4 10.8 10.2	9.7
	Nu cv	4.7	4.9	5.0	4.8

Table 6 :

 6 Numerical Nusselt numbers for the enclosure D (N=11, T h -T c =29.7 K)

	Alveolus	1	2	3	4	5	6	7	8	9	10	11	12
	Nu rad	4.2	4.1	4.0	3.9	3.9	3.7	3.7	3.6	3.5	3.4	3.3	3.2
	Nu cv	1.2	1.3	1.4	1.4	1.5	1.5	1.5	1.5	1.5	1.5	1.4	1.3


Björn Widarsson
email: bjorn.widarsson@mdh.se

Erik Dotzauer

Bayesian network based early-warning for leakage in recovery boilers

Keywords: Early-warning, Bayesian networks, Recovery boiler

Early-warning for leakage in a recovery boiler can help the process operator to detect faults and take action when a dangerous situation is developing. By analysing the massbalances on both the steam and combustion side of the boiler in a Bayesian network, the probability of leakage can be determined and used as an early-warning. The method is tested with real plant data combined with leakage simulations. The results show that it is possible to detect considerably smaller leakages using this method than using the type of simple steam-side mass-balance method that is in use today. Bayesian network is an efficient tool to combine information from measurement-signals and calculations giving an early-warning system that is robust to signal faults.

Introduction

Operating and diagnosing complex industrial processes are usually difficult tasks.

Faults that cause only minor disturbances in the production can be frequent. It is mainly the operator who detects, isolate and take action when a disturbance or fault appears. It can be very hard to distinguish between normal process operation with its minor disturbances and a developing fault.

Early-warning systems can help the process operator to detect faults at an early stage and thus prevent further fault development. Early warning systems have been developed for safety critical events, e.g. fire on ships [START_REF] Lee | Early Warning of Ship Fires Using Bayesian Probability Estimation Model[END_REF], disconnections in electrical systems [START_REF] Tsai | Development of Islanding Early Warning Mechanism for Power Systems[END_REF] and landslides [START_REF] Zan | Landslides Early Warning Monitoring System[END_REF]. Common for early-warning systems are that they are designed to make a person aware of that a critical situation is or might be developing.

Industrial processes normally have alarms on signal levels to pay attention to large deviations from a normal value. The alarm trigging levels are set to protect the process equipment, maintain product quality, etc. A diagnostic system can give the process operator more information than just an alarm that a signal is deviating. Fault diagnosis is commonly divided in the two steps; detection and isolation [START_REF] Venkatasubramanian | A review of process fault detection and diagnosis Part I: Quantitative model-based methods[END_REF]. Fault detection is to determine that there is a fault and isolation is to determine where the fault is located.

An early-warning system is a diagnosis system aimed to help the process operator to detect and isolate faults as early as possible.

Leakage detection in recovery boilers is important to avoid severe damages on equipment. The walls of the furnace are containing evaporating water with high pressure. Fireside corrosion and thermal stress can cause leakages, implying that water or steam comes in contact with the smelt. Water in the smelt can cause an explosion with total destruction of the boiler as result. There are a number of commercial systems for detection of leakage flows in both conventional boilers and recovery boilers [START_REF] Buckner | Early Steam Leak Detection Can Help Limit Boiler Damage, Avert Explosion[END_REF] [START_REF] Studdard | Early boiler leak detection cuts costs[END_REF].

These systems are detecting the leakage flow or cracking with acoustic sensors.

In this work, a new method to detect boiler leakage flows by a mass balance in a Bayesian network is evaluated by study of a typical recovery boiler. Data from real process operation are combined with a process model to simulate leakage flows.

Bayesian networks [START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF] [START_REF] Pearl | Probabilistic reasoning in intelligent systems[END_REF] have been used for diagnostics in many different areas e.g. medicine, electronics and mechanics [START_REF] Przytula | Construction of Bayesian networks for diagnostics[END_REF]. An advantage over other diagnostic methods is the handling of uncertainty. An alternative would be to use e.g. fuzzy logic [START_REF] Venkatasubramanian | A review of process fault detection and diagnosis Part II: Qualitative models and search strategies[END_REF]. This methodology will be investigated in future research on the current application.

Chemical recovery

Black-liquor is a residue from the cooking-plant in sulphate pulp mills. It consists of both organic and inorganic materials and has a dry solids content of 14-18 %. The organic material contains energy and the inorganic material gives valuable process chemicals. Before it can be combusted, the dry solid content has to be increased in an evaporation plant. Content of dry solids is 65-80 % after the evaporation plant. 4

Black-liquor combustion

The recovery boiler burns the organic material in the concentrated black-liquor and produces high pressure steam [START_REF] Gullichsen | Papermaking Science and Technology -Chemical Pulping[END_REF]. It also recycles and regenerates chemicals through reduction, mainly from potassium-sulphate to potassium-sulphide. The recovery boiler is an important component in the process because it both regenerates the cooking chemicals and produces steam to processes. Concentrated black-liquor is fed in a furnace where volatile organic material is combusted. The remaining material falls down to the bottom and forms a smelt.

Reducing conditions in the lower part of the boiler reduces the sulphates to sulphide with some of the energy generated in the combustion. The smelt is also containing a number of other components, mainly sodium carbonate, which are not active in the process. Smelt spouts bleeds the smelt out of the furnace to be solved in weak-liquor.

Steam generation

Water is fed from a feed-water deaerator thorough the economiser. The economiser is located as the last component in the flue gas path and utilises the reminding heat in the flue gas to preheat the feed water before it enters the steam drum. Furnace walls and a convection section in the flue gas path are cooled by evaporating water. From the steam drum, the furnace walls and the convection section are connected in a natural circulation system, see Figure 1.

Steam is separated in the upper part of the steam drum and lead to the super-heaters.

Super-heaters are located first in the flue gas path, close to the furnace. To protect the super-heater closest to the furnace from radiation, water or steam cooled screen-tubes can be located in front of the super-heater. Steam cooling to control the outgoing steam temperature of each super-heater section is done with injection of water or by a surfacecooler. Injection water can be feed-water or, to prevent salts and other contaminations in the steam, water produced by condensing steam from the drum thorough cooling with feed-water [START_REF] Gullichsen | Papermaking Science and Technology -Chemical Pulping[END_REF].

Leakage

Corrosion, erosion and thermal stress can cause holes and cracks resulting in a flow of water or steam into the combustion side of the boiler. The leak can be located in either the water or steam tubes, but a leak in the furnace wall-tube is the most severe due to the risk of water to come in contact with the smelt. A leak in the economiser results mainly in a water flow, with some flashing steam, while a leak in the super-heaters results in a pure steam-flow.

A smaller leak from a corrosion hole can cause erosion on tubes close to the leaking tube. The erosion can then lead to a larger tube rupture with an extensive leakage flow as result. The magnitude of a leakage can be varying. Leakages have in other models been simulated with flows from 0.2 kg/s [START_REF] Sun | Efficient model-based leak detection in boiler steam-water systems[END_REF] [13].

Early-warning

Leakage detection by analysing the mass balance on the steam side of the boiler is only possible when a large leakage flow is present. This due to the relatively low measurement precision on feed-water and steam-flow in comparison to a leakage flow.

A leakage can also be detected on the combustion side by indirect calculations of mass flows, but the precision is lower with this method compared to the method exploiting the mass balance on the steam side. By combining the two balances, indications on a leakage can be considered from both the steam side and the combustion side. This is exploited here to generate an early-warning to the process operator.

Bayesian networks

In this paper, we use a Bayesian network for diagnostics. A Bayesian network is defined from [START_REF] Jensen | Bayesian Networks and Decision Graphs[END_REF];

• a set of nodes representing variables,

• directed arcs between the nodes, and

• a probability table for each node.

If there is an arc from node A to node B, B is a child of A and A is a parent to B. The nodes are originally discrete, but handling of continuous variables is possible if they are not parents to discrete nodes. Discrete nodes have two or more states. A conditional probability table (CPT) describes the probability for the states given all possible combinations of the parent states. This means that a combination of many parent nodes with a lot of states gives a large CPT. An example of a CPT is shown in Figure 2, which is describing the conditional probabilities for node B depending on the states of node A and C. Continuous nodes have instead of discrete states a Gaussian distribution with a mean and standard deviation. Analogous to the CPT for a discrete node, the mean and standard deviation of the continuous node are computed as functions of its parents.

The nodes in a Bayesian network for diagnostics can be classified as fault nodes, intermediate nodes and symptom nodes, see the illustrating example in Figure 3. Input data are inserted as evidences of a node state in the symptom nodes. The network is then propagated, which means that the probabilities are computed for all nodes. After the propagation, the probabilities of the fault nodes form the output.

Calculation of mass flows

The mass-balance on the combustion side is defined as (

) 2
The feed-water flow (

) 9
The specific heat capacity is determined from flue gas temperature and composition.

Signal diagnostics

Uncertainty in process signals and calculated values are represented in the Bayesian network as a signal or calculation node. The node is connected as a child to the node representing the real value together with a node representing the signal or calculation status [START_REF] Weidl | Root cause analysis for pulp and paper applications[END_REF], see Figure 4.

The signal status node has the states "normal" and "faulty". If the state of the signal status node is "normal", the process signal node is Gaussian distributed with a mean of the real value and a standard deviation representing the measurement precision. When the signal status is "faulty", the standard deviation in the signal node is changed to a lager value, covering the entire signal range. For the current process signals it was found suitable to increase the standard deviation with a factor of ten.

Mass balance in Bayesian network

The mass balances on both the combustion side (1) and the steam side (2) of the boiler are defined in a Bayesian network. All mass flows are represented with continuous nodes and forms a structure of intermediate nodes.

A discrete node with the states "none" and "leakage" connects the two mass balances, see Figure 5. The upper part is representing the steam-side balance and the lower part the combustion side balance. The soot-blow flow is also connecting the two mass balances. When the state of leakage is "none", the node has no impact on flue-gas and high-pressure steam. When the state is "leakage", a flow is added to the mean value of the flue-gas and subtracted from the mean value of high-pressure steam.

In particular (see Figure 5), equation ( 1) is modelled in the node calculated values are interpreted as symptoms. Leakage and status nodes are defined as fault nodes. If the signal and calculated values are added into the network, and the network is propagated, the probabilities of leakage and signal status are obtained. The leakage probability is then used to assess if there is a possible leakage or not.

If available, other observations such as acoustic, blow-down analysing or image leak detection can be inserted as nodes connected to the leakage node. The studied recovery boiler has none of these leak detection systems installed.

Results

The Bayesian network defined in Section 3.4 is applied to process history data from a recovery boiler at Vallviks Bruk, Sweden. The data is from operation at varying load conditions during a time period of three months.

Three different cases have been evaluated;

• normal operation,

• test with simulated leakage, and

• test with fault in signals.

Since data from operation with a leaking boiler is not available, a simulation model has been developed. The software Prosim [START_REF] Endat | [END_REF] is used to build a model with possibility to simulate leakage in four components; furnace wall, super-heater, convection section and economiser. Prosim is a thermal power plant simulator. The software includes general modules for conventional thermal power plants and specific modules for recovery burning.

Normal operation

A test with data from normal process operation is verifying the indirect calculation of liquor flow, etc. The results are shown in Figure 6. Case number 1-29 are from normal load operation and case number 30-34 are from low load operation.

The leakage probability is below 7 % during normal load conditions. In the low load cases, the probability is between 13 and 18 %. This indicates that the indirect fuel mass flow calculation is not enough accurate under low load conditions.

Simulated leakage

Real plant operation data are used as input to the model and a leakage is simulated.

Output from the model is used to determine the probabilities of lakage for a leakage flow of 0.25, 0.5 and 1 kg/s, respectively. Figure 7 depicts the results.

A calculation without leakage gives a leakage probability of about 5 %. Simulation with 0.25 kg/s leakage results in a significant leakage probability for all components (above 20 %). The probabilities increase when the leakage flow increases to 0.5 kg/s. If a leakage of 1 kg/s is simulated in the economiser or convection section, the leakage probabilities are continuing to increase. A simulated leakage of 1 kg/s in the furnace or super-heater is resulting in a decreasing leakage probability. The reason for this may be that the network finds it more probable that a signal is faulty when the leakage is in the furnace or super-heater. The performance of the system is shown in Table 1, where the average leakage probability of normal operation is compared to cases with 0.25 kg/s leakage flow and cases with signal faults of ±2 %. In the table, we see that the leakage probability is below 15 % during normal operation for both normal and low load. A leakage of 0.25 kg/s results in a leakage probability of about 25-30 %. A signal fault may cause a high leakage probability when the current sensor is measuring feed-water flow or blackliquor dry substance.

Fault in signals

Conclusions

Early-warning for fault detection in a recovery boiler was considered. A methodology based on mass balances in a Bayesian network was developed. The early-warning system was tested on real plant data from normal operation combined with simulations to derive leakage cases.

The results show that normal load operation gives a leakage probability below 10 %. A simulated leakage of 0.25 kg/s gives a leakage probability of about 25 %. This means, during normal operation, that an early-warning can be generated when the leakage probability exceeds 10 % without generating false alarms. A fault in measurement signals can cause a leakage probability of more than 10 %, but only when the relative fault in the signal is large. The detection level and reliability are depending on the measurements on feed-water flow and black-liquor dry substance.

A limit for early-warning when the leakage probability exceeds 20 % can help the operator to be aware of that something is wrong, possibly a leak in the boiler. Further investigations can then be done to assess if it is a leakage or not. The ability to combine the mass-flow balance with other leakage detection methods in the Bayesian network is a possible way to further improve the warning sensitivity and reliability.

The method with mass balances formulated in a Bayesian network is an efficient tool to analyse the balances and warn the process operator if the leakage probability becomes too high. Integrating signal diagnosis in the network structure makes the system robust to minor signal faults. Early-warning for a recovery boiler is one application in which the method is applicable, but it can also be generalized to other applications, e.g. screen clogging, hang-ups or sintering.

  the sum of fuel flow ( fuel m ), air flow ( air m ) and soot-blow steam-flow ( blow soot m - ) equals to the sum of flue-gas flow ( gas flue m -) and smelt flow ( smelt m). On the steam-side, the mass balance is

  flows in and out from the boiler on the steam side, while only fuel is measured on the combustion side. This means that gas flue m -, air m and smelt m has to be computed by an indirect method. By a combustion and smelt reduction calculation on the combustion side, it is possible to determine flue-gas (g), combustion-air (a) and smelt quantity (s) per mass unit of fuel. The calculation is based on the ultimate analysis of the fuel, reduction rate (R), fraction of sulphur-oxide, [SO 2 ] flue-gas , in flue-gas and fraction of oxygen, [O 2 ] flue-gas , in flue-gas [11]. Solid fuel fired boilers utilise an indirect combustion calculation to determine flow of fuel, flue-gas, air and slag. This method is also applicable on boilers where the fuel flow is measured, to get a calculated fuel mass-flow. With a measured thermal power (P th ), a known lower heating value (LHV) of the fuel and the boiler efficiency ( boiler ), the fuel mass-flow is calculated as is calculated by determining the relative losses in the flue gas (l flue-gas ), unburnt in flue gas (l CO ) and heat in the slag (l slag ). Heat losses to the surroundings (P radiation ) are constant heat flows determined by the size and type of boiler. Nominal boiler power (P nominal ) is used to get the relative loss. The boiler efficiency is then given as similar to conventional boilers in the sense that the smelt is equivalent to the slag. Slag losses are often negligible in conventional boilers since the ash content in the fuel is low. The flow of smelt in a recovery boiler is considerably higher and it leaves the boiler at a high temperature. Instead of an effective heating value, a net heat value (NHV) is used for recovery boilers in (3). The net heat value is the lower heating value for the black-liquor corrected with the energy to the reduction process [11], for reduction (Dh R ) is 13.1 MJ/kg Na2S and C S is the sulphur content. The efficiency of reduction (R) is expressed as [11], of combustion-air ( air m ) are then calculated using the fuel mass-flow and air quantity from the combustion calculation, way, the mass-flow of smelt ( smelt m ) is computed from smelt quantity and fuel mass-flow, application, the combustion-air ratio (a) is about 4.1 kg air per kg fuel and the smelt ratio (s) is about 0.26 kg smelt per kg fuel. The flue-gas flow is calculated using an energy balance over the economisers. Heat transferred (P economiser ) is determined on the feed-water side. Specific heat capacity (c p flue-gas ) and temperature difference of the flue gas (t flue-gas in -t flue-gas out ) gives the flue-gas mass

  3) -(9) are parts of the nodes representing calculated values. For all nodes representing signals and calculated values, the status handling described in Section 3.3 is applied.

A

  bias is added to the process signals one at the time to simulate a fault in the signal, see the results in Figure 8. The feed-water flow and black-liquor dry substance is the most important signals. A fault in one of these signals can make the early-warning system indicating leakage probabilities over 15 % and thereby generating false warnings. Faults in the other signals result in a leakage probability below 10 %.

Figure 1 .

 1 Figure 1. Recovery boiler.

Figure 2 .

 2 Figure 2. Bayesian network and the conditional probability table for node B.

Figure 3 .

 3 Figure 3. Principal construction of Bayesian network for diagnostics.

Figure 4 .

 4 Figure 4. Real value node and signal status node connected to process signal node.

Figure 5 .

 5 Figure 5. Bayesian network for diagnosis of a recovery boiler.

Figure 6 .

 6 Figure 6. Test with data from real cases.

Figure 7 .

 7 Figure 7. Test with simulated leakage flows in four positions.

Figure 8 .

 8 Figure 8. Test with simulated fault in signals.

Leakage


Anastasios I Stamou
email: stamou@central.ntua.gr

Ioannis Katsiris

Alois Schaelin

EVALUATION OF THERMAL COMFORT IN GALATSI ARENA OF THE OLYMPICS "ATHENS 2004" USING A CFD MODEL

Keywords: Computational Fluid Dynamics (CFD), indoor stadiums, mathematical models, thermal comfort

A Computational Fluid Dynamics (CFD) model was used to evaluate the thermal comfort conditions in the indoor stadium of the Galatsi Arena, which hosted the sports of rhythmic gymnastics and table tennis during the Olympic Games "Athens 2004". The CFD code CFX was applied to calculate the 3-D airflow and temperature fields in the Arena for various values of temperatures (T in ) of conditioned inlet air. Calculated mean velocities and temperatures were used to determine the thermal comfort indices PMV (Predicted Mean Vote) and PPD (Predicted Percentage of Dissatisfied) and to evaluate the thermal conditions in the various regions of the Arena. Calculated PMV and PPD

values showed that thermal conditions in the Galatsi Arena were very satisfactory for T in =16 o C; only a small percentage (less than 7 %) of the spectators was expected to be slightly uncomfortable.

Introduction

In the last two decades Computational Fluid Dynamics (CFD) models have been increasingly used for the calculation of airflow velocities and temperatures in indoor environments for the evaluation of comfort conditions, smoke conditions and air quality. A significant number or scientific papers exists dealing with the application of CFD models in various indoor environments, such as apartments [START_REF] Bojic | Locating air-conditioners and furniture inside residential flats to obtain good thermal comfort[END_REF], offices [START_REF] Hien | Effects of double glazed facade on energy consumption, thermal comfort and condensation for a typical office building in Singapore[END_REF][START_REF] Stamou | Verification and application of a CFD model for the evaluation of thermal comfort in office spaces[END_REF], museums [START_REF] Papakonstantinou | Computational analysis of thermal comfort: the case of the archaeological museum of Athens[END_REF], lecture theatres [START_REF] Cheong | Thermal comfort study of an air-conditioned lecture theatre in the tropics[END_REF], classrooms [START_REF] Karimipanah | Investigation of air quality, comfort parameters and effectiveness for two floor-level air supply systems in classrooms[END_REF], clinics [START_REF] Helmis | Indoor air quality in a dentistry clinic[END_REF], industrial premises [START_REF] Rohdin | Numerical predictions of indoor climate in large industrial premises, A comparison between different k-models supported by field measurements[END_REF] and car parks [START_REF] Chow | Numerical studies of airflows induced by mechanical ventilation and air-conditioning (MVAC) systems[END_REF]. Moreover, CFD models have been combined with other tools or methods for the design of indoor spaces. For example; Desta et al [START_REF] Desta | Combining CFD and data-based mechanistic (DBM) modeling approaches[END_REF] combined a CFD model with a low order Data Based Mechanistic (DBM) model to predict air temperatures in a full scale, adiabatic walled ventilated test chamber. The CFD model supplied data and the DBM model performed model order reduction; in this way the two modeling approaches supplemented each other symbiotically to eliminate their disadvantages.

Zhai and Chen [START_REF] Zhai | Sensitivity analysis and application guides for integrated building energy and CFD simulation[END_REF] analyzed the potential building and environmental characteristics (such as environmental conditions, HVAC systems, building occupying and operating conditions, envelope properties and building sizes) that may influence the necessity and effectiveness of applying Energy Simulation-CFD coupling simulation. They conducted a sensitivity analysis of the coupling simulation in a representative office building in Boston and provided general suggestions on appropriate development and usage of the coupling simulation. Kim et al [START_REF] Kim | Two-step optimal design method using genetic algorithms and CFD-coupled simulation for indoor thermal environments[END_REF] presented a two-step optimal design investigation method using Genetic Algorithms for CFD indoor thermal environments and applied the optimal design in an office space.

In spite of the extensive application of CFD models in various indoor environments, there are only a few studies in indoor stadia. In 2001 a CFD model was applied for the Yoyogi National Stadium in Tokyo [START_REF]FLUENT News[END_REF]. The Stadium was constructed for the Olympics of 1964. Since then, it has continued to be used for the performance of swimming competitions or concerts; the CFD model was used to evaluate indoor conditions for both modes of stadium operation. Very recently, a CFD study was performed in the Olympic ice hockey stadium in Turin to evaluate fire hazards; smoke conditions were reviewed in the stadium for various scenarios and uses [START_REF] Chown | The Olympic ice hockey stadium[END_REF].

The present work demonstrates the application of a CFD model to calculate airflow velocities and temperatures in the Galatsi Arena, shown in Figure 1; the Arena hosted the sports of rhythmic gymnastics and table tennis during the Olympics of 2004.

The calculated velocities and temperatures were used to determine the main thermal comfort indices Predicted Mean Vote (PMV) and Predicted Percentage of Dissatisfied (PPD) [START_REF]Moderate Thermal Environments -Determination of the PMV and PPD Indices and Specification of the Conditions for Thermal Comfort[END_REF]; these were subsequently used for the evaluation of the thermal conditions in the Arena.

The mathematical model 2.1 The CFD code

There are three main types of CFD methods: (1) Direct Numerical Simulation (DNS),

(2) Large Eddy Simulation (LES) and (3) Reynolds Averaged Navier-Stokes (RANS).

A short presentation of these methods can be found in Stamou and Katsiris [START_REF] Stamou | Verification and application of a CFD model for the evaluation of thermal comfort in office spaces[END_REF]. These methods are employed in various efficient computer codes, which are most frequently used for indoor CFD calculations, such as PHOENICS [START_REF]PHOENICS 3.2 (software package)[END_REF], FLUENT [START_REF]Fluent user' s guide, Version 5[END_REF] and CFX [START_REF]CFX-5[END_REF]. In the present work, the latest version of the computer code CFX [START_REF]CFX-5[END_REF] was used.

The CFX code calculates the 3-D flow field and heat transfer using the continuity, momentum and energy equations, which are written as follows:

j j U = 0 x ρ ∂ ∂ (1)
j i i i i j i j j P U U U U + = - + ( RTS) + g t x x x x ρ ρ µ ρ ∂ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ∂ (2) j j j j T T T U + = ( TTS) t x x x ρ µ ∂ ∂ ∂ ∂ + ∂ ∂ ∂ ∂ ( 3
)
where t is the time, x i is the Cartesian coordinate in the i-direction, U i is the flow velocity in the i-direction, is the density of water, T is the temperature, P is the pressure, is the molecular viscosity of the water and g i is the acceleration of gravity.

For the calculation of the Reynolds (turbulent) stresses (RTS), the assumption of the isotropic turbulence is applied combined with the Bussinesq approximation, i.e.

j i ij t j i 2 U U RTS = ( + ) - k 3 x x ρ µ δ ∂ ∂ ∂ ∂ ( 4
)
where t is the eddy viscosity, ij is the Kronecker delta ( ij =1 for i=j and ij =0 for i≠ j )

and k is the average turbulent kinetic energy per unit mass, given by

2 2 2 3 1 2 1 k = ( + + ) u u u 2 (5)
Similarly, the Turbulent Thermal Stresses (TTS) are calculated by the following equation

t j TTS x µ σ Τ ∂ Τ = ∂ ( 6
)
where T is the turbulent Schmidt number for T.

The turbulence model

In the present work, turbulence is modeled with the Shear Stress Transport (SST) kbased model [START_REF] Menter | Two-equation eddy viscosity turbulence models for engineering applications[END_REF]. This model combines the standard k-model [START_REF] Rodi | Turbulence Models and Their Application in Hydraulics -A State of the Art Review[END_REF] and the k-model [START_REF] Wilcox | Reassessment of the scale determining equation for advanced turbulence models[END_REF]; virtually, it is a transformation of the k-to a k-formulation and a subsequent addition of certain equations.

The standard k-turbulence model [START_REF] Rodi | Turbulence Models and Their Application in Hydraulics -A State of the Art Review[END_REF] relates the eddy viscosity to k and the rate of its dissipation ( ) via the following equation.

2 t k = cµ ρ µ ε (7)
where c is an empirical constant and is given by the following equation

i i j j u u = x x µ ε ρ ∂ ∂ ∂ ∂ (8)
The distributions of k and are calculated from the following semi-empirical modeled transport equations

j t j j k j k k k U + = ( ) + G - t x x x ρ µ ρ ρ ε σ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (9)
1 2

2 j t j j j U + = ( ) + c G -c t k k x x x ε ρ ε µ ρ ε ε ε ε ρ σ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ( 10
)
where G is the production term of k by the mean velocity gradients given as follows

j i i t j i j U U U G = ( + ) x x x µ ∂ ∂ ∂ ∂ ∂ ∂ (11)
The standard values of the constants c =0.09, c 1 =1.44, c 2 =1.92 and the turbulent Schmidt numbers for k and , k =1.0 and =1.3, respectively, are used in the present computations.

The standard k-model [START_REF] Wilcox | Reassessment of the scale determining equation for advanced turbulence models[END_REF] relates the eddy viscosity to k and to the turbulence frequency ( ).

t k = µ ω (12)
The k and are calculated from the following model transport equations

j ' t j j k j k k k U + = ( ) + G - k t x x x ρ µ ρ ρβ ω σ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ( 13
) 2 j t j j j U + = ( ) + G - t k x x x ε ρ ω µ ρ ω ω ω α ρβω σ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ (14)
The standard values of the constants are '

Discretization and solution of the equations

The finite control-volume method is implemented for the spatial discretisation of the domain. The continuity, momentum and energy equations are integrated over each control volume, such that the relevant quantity (mass, momentum, energy, k, and ) is conserved, in a discrete sense, for each control volume. For the continuity equation (pressure-velocity coupling) a second order central difference approximation is used, modified by a fourth order derivative in pressure, which redistributes the influence of pressure. The second order upwind Euler scheme approximates the transient term.

The code employs an automatic, unstructured hybrid element mesh generator with an adaptive mesh refinement algorithm, which permits a very accurate representation of the boundaries. The advantage of using unstructured mesh is the minimization of numerical errors and the consistency of the solution throughout the domain. For the solution of the equations, a scalable and fully implicit coupled solver is used.

The Galatsi Arena

The Galatsi Arena is located in the district of Galatsi, about 8 km north of the Athens city centre and has a capacity of 6,000 spectators. The Arena, which measures 122.0 m long, 100.0 m wide and a maximum roof height of 30.0 m, features a very interesting design. Its roof takes the dynamic shape of the wing of an airplane (see Figure 1), while the large glass surfaces make sure that its interior is filled with sunlight all day. The Arena was designed by A. Tobazis and built by EMPEDOS SA at a cost of 43 million, while construction works commenced in November 2001. Heat fluxes were modeled to represent the amount of heat generated by the occupants, e.g. spectators and athletes, and the furniture in Galatsi Arena. Each simulated occupant was assumed to have a surface area of 1.6 m 2 and was modeled to generate convective heat equal to 70.0 W. This value corresponds to 43.7 W/m 2 [START_REF] Ashrae | ASHRAE Standard 55-1992R, Thermal environmental conditions for human occupancy[END_REF]. At the walls, floor and ceiling of the Galatsi Arena the "radiation" condition has been applied, which does not require the specification of temperatures at these boundaries.

In Figure 4 a view of the numerical grid is shown; due to symmetry half of the Arena was modeled. The grid consisted of approximately 4,000,000 unstructured elements (tetrahedrals, prisms and pyramids) with grid refinement on the surface of the heat sources and in the inlet and outlet regions. A sensitivity analysis was conducted to ensure that the resolution of the mesh was not influencing the results. The construction of the 3-D geometry and the relevant grid was a difficult and time-consuming task; it was performed with the CAD software Rhinoceros [START_REF] Rhinoceros | [END_REF] using as inputs the available 2-D architectural and mechanical engineering CAD drawings of the Arena.

To reduce the total computation time the initial air temperatures were set equal to the values determined by a simple heat balance assuming that the Arena is a completely mixed reactor; these values were calculated equal to 23.52 C and 25.08 C for T in =14 C and 16 C, respectively.

Results and discussion

Air velocities and temperatures

In Figure 5 this is due mainly to the interaction of the incoming airflow from the 256 diffusers with the buoyant plumes originating from the 6,000 occupants, who constitute a significant heat source in the Arena. Figure 5 depicts that the heavy airflow from the diffusers in the roof of the Arena is directed downwards, where it joins the airflows from the diffusers above the entrances; the resulting flow is subsequently "lifted"; then, is descends to lower levels and finally it is directed to the outlet openings. The air velocities decrease as the flow progresses downwards. Practically, the above-mentioned flow forms a clockwise recirculation region, which occupies a large part of the Arena; in the area of athletes almost completely mixed conditions are created. Two recirculation regions are also formed in the areas of corridors from the airflows of the diffusers in the roof of the Arena, which are directed towards the inner walls.

Figures 6, 7, 8 and 9 depict that the temperature field depends on the flow field and vice-versa. In spite of the approximately complete mixed regime in a large part of the Arena, the distributions of temperature in Figures 7 and8 show a general tendency to stratify; this tendency is interrupted locally by the airflows originating from the diffusers (mainly from these above the entrances), which create "islands" of regions with low air temperatures.

In the areas of the spectators air velocities range from 0.08 m/sec up to 0.64 m/sec for T in =14 C and from 0.08 m/sec up to 0.48 m/sec for T in =16 C; in some regions air velocities are much higher than the recommended limit of less than 0.25 m/s by the ISO Standard 7730 [START_REF]Moderate thermal environments-determination of the PMV and PPD indices and specifications of the conditions for thermal comfort[END_REF]. The ranges of air temperatures are calculated equal to 22-26 C and 24-27 C for T in =14 C and T in =16 C, respectively; low values of air temperatures are noticed in the eastern sector of the Arena.

In the area of the athletes air velocities range from 0.10 m/sec to 0.80 m/sec. The regions with high air drafts were identified and indicated to the organizers of the Games, because these may affect the motion of (a) the ball of the table tennis or (b) the ribbon of the rhythm gymnastics. The air temperatures are practically uniformly distributed in the area of the athletes due to the relatively high level of mixing; they range from 22 C to 24 C and from 24 C to 26 C for T in =14 C and T in =16 C, respectively.

Thermal comfort

ISO 7730 [START_REF]Moderate thermal environments-determination of the PMV and PPD indices and specifications of the conditions for thermal comfort[END_REF] defines thermal comfort as "the condition of mind that expresses satisfaction with the thermal environment". The reference to "mind" emphasises that comfort is a psychological phenomenon; therefore, it is "measured" using subjective methods. Over many years empirical research has related environmental conditions to physiological and subjective responses of subjects. Nowadays, the thermal comfort indices Predicted Mean Vote (PMV) and Predicted Percentage of Dissatisfied (PPD) of Fanger [START_REF] Fanger | Thermal Comfort, 1st edn[END_REF] are accepted as ISO 7730 [START_REF]Moderate thermal environments-determination of the PMV and PPD indices and specifications of the conditions for thermal comfort[END_REF]; these are calculated by two empirical equations. The PMV equation uses a steady-state heat balance for the human body and postulates a link between the deviation from the minimum load on heat balance effector's mechanisms, e.g. sweating, vaso-constriction, vaso-dilation, and thermal comfort vote; the greater the load, the more the comfort vote deviates from zero. The basis of this equation was obtained from experiments, in which the thermal sensation vote indicated the personally experienced deviation to the heat balance (-3 [cold] to +3

[hot]; seven point scale, 0 = neutral (optimum)). The PPD equation indicates the variance in the thermal sensation of the group of persons exposed to the same conditions. Dissatisfaction with the thermal environment, discomfort, was defined for those who voted cool (-2), cold (-3), warm (+2) or hot (+3). Under optimal thermal conditions (PMV = 0) a minimum of 5% dissatisfied is found, assuming identical activity levels, clothing and environmental conditions.

The PMV and PPD are calculated from knowledge of the so-called six basic variables: (1) activity, (2) clothing and (3) four environmental variables: air temperature, air velocity, mean radiant temperature and air humidity. Generally, the value of the PMV is very sensitive to the air temperature; the effect of the metabolism on PMV is also pronounced. In the present work the values of the variables for the activity (metabolic rate) and the clothing (ensemble insulation) were determined using ASHRAE Fundamentals [START_REF]ASHRAE Handbook Fundamentals[END_REF]; more specifically, the total thermal resistance of clothing was taken equal to 0. In the present work the factor of humidity was not considered, because people are quite insensitive to humidity levels over a wide humidity range -at the temperatures, which are normally found in air-conditioned spaces. Furthermore, the effect on comfort of a shift in humidity may be compensated for by a small adjustment of air temperature;

for example, higher humidity makes a person feel warmer, thus a slight lowering of temperature will compensate for the comfort effect of this higher humidity.

Conclusions

CFD calculations in the Galatsi Arena demonstrated that that the variable air volume conditioning system serving the Galatsi Arena is capable to ensure satisfactory thermal conditions for the spectators. Calculated PMV and PPD values showed that thermal conditions in the Arena are very satisfactory for T in =16 o C; only a small percentage of less than 7% of the spectators is expected to be slightly uncomfortable.

The temperature of T in =16 o C was proposed for application during the Olympic Games "Athens 2004". Moreover, the regions with high air drafts in the areas of athletes were identified and indicated to the organizers of the Games.

  = 0.09, =5/9, =3/40 and the turbulent Schmidt numbers for k and , k =2 and =2. The model is numerically stable, especially the low-Reynolds number version, as it tends to produce converged solutions more rapidly than the k-model. Furthermore, the Low-Reynolds (LR) number version is more efficient than the LR k-model in that it does not require calculations of wall distances, additional source term and/or damping functions based on the friction velocity. The main weakness of the k-model is its strong sensitivity to free-stream conditions; Gebremedhin and Wu [22] applied the standard k-model in an indoor environment together with four other turbulence models: namely (1) the standard kmodel, (2) the RNG k-model, (3) the Low-Reynolds Number k-model and (4) the Reynolds Stress Model (RSM). Based on convergence and computational stability criteria, they concluded that (a) the application of the standard k-model was not successful and (b) the RNG k-model was found to be the most appropriate model to characterize the flow field of the ventilated space studied. The SST k-based model [19] uses a blending function to combine the k-and the k-models; in the near-wall region the k-model is activated, while for the rest of the flow the k-model is used. By this approach, the attractive near-wall performance of the k-model is utilized without the potential errors resulting from the free stream sensitivity of that model. Stamou and Katsiris [3] applied the SST k-based model for the first time in indoor environment; they calculated air -flow velocities and temperatures in a model office room. Calculations were compared with experiments and with the results of the standard k-, the RNG k-model and the laminar model. It was concluded that (a) all the three tested turbulent models predict satisfactorily the main qualitative features of the flow and temperature fields and (b) computations with the SST k-based model showed the best agreement with measurements.

  The premises are served by a variable air volume (VAV) air conditioning system; the temperature of the incoming air could be automatically regulated. The main airsupply and exhaust devices are shown in Figure2. The air is delivered to the Arena via 256 jet nozzles, which are organized as follows: (1) 32 diffusers are mounted above the four entrances in the area of athletes in four series of 8 having a total capacity equal to 32X1,675= 53,600 m 3 /h, and (2) 224 diffusers are mounted in groups of 56 on the 4 air pipes in the roof of the Arena having a total capacity equal to 224X1,005= 225,120 m 3 /h. The angle of all the diffusers above the entrances was set equal to 28 with the horizontal; a positive angle means that the jet faces upwards. The jets from the diffusers in the roof were facing towards the interior of the Arena or towards the outer walls for the ventilation of the space of the corridors; the angles of these jets were variable ranging from -30 o to +18 o .The air is extracted to the plenum via: (1) 32 square extract grilles on walls having a total capacity equal to 36,100 m 3 /h, (2) 2,046 openings under 20 groups of the seats of a total capacity equal to 160,774 m 3 /h and (3) 57 floor openings (total capacity= 57,351 m 3 /h). The total outgoing air flow-rate is equal to 254,225 m 3 /h. In the calculations the outlet flow was set equal to the incoming flow rate.4 Application of the modelCalculations were performed for two scenarios of operation of the HVAC system, which corresponded to two values of temperatures of conditioned inlet air (a) T in =14 o C and (b) T in =16 o C.4.1 Boundary conditions and the numerical gridIn Figure3a view of (a) the real geometry and (b) the simulated geometry of the interior of the Arena are shown. Figure3depicts that the simulated model has almost the same geometrical configuration as the real Galatsi Arena.Galatsi Arena was modeled with a normal occupancy of 6,000 persons. Existing furniture, supply-exhaust devices, lights etc. were included in the model with their approximated dimensions to simulate the real environment. The supply and exhaust boundary conditions were specified to summer conditions; the outdoor air temperature was set equal to 36.5 o C. The air flow-rates and the temperatures of the air supply diffusers and jets were set equal to their actual values, while at the extract grilles and openings the pressure was set equal to the atmospheric pressure.

Figure 6 ,

 6 Figure 6, are plotted. In Table1the values (ranges and average values) of the calculated

Figure 1 .Figure 2 .Figure 3 .Figure 4 .Figure 5 .Figure 6 .Figure 7 .Figure 8 .Figure 9 .

 123456789 Figure 1. View of the Galatsi Arena (www.stadia.gr)

Figure 2 Figure 3 (

 23 Figure 2

Figure 9

 9 Figure 4

Table 2 .

 2 5 clo (i.e. 0.078 m 2 °K/W) and the metabolic rate of the spectators equal to 1.2 met (i.e. 70 W/m 2 ). The environmental variables were determined from the CFD calculations. The mean radiant temperature was calculated equal to 24.62 C and 26.11 C for T in =14 C and 16 C, respectively. The mean values of the air temperature and velocity in the areas of interest, which are shown in Table 1, were used in the calculation of the comfort indices. An EXCEL spreadsheet was formulated for the calculation of PMV and PPD indices in the areas of interest; these values are shown in The figures of Table2show that for T=14 C the thermal conditions are not

	satisfactory. The low-negative values of PMV (almost in all rows of seats PMV<-0.5)
	indicate that the air temperatures are lower than the desired values; the percentage of the
	occupants who are expected to express dissatisfaction ranges from 8.4% to 27.4%. For
	T=16 C the relatively high air velocities (up to 0.40 m/sec) are combined with the
	relatively high temperatures (up to 26.5 C); the desired range of temperatures is 24-
	26 C [25]) to create thermal conditions which are generally pleasant; PPD values show

that less than 7% of the spectators are expected to express dissatisfaction. In the area of athletes air temperatures for T=14 C are low (22-24 C), while for T=16 C are very satisfactory (24-26 C); moreover, the relatively high air drafts generally create to the athletes a feeling of freshness.

Table 2 .

 2 Calculated values of PMV and PPD

	T in	14 o C	14 o C 16 o C	16 o C
	Area-Sector	PMV PPD	PMV PPD
	Lower-Southern	-0.71	15.5	-0.05 5.1
	Middle-Southern	-0.66	14.1	+0.26 6.5
	Upper-Southern	-0.59	12.4	+0.07 5.1
	Lower-Western	-0.64	13.6	-0.22 6.0
	Middle-Western	-0.90	22.2	-0.24 6.1
	Upper-Western	-0.40	8.4	+0.19 5.7
	Lower-Eastern	-1.03	27.2	-0.04 5.0
	Middle-Eastern	-1.03	27.4	-0.10 5.2
	Upper-Eastern	-0.52	10.6	+0.16 5.5

Acknowledgments

This work is a part of a Research Project: "Evaluation and improvement of the Indoor Conditions (safety, air quality and comfort) in the Indoor Stadiums and Amphitheaters of the Olympic Games Athens-2004 with the use of Mathematical Models", 2003-2004, which was financed by the Ministry of Culture (General Secretariat for the Olympic Games 2004). The authors would like to thank the Ministry of Culture; especially the former General Secretary Prof. C. Kartalis. Thanks are also due to all the scientists, who provided information on the Galatsi Arena.