Thermal stratification in an underfloor-ventilated room: comparing code_saturne against the Zhang and Chen experiment

Réalisé via code_saturne v9.1.0 - le code CFD open-source d'EDF

At a glance

The experiment of Zhang & Chen (2006) reproduces an occupied room ventilated through the floor. Fresh air is supplied by two inlets located at floor level and extracted through an opening in the ceiling. The internal heat loads consist of four heated manikins dissipating $100$ W each and six $64$ W lamps. The authors publish vertical temperature and velocity profiles at seven positions, together with the full set of measured boundary conditions. This type of mixed-convection flow, resulting from the interaction between mechanical ventilation and buoyancy effects, exhibits strong thermal stratification and constitutes a demanding validation case for CFD applied to buildings.

The case is simulated with code_saturne v9.1 (Boussinesq approximation, k-ε, steady algorithm) and compared with an OpenFOAM v2406 computation using an identical setup. The main findings:

• On the same 2.7-million-cell mesh, code_saturne and OpenFOAM give nearly indistinguishable fields and profiles.
• The flow structure is correctly reproduced: supply jets spreading over the floor, plumes above the manikins, temperature increasing continuously from floor to ceiling. Velocities are reproduced without systematic bias (mean deviation of about 0.02 m/s). The occupied-zone temperature is however underestimated by about 0.6 K by both codes.
• By imposing the thermal load powers as fluxes on their surfaces, only convection is modelled, whereas in practice part of this heat is exchanged by radiation. This hypothesis, the most likely explanation for the cold bias, will be examined in a dedicated article.

The reference case

The Zhang and Chen experiment reproduces an occupied room ventilated through the floor, in which four heated manikins and six lamps constitute the internal heat loads. The experimental cell measures $4.9 \times 2.4 \times 4.3$ m and is equipped with an UFAD (underfloor air distribution) system. Air at $20^\circ\text{C}$ is introduced through two square inlets of $0.26$ m side length, located at floor level, and extracted through an opening in the ceiling. Each inlet delivers a flow rate of $0.0472$ m³/s, corresponding to a mean discharge velocity of about $0.70$ m/s.

The four manikins, of dimensions $0.38 \times 0.38 \times 0.90$ m, dissipate $100$ W each. The six lamps, flush-mounted in the ceiling, dissipate $64$ W each, giving a total internal heat load of $784$ W. A uniform temperature, taken from the experimental measurements, is imposed on each of the six walls. The prescribed values differ from one wall to another and range between $23.9$ and $25.6^\circ\text{C}$.

The validation relies on vertical profiles of temperature and velocity magnitude measured on seven poles, labelled V1 to V7, placed in the mid-plane of the room. Each profile comprises seven measurement points, distributed between $0.1$ and $2.2$ m above the floor. The temperature is published in dimensionless form:

$$T^* = \frac{T-T_{\text{supply}}} {T_{\text{exhaust}}-T_{\text{supply}}}.$$

In the present modelling, the powers of the manikins and lamps are converted into uniform heat fluxes applied to their surfaces, rather than imposing the measured surface temperatures. This choice corresponds to a common situation at the design stage, where the powers dissipated by the equipment are known before their operating temperatures.

Modelling

The case is simulated with two open-source CFD codes. The reference computation is carried out with code_saturne v9.1, while an independent computation with OpenFOAM v2406 and the buoyantBoussinesqSimpleFoam solver is used for cross-verification. Both simulations rely on the same mesh, the same physical properties, the same boundary conditions and the same thermal powers. The discretisation schemes and under-relaxation parameters are chosen as close as the formulations specific to each code allow.

In both models, density variations are accounted for through the Boussinesq approximation. The density is assumed constant in the conservation equations, except in the buoyancy term:

$$-\rho_0 \beta \left(T-T_0\right)\boldsymbol{g},$$

with $\rho_0 = 1.1965$ kg/m³, $\beta = 3.43 \times 10^{-3}$ K⁻¹ and $T_0 = 298.15$ K. Turbulence is modelled with the standard $k$–$\varepsilon$ model and the turbulent Prandtl number is set to $0.85$ in both codes. The near-wall treatments follow the formulations available in each of them.

Parameter	code_saturne v9.1	OpenFOAM v2406
Density	$1.1965$ kg/m³	$1.1965$ kg/m³
Dynamic viscosity	$1.83 \times 10^{-5}$ Pa·s	$1.83 \times 10^{-5}$ Pa·s
Buoyancy	Boussinesq, $\beta = 3.43 \times 10^{-3}$ K⁻¹, $T_0 = 298.15$ K	same parameters
Turbulence	standard $k$–$\varepsilon$	standard $k$–$\varepsilon$
Wall treatment	two-scale wall law	corresponding wall functions, including the Jayatilleke model for thermal transfer
Turbulent Prandtl number	$0.85$	$0.85$
Pressure-velocity coupling	steady SIMPLEC	steady SIMPLE algorithm
Heat sources	uniform fluxes of $427$ W/m² for the lamps and $66$ W/m² for the manikins	same surface fluxes
Room walls	imposed uniform temperatures, between $297.0$ and $298.7$ K	same temperatures
Supply	imposed temperature and flow rate	same conditions
Exhaust	outlet condition at the ceiling	equivalent condition

The hexa-dominant mesh is generated with snappyHexMesh. It is used directly in OpenFOAM, then converted to the MED format for import into code_saturne using meshlane. The geometry, connectivity and number of cells are thus preserved between the two computations.

In code_saturne, the Boussinesq approximation is activated through the algo_density_variation parameter. The thermal expansion coefficient and the turbulent Prandtl number are defined through short user functions compiled at run launch.

A preliminary mesh sensitivity study was carried out with two grids of $1.1$ and $2.7$ million cells respectively. Since the differences observed on the temperature and velocity profiles at the measurement positions are small, only the results obtained with the $2.7$-million-cell mesh are presented in what follows.

Results and discussion

The following figure shows the temperature fields (top) and the velocity fields with streamlines (bottom) in the mid-plane $z = 0$, for code_saturne (left) and OpenFOAM (right).

These fields illustrate the principle of underfloor ventilation. The two dark columns in the temperature field correspond to the fresh-air jets issued from the supply inlets, located at $x = \pm 1.5$ m in this plane. Introduced at 0.70 m/s, these jets retain enough momentum to rise to about 1.5 m, a height at which their negative buoyancy progressively opposes their upward motion. They reach a height of about $1.5$ m before deflecting and feeding the lower part of the room.

The result is cooler air at floor level, then a temperature that increases continuously with height, the top of the room being fed with warm air by the plumes of the heat sources. The warmest film, visible under the ceiling, comes from the lamps, two strips of which cross the plane shown. The velocity field confirms this organisation: outside the jets, velocities remain below 0.1 m/s, and the streamlines show the upper layer converging towards the exhaust opening at the centre of the ceiling. The two codes produce fields that are difficult to tell apart.

The mid-plane, however, contains no heat source: the manikins are located at $z = \pm 1.5$ m. The next figure shows the same fields in the plane $z = 1.5$ m, which crosses two of the four manikins.

Above each manikin, the heated air rises in a plume clearly visible in the temperature field. The velocity field shows that this plume reaches 0.25 to 0.30 m/s, more than the supply jets at the same height. The streamlines converge towards the plume from the ambient air: as it rises, the plume draws in and entrains the surrounding air, so that the flow rate it carries increases with height. Once under the ceiling, this flow discharges into the upper part of the room and feeds it continuously. The observed stratification results from this balance: the plumes carry the heated air upwards, and the ceiling exhaust removes part of it. The two codes reproduce this mechanism in a very similar way.

The quantitative comparison with the measurements is based on the vertical profiles at the seven poles.

On all seven poles, the computed profiles show a temperature that increases continuously from floor to ceiling, with no homogeneous zone. By 0.1 m above the floor, the air supplied at 20 °C has already warmed rapidly in contact with the floor, maintained at nearly 24 °C, and through mixing with the ambient air around the jets. Higher up, the temperature increases steadily, by about 1.1 to 1.7 K per metre, and the profile steepens slightly towards the ceiling. The value $T^{*} = 1$, which by definition corresponds to the temperature of the extracted air, is only reached at about 2.2 m, just below the ceiling where the exhaust opening is located. The stratification thus occupies the whole height of the room rather than organising itself into two layers separated by a sharp interface, and the experimental points show the same continuous behaviour.

Over the last few centimetres below the ceiling, poles V1, V2, V3, V5, V6 and V7 exhibit a sharp temperature rise, well beyond $T^{*} = 1$: these poles end in the warm-air film maintained by the lamps. The rise is much weaker on pole V4, the farthest from the lamp strips, whose top is moreover located 25 cm from the exhaust opening, which locally draws in this warm film.

The code_saturne and OpenFOAM curves are almost superimposed on all poles, including close to the ceiling, and the differences between codes remain far smaller than the deviation from the measurements. The latter therefore stems from the modelling rather than from the choice of solver. The deviation from the measurements is visible in the intermediate zone: between 0.5 and 1.7 m, the experimental points are systematically warmer than both computations.

The velocity profiles can be read in three zones. Near the floor, the computations show a thin current spreading from the supply inlets, with a maximum of 0.10 to 0.13 m/s located within the first two centimetres, below the first measurement point. At 0.1 m, the computed velocities (0.05 to 0.10 m/s) remain higher than the measured ones (0.02 to 0.07 m/s). Pole V4, equidistant from the two inlets, is an exception: the floor current reaches it weakened, and the computed maximum, about 0.08 m/s, is located around 0.3 m. In the intermediate zone, between 0.9 and 1.9 m, the central poles V3 to V5 cross nearly quiescent air (less than 0.02 m/s), while poles V1, V2, V6 and V7 retain a bulk motion of 0.05 to 0.08 m/s; the measured velocities there are of the same order. Near the ceiling, finally, the computed velocity rises again on the central poles, strongly on V4 whose top lies 25 cm from the exhaust opening: the computations give about 0.14 m/s at 2.2 m against a measured 0.06 m/s, the largest local deviation of the comparison, in the suction zone where the velocity varies strongly over a few centimetres. Over all points, the bias remains nearly zero (MBE below 0.005 m/s) and the curves of the two codes are once again almost superimposed.

The agreement with the measurements is quantified by three metrics computed over the $N$ experimental points of the seven measurement poles ($N = 49$ measurement points), where $\phi$ denotes the compared quantity ($T^{*}$ or $|V|$):

$$\mathrm{MAE} = \frac{1}{N}\sum_{i}\left|\phi_i^{\mathrm{CFD}} - \phi_i^{\mathrm{exp}}\right|, \qquad \mathrm{MBE} = \frac{1}{N}\sum_{i}\left(\phi_i^{\mathrm{CFD}} - \phi_i^{\mathrm{exp}}\right), \qquad \mathrm{RMSE} = \sqrt{\frac{1}{N}\sum_{i}\left(\phi_i^{\mathrm{CFD}} - \phi_i^{\mathrm{exp}}\right)^{2}}.$$

The MAE (Mean Absolute Error) measures the average deviation, the signed bias MBE (Mean Bias Error) indicates whether the computation over- or underestimates on average, and the RMSE (Root Mean Square Error) penalises the largest deviations. The computed profile values are interpolated at the heights of the measurement points.

Metric	code_saturne	OpenFOAM
MAE on $T^{*}$	0.121	0.109
MBE on $T^{*}$	−0.105	−0.088
RMSE on $T^{*}$	0.136	0.123
MAE on $\lvert V \rvert$ (m/s)	0.024	0.022
MBE on $\lvert V \rvert$ (m/s)	+0.003	+0.002
RMSE on $\lvert V \rvert$ (m/s)	0.030	0.029

The metrics table sharpens the observation made on the temperature profiles: the MBE of −0.105 accounts for nearly all of the MAE (0.121), which means that the error is almost entirely a systematic cold bias, not scatter. Between 0.5 and 1.7 m, that is, in the occupied zone, the computation is too cold by about 0.6 K on average. This bias is similar for both codes, which points to a missing physical ingredient rather than a numerical artefact.

By imposing the load powers as heat fluxes on their surfaces, with no radiation model, the energy is transferred to the air by convection only. The heat dissipated by the lamps thus mainly feeds rising plumes directed towards the upper part of the room, with no direct radiative transfer towards the floor or the occupants. In the experiment, a fraction of this energy is most likely exchanged by radiation with the surfaces located in the lower part of the room. The absence of this mechanism in both computations is therefore a plausible explanation for the cold bias observed in the occupied zone, all the more so as tests carried out with another turbulence model and another near-wall treatment did not significantly reduce it. This hypothesis will be examined in an article devoted to the explicit modelling of radiative exchanges.

Conclusion

The underfloor-ventilation experiment of Zhang and Chen was simulated with two open-source CFD codes, code_saturne v9.1 and OpenFOAM v2406, with identical mesh, physical properties and boundary conditions. Three lessons emerge.

The two codes give nearly indistinguishable results, on the fields as well as on the profiles at the seven measurement poles. On this mixed-convection case strongly influenced by buoyancy effects, code_saturne constitutes a credible alternative to OpenFOAM, and the toolchain used, from the snappyHexMesh mesh converted to the MED format, thanks to meshlane to the Boussinesq setup with a steady algorithm, is entirely open-source.

The physics of underfloor ventilation is correctly reproduced: fresh-air jets spreading over the floor, plumes above the manikins, temperature increasing continuously from floor to ceiling. Velocities are reproduced without systematic bias, with localised deviations in the floor current and in the suction zone of the exhaust.

The occupied-zone temperature, however, remains underestimated by about 0.6 K, in a similar way for both codes. The closeness of the two solutions suggests that the main deviation is associated with a shared modelling assumption rather than with the implementation specific to either code. The absence of radiative transfer between the heat loads and the lower surfaces of the room is its most likely explanation. Taking it into account will be the subject of a forthcoming article.

References

• Z. Zhang, Q. Chen, Experimental measurements and numerical simulations of particle transport and distribution in ventilated rooms, Atmospheric Environment 40 (2006) 3396–3408. [doi:10.1016/j.atmosenv.2006.01.014]
• code_saturne v9.1, EDF. [code-saturne.org]
• OpenFOAM v2406, OpenCFD Ltd. [openfoam.com]
• meshlane, go to https://github.com/simvia-tech/meshlane