Building and Environment 123 (2017) 37e49

Contents lists available at ScienceDirect

Building and Environment

journal homepage: www.elsevier.com/locate/buildenv

Transient three-dimensional CFD modelling of ceiling fans

Francesco Babich a, *, Malcolm Cook a, Dennis Loveday a, Rajan Rawal b, Yash Shukla b
a
School of Civil and Building Engineering, Loughborough University, Loughborough, UK
b
CEPT University, Ahmedabad, India

a r t i c l e i n f o a b s t r a c t

Article history: Ceiling fans have been used for decades as a means of providing thermal comfort in tropical countries
Received 15 April 2017 such as India. However, recent years have witnessed a significant increase in the use of air conditioning
Received in revised form as a means to achieve comfort, and therefore in the total energy consumption and related CO2 emissions.
15 June 2017
Ceiling fans are still viable options to limit use of air conditioners or in combination with air conditioners
Accepted 21 June 2017
Available online 23 June 2017
without compromising on thermal comfort and still achieving energy savings. Ceiling fans generate non-
uniform velocity profiles, and therefore relatively non-uniform thermal environments, whose charac-
teristics may be tough to analyse with simple modelling methods. This issue can be investigated using
Keywords:
Ceiling fan
CFD. However, to date, there are few works on ceiling fans, CFD and thermal comfort. More accurate
Thermal comfort models are therefore required to predict their performance. The research presented in this paper aimed
CFD validation to develop and validate a three-dimensional transient implicit CFD model of a typical ceiling fan available
Turbulence modelling in India by comparing simulation results obtained using different URANS turbulence models with
India measured data collected in controlled environment. The results highlight that this ceiling fan model is
Environmental chamber able to replicate the predominant characteristics of the air flow generated by the fan such as the
meandering plume and the local fine free shear layers. The best results are achieved when the SST k-u
turbulence model is used, with 83% of the simulated values being within the error bars of the respective
measured value.
© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license
(http://creativecommons.org/licenses/by/4.0/).

1. Introduction energy-intensive air conditioning is limited only to periods of

extremely hot weather, then overall Indian energy consumption
Ceiling fans have been used for decades as a means of providing and related CO2 emissions will significantly increase, leading to
thermal comfort in tropical climates. In Indian residences, ceiling severe implications for the global climate and also challenging the
fans are as common as electric light bulbs, being present in almost reliability of the Indian electricity grid.
every habitable space. They are part of most of Indian residences, Air conditioning usually provides uniform thermal environ-
and they are widely used in both old and more recent buildings mental conditions, therefore designers can predict with confidence
[1,2]. In the event of growing economy and rising percentage of its performance using traditional thermal comfort models [6]. On
population which can afford purchase and operation of air condi- the other hand, ceiling fans generate non-uniform velocity profiles,
tioner for higher want of thermal comfort India has experienced and therefore relatively non-uniform thermal environments. This
rise in sales of air conditioners [3]. Cities such as Delhi, Mumbai and does not imply a lower thermal comfort for the occupants, but
Kolkata, having Cooling degree days in range of 3000 and 3500, are could possibly lead to more thermally comfortable environments
replacing fans with air conditioners despite possibility of use of fan with lower energy cost due to air velocity, as research on allies-
during certain part of year. It is important to note that most air thesia suggests [7]. The positive effect of air movement on thermal
conditioners have low-efficiency and use high-GWP refrigerants comfort in warm and hot conditions has also been included in in-
[4]. The electricity demand for space cooling comprises up to 60% of ternational standard such as the current version of ASHRAE 55 [8]:
summer peak load in large cities such as Delhi [5]. Unless the use of for operative temperature above 25.5
C, the air speed limits have
been raised to 1.2 m/s with occupant control and 0.8 m/s without
occupant control. However, accurate models of the ceiling fans are
required to predict their ability to generate air velocity and their
* Corresponding author. effect on occupant comfort using newer and more advanced
E-mail address: f.babich@lboro.ac.uk (F. Babich).

http://dx.doi.org/10.1016/j.buildenv.2017.06.039
0360-1323/© 2017 The Authors. Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/).
38 F. Babich et al. / Building and Environment 123 (2017) 37e49

thermal comfort models such as the IESD-Fiala model [9,10] set-up has then been modelled using a commercial CFD program,
coupled with a commercial CFD software [11]. Traditional thermal ANSYS CFX [26]. In this model, the fan is modelled as a momentum
comfort models, namely the PMV-PPD model developed by Fanger source that is specified by radial components and applied to a thin
[6] and the more recent adaptive approach [12,13], are indeed not cylindrical sub-domain of the CFD model which has the same
suitable to investigate thermal comfort in non-uniform thermal diameter as the fan used in chamber experiment. Measured values
environments. The former is applicable only to uniform steady- and computer-generated predictions have been compared ana-
state conditions, while the adaptive model uses only the outdoor lysed, and the reasons for any deviation discussed.
temperature as input parameter, and therefore it cannot be used to
study the effect of a specific device such as a fan on thermal 2.1. Experimental set-up
comfort.
The IESD-Fiala model is a model of human thermal physiology Environmental chambers have been extensively used to collect
and comfort that predicts passive [9] and active [10] reactions of a data for validating CFD models [27]. Within an environmental
human body in response to certain environmental conditions. It has chamber, most of the variables can be controlled, and the state of
been successfully coupled in real-time with a commercial CFD the uncontrolled variables ought to be accurately determined. Thus,
software [11], where real-time means that at each time step of a the number of required assumptions due to the lack of measure-
transient simulation there is an exchange of information between ments or better information, and therefore the number of potential
the two parts of the system. Thus, this coupled system can sources of uncertainties, are minimised.
accommodate transient and asymmetrical environments [14], and For this research, an environmental chamber located at CEPT
can therefore be used to model the usage of ceiling fans, for University, Ahmedabad, in India was used to generate validation
instance identifying the best location for their installation. In the data (see Fig. 1). Chamber is located in basement with only one side
literature, there are other advanced models such as the model exposed to outdoor environment. Internal (thickness 200 mm) and
developed at UCB by Huizenga [15], but their capabilities are lower. external (thickness 350 mm) walls have U-value of 0.29 and
Like the IESD-Fiala model, the UCB model is based on the multi- 0.28 W/m2K, respectively. Roof, false ceiling, and floor have U-value
node model developed by Stolwijk for aerospace applications of 0.29, 0.55 and 2.80 W/m2K, respectively. External window facing
[16,17], and it has been used together with CFD software [18e20]. West direction has U-value of 1.79 W/m2K. At no point in time
However, in these studies there is only a manual coupling, with no direct solar radiation hit this window as it is well protected outside.
real-time automated procedure. The mechanical ventilation system comprises four supply grilles
To date, there are few existing reports on the use of CFD to and two return grilles, which are located in the ceiling. The
model ceiling fans. Bassiouny et al. [21] implemented simple 2D chamber ventilation system was turned off during any experiment,
models. Thus, the truly three-dimensional behaviour of air flow and the diffusers were sealed during both the second and third
generated by a ceiling fan is not captured. Momoi et al. [22,23] repetition of the measurements. This did not produce any signifi-
developed too complex modelling approaches that are limited in cant difference in the measured values, but allowed to more
their application because many measurements are needed due to accurately set the boundary conditions in CFD. Moreover, since the
the required input data. Adeeb et al. [24] focused on the effect of the objective was to obtain data of velocity profiles under uniform and
different numbers of blades, but their work did not provide suffi- stable environmental conditions, the experiment was conducted
cient information to build a 3D model of a ceiling fan in CFD to be once the chamber stabilized at constant air and surface tempera-
used for thermal comfort studies within the Fiala-IESD and CFD tures. The stability of chamber was determined based on less than
coupled system. To do so, the most important thing is accurate 0.1
C change in air and surface temperature for more than 12 h. A
modelling of the environment around the occupants, but keeping typical 3-blade 1200 mm sweep Indian ceiling fan [28] with a 4-
the CFD model as simple as possible for two reasons: reducing the step regulator was installed on the ceiling of the chamber. The
time and computational power required to achieve a converged diameter of central part is 190 mm and the length of the blade is
solution, and avoiding potential sources of error due to the use of an 500 mm. Other than measuring equipment, no other objects were
unnecessary complex model. For instance, the use of a moving present in the room.
mesh used by Zhu et al. [25] to study the effect of the ceiling fans on Air speed measurements were recorded in three horizontal
air mixing and UR-UVGI disinfection efficacy [25] would not planes in the chamber: 0.1 m, 0.7 m and 1.3 m above the floor. In
guarantee more accurate thermal comfort predictions, but would each plane, 12 measurements were recorded. These were located
certainly increase the computational time, especially when tran- (see Fig. 2) below the centre of the fan (“centre”), below the
sient and steady-state simulations are used. perimeter of the fan (“north, east, south, west”), and on a radius at
Thus, the research presented in this paper aimed to develop and increasing distance from the centre of the fan. For instance, “r 800”
validate a three-dimensional transient CFD model of a typical means 800 mm away from the fan centre. Thus, in total, mea-
ceiling fan predominantly used in India by comparing simulation surements were taken at 36 points (see Table 1).
results and measured data. This model combines accuracy with Measurements were recorded simultaneously at these three
efficient computation, and can be used for accurate thermal com- heights in each of the 12 horizontal locations using three air speed
fort studies. probes (see Table 2 and Fig. 1). After the recording was completed in
This paper is divided in to five parts: (1) Introduction (2) one location, the measurement equipment was moved to the next
Methods, this has two sub sections one deals with experiment and spot, but the new recording period always started a few minutes
second modelling (3) Results from experiment and modelling (4) after having moved the equipment. The logging duration per
Discussion (5) Summary, including key findings and limitation of location was 300 s and the logging period 1 s. It is important to
work. point out that both types of probe used for measuring air speed are
omnidirectional, which means that only air speed is recorded,
2. Methods rather than the three components of air velocity separately. Their
measuring range is 0e10 m/s, in 20 to þ70
C range. Instanta-
Detailed measurements of the air movement generated by a neous measurements of room air temperature and relative hu-
typical Indian ceiling fan have been collected in an environmental midity, room surface temperatures, and fan rotational speed have
chamber in controlled environment, and the same experimental also been taken (see Table 3). In this study, the highest available fan
 F. Babich et al. / Building and Environment 123 (2017) 37e49 39

Fig. 1. Environmental chamber at CEPT University, Ahmedabad, India (on the left), and panoramic view of the air speed measurement equipment (on the right).

Table 3
Instantaneous measurements of the boundary conditions.

Parameter Value

Air temperature 31.9
C

Relative humidity 66.6%
Mean radian temperature 31.8
C
Fan rotational speed 290 RPM

A complete set of measurements was taken three times to

ensure the quality of the data and repeatability of the experiment.

2.2. Modelling approach and assumptions

The identical configuration used during the experiments was

then recreated within the CFD program. In this study, the chosen
Fig. 2. Plan view of the environmental chamber with the measurement locations. software is ANSYS CFX [26] primarily because the coupling be-
tween the IESD-Fiala model and CFD was completed using ANSYS
CFX due to its customization features which facilitate the connec-
rotational speed setting has been chosen based on an on-going field tion with other software [11]. Thus, this facilitates the use of this
study in India and in the UK on air movement in domestic buildings ceiling fan modelling approach in thermal comfort applications.
[29]. However, the modelling approach presented in this paper could be

Table 1
Measurement points.

Height [mm] Horizontal location

Centre North East South West r200 r500 r600 r800 r1200 r1700 r2000

1300 Centre 1300 North 1300 East South West r200 r500 r600 r800 r1200 r1700 r2000
1300 1300 1300 1300 1300 1300 1300 1300 1300 1300
700 Centre 700 North East South West r200 r500 700 r600 r800 r1200 r1700 r2000
700 700 700 700 700 700 700 700 700 700
100 Centre 100 North East South West r200 r500 100 r600 r800 r1200 r1700 r2000
100 100 100 100 100 100 100 100 100 100

Table 2
Measurement equipment characteristics.

Instrument Parameter measured Accuracy

Testo 480 Climate e hot bulb probe type 1 Air speed (at 0.1 m and 1.3 m height) ±(0.03 m/s þ 5% of reading)
Testo 480 Climate e hot bulb probe type 2 Air speed (at 0.7 m height) ±(0.03 m/s þ 4% of reading)
TSI Veloci-calc Multifunction Ventilation meter 9565 Dry-Bulb Temperature ±0.3
C, ±3% RH
and RH
Fluke-561 Infrared Thermometer Surface temperature ±1% or ±1
C whichever is greater
Lutron DT-2236 Digital contact tachometer RPM ±0.05% or 1 digit
40 F. Babich et al. / Building and Environment 123 (2017) 37e49

implemented in most CFD programs.

The geometry (see Fig. 3) was created in ICEM CFD [30], a gen-
eral purpose geometry and mesh generator for CFD. Air diffusers,
windows and door were not included since they were all closed
during the experimental work, and because only isothermal sim-
ulations were conducted to reduce the required computational
power. This has been done because the surface temperatures and
the air temperature were similar, being 31.8
C and 31.9
C,
respectively (see Table 3). Thus, the magnitude of any buoyancy-
driven air flow would have been negligible compared with the air
speed generated by the fan.
The ceiling fan has been modelled as a momentum source.
Modelling the actual blades would require very detailed informa-
tion about their geometry, and would lead to a much higher
number of mesh elements and to the use of a moving mesh. Both
these features would significantly increase the demand for
computational power, and the possible sources of uncertainty,
without guaranteeing better results, but limiting the usability and
applicability of the model. Thus, the fan has been modelled as a ring
Fig. 4. Detail of the fan simplified geometry (ring in yellow, central element in purple).
with the same diameter, 120 cm, and distance from the ceiling,
(For interpretation of the references to colour in this figure legend, the reader is
30 cm, as the actual fan, and with a central cylindrical solid referred to the web version of this article.)
element, since in a real ceiling fan no air emanates from the centre
[31]. The momentum source has been applied to this ring (see
Figs. 3 and 4, in which the ring is highlighted in yellow). fan model as a component that can be used with advanced tran-
In this study, an unstructured mesh was used because it is more sient thermal comfort models, running transient simulations from
flexible and it can be used with any complex shape such as a human the very beginning makes its applicability easier and more reliable.
body, easing the wider applicability of the proposed modelling A second order backward Euler scheme, and high resolution
approach. Similarly to previous studies [14], both the surface and advection scheme and turbulence numerics have been used. The
the volume mesh were initially generated in ICEM CFD using the double precision option was also used to enhance accuracy and
Octree algorithm, then the volume mesh was deleted, and finally robustness of the results. Convergence criteria have been set equal
regenerated using the Delaunay algorithm (see Fig. 5). This pro- to 1e-05 for the RMS residuals and 0.01 for the conservation target.
cedure ensures the robustness of the mesh and a smoother tran- Since there are no inlets and outlets in this model, the conservation
sition from smaller elements to larger ones. Ten prism layers were target does not have any influence. Furthermore, the sub-domain
added adjacent to the walls to accurately model the boundary layer used to model the fan simply acts as a momentum source, but
near surfaces. cannot generate any mass imbalance as there are no physical bar-
Transient simulations have been performed to better model the riers (real surfaces) between this sub-domain and the remaining
real behaviour of the ceiling fan. Moreover, considering this ceiling part of the room, and there is no mass generation or dissipation
within this sub-domain. Moreover, an adaptive time step as a
function of RMS Courant number was chosen, with the limit for the
RMS Courant number set equal to 1. Since it is an implicit code, CFX
does not require the Courant number to be small for stability, but
this still leads to more accurate results. The initial time step was
0.1s, with the maximum and minimum values set equal to 1s and
0.01s, respectively.
The momentum source that simulates the actual fan has been
applied to the subdomain using cylindrical components: axial
component 55 kgm2s2, radial component 0 kgm2s2 and theta
component 8 kgm2s2. The axial component pushes air down-
wards while the theta component generates rotational movement.
The sign of each component depends on the orientation of the
defined local system of coordinates (see Fig. 4), while the absolute
figures have been defined with an iterative process by comparing
simulation results with measured values and by comparing the
qualitative features of the flow field with qualitative previous
studies. A study by Jain [31] provided with very useful qualitative
information about the flow field generated by a ceiling fan. Using
smoke from thick incense sticks, the flow field created by the
ceiling fan was visualised, identifying the main regions of the flow
and its key behaviours such as the swirling movement.
Due to the intrinsic turbulent nature of the air flow generated by
the fan, the correct choice of the turbulence model is fundamental.
In this study, four of the most widely used Unsteady Reynolds
Averaged Navier-Stokes (URANS) eddy-viscosity turbulence models
have been fully tested: the SST (Shear Stress Transport) k-u
Fig. 3. CFD model of the environmental chamber.
 F. Babich et al. / Building and Environment 123 (2017) 37e49 41

Fig. 5. Volume mesh generated using Delaunay algorithm.

[32e36], the RNG (Re-Normalisation Group) k-ε [37], the standard (see the yellow marks in Fig. 3) and considered to be in agreement
k-ε [38] and the standard k-u [39e42]. Two Reynolds-stress models when the respective error bars overlapped. For the simulated
have also been considered, namely the SSG (Speziale-Sarkar-Gat- values, the upper and lower limits of the error bars are the mean
ski) [43] and the BSL (Menter Baseline Model) [44], for their po- over the simulated period plus and minus the discretization error,
tential higher accuracy for simulating swirling flows and boundary respectively. For the measurements, since the experiment has been
layers [45]. repeated three times, the upper limit of the error bars is the highest
Before running the main simulations, a mesh sensitivity analysis of the three measurements plus the measurement error, while the
has been completed, with the purpose of estimating the dis- lower limit of the error bars is the lowest of the three measure-
cretization error and calculating the numerical uncertainty in the ments minus the instrument error.
fine-grid solution, and thus identifying the most suitable mesh for Statistics such as root mean square error (RMSE) and the mean
the main simulations. To do so, three meshes of varying density absolute error (MAE) have not been employed in this research for
have been used (see number of elements ni in Table 4), a repre- two key reasons. Firstly, using RMSE and MAE with n ¼ 36 might be
senting dimension hi to be calculated, and a significant simulated misleading as they would spread the error across all points. How-
value 4i has to be used to quantify this uncertainty [46]. In this ever, in this study, it is more important to evaluate the accuracy of
research, the air speed at a point below the centre of the fan at 1.5 m the model in the different regions (centre, perimeter, along a
above the floor has been chosen. The mesh chosen for the main radius), separately. Secondly, all statistics are less useful when there
simulations comprises 1,996,278 elements and a numerical un- are only a limited number of error samples [47]. Thus, in this
certainty equal to 15.80%. research, presenting the values of the errors themselves (as error
bars) is more appropriate.
2.3. Model validation - criteria for agreement
3. Results
Measured and simulated values were compared in 36 points
In this section, qualitative and quantitative simulation results
Table 4
are presented. The former are compared with previous qualitative
Discretization error calculation [46]. research, while the latter are compared with the measurements
taken in this study.
Parameter Value

n1 1,996,278
n2 678,834 3.1. Qualitative analysis of the air flow generated by the ceiling fan
n3 76,539
h1 0.03 Previous research [31] investigated the major flow regions
h2 0.05
observed in a room with an operating ceiling fan, identifying eight
h3 0.09
r21 1.43 regions (Fig. 6). The comparison between previous research and the
r32 2.07 developed CFD model is summarised in Table 5. Region one is the
41 1.56 area below the fan, and, in good agreement with previous research,
42 1.53 in the CFD model presented in this paper air speed reaches the
43 1.50
ε21 0.0298342
highest values in this region (Fig. 7). The diameter of downward
ε32 0.029464 flow immediately below the fan is smaller than the diameter of the
s 1 fan, the flow then starts to diverge, and there is a significant
iteration 0.00 swirling component (see Fig. 8). In the qualitative study [31], the
p 0.8947256
half cone angle is about 10
, while in the CFD results this varies
421
ext
1.6
depending on which time-step is analysed. This is due to the fact
e21 0.0191239
a
0.0479778
that in this region the flow is highly turbulent, and therefore the
e21
ext
21 15.80% angle of the meandering plume and the characteristic of the local
GCIfine
fine free shear layers vary over time (see Fig. 9).
42 F. Babich et al. / Building and Environment 123 (2017) 37e49

Similarly to previous research [31], in region two, the air speed

near walls and ceiling is very low, air is moving upwards near the
walls, and as soon as it is again in proximity of the fan, in region
three, its speed increases and it starts to redevelop a swirling
component. Moreover, the CFD model accurately predicts the small
region below the ceiling fan where air is not driven downwards due
to the blockage caused by the motor of the fan. In agreement with
the qualitative study [31], the wide region (number four) between
the surfaces of the room and region one is characterised by very low
air speed and negligible effectiveness of the fan. Thus, the CFD
model is able to accurately replicate all the occupied regions of the
room.
The major differences between the CFD results and previous
qualitative studies are located in the two small regions near the
ends of the blades, where the local air recirculation is under-
estimated by the CFD model. This is due to the fact that the actual
blades are not modelled, but it is not a significant limitation of the
model in this case since the aim is to model the air movement
generated by the ceiling fan in the room with particular interest in
the occupied zone.
Fig. 6. Flow regions identified by Jain [31].

3.2. Quantitative validation of the CFD model

Table 5
Comparison between previous qualitative research and the developed CFD model.
In this study, four turbulence models have been fully tested and
Region Agreement Key characteristics the CFD results significantly vary depending on the chosen model
1 Good Highest higher speed, significant swirling component, (see Table 6).
divergent flow Taking into account the uncertainties in measurement and the
2 Good Very low air speed near walls (moving upward) and ceiling
discretization error in CFD, measured and simulated values are in
3 Good Increasing air speed and development of swirling
component excellent agreement when the SST k-u was used, with 83% of the
4 Good Very low air speed and negligible effectiveness of the fan simulated values being within the error bars of the respective
5 Weak Local air recirculation underestimated by the CFD model measured value. There is an excellent agreement in the three points
6 Good Low air speed, recirculation area below the centre of the fan (see Fig. 10, three “centre” points) and in
7 Good Air not driven downwards due to the blockage caused by
the motor of the fan
the 21 points located on a radius at increasing distance (see Fig. 10,
8 Weak Local air recirculation underestimated by the CFD model points from “r200” rightwards). Only in six perimeter points (see
Fig. 14), namely “Est 700” and “Est 1300”, “North 700” and “North
1300”, and “West 700” and “West 1300”, the respective error bars
do not overlap. In three of them, namely “North 700” and “North

Fig. 7. Air flow generated by the ceiling fan (SST k-u turbulence model) e side view.
 F. Babich et al. / Building and Environment 123 (2017) 37e49 43

Fig. 8. The swirling air movement generated by the ceiling fan 10 cm below the blades (SST k-u turbulence model) e plan view.

1300”, and “West 700”, the difference between measured and [14]. This can be explained by analysing the nature of the four eddy-
simulated values is wide, while it is significantly smaller in the viscosity turbulence models in relation to the application presented
others. This is likely to be due to the fact that these six points are in in this paper.
the most turbulent and critical region of the flow, where air speed is Turbulence occurs when the inertia forces in the fluid become
characterised by rapid spatial and temporal variations (see Fig. 9). significant compared to the viscous forces. The Re is defined as the
The other three turbulence models did not reach the same high ratio between inertia and viscous forces, and therefore the higher
level of agreement with the measured values (see Table 6). In the the Re, the more turbulent the fluid flow becomes. Navier-Stokes
three points below the centre of the fan, the RNG k-ε (see Fig. 11) equations could be used to describe any fluid flow, regardless of
and the k-u (see Fig. 12) turbulence models are as good as the SST k- its turbulence level, but the required computational power to fully
u, while the k-ε fails to predict the air speed at 0.1 m above the floor, solve a given problem significantly exceeds the power that is
being 0.1 m/s the difference between the error bars (see Fig. 13). In currently available. For this reason, several turbulence models have
the 12 perimeter points, 50% or less of the simulated values are been developed which simplify the original unsteady Navier-Stokes
within the error bars of the respective measured value (see Table 6, equations by the introduction of time-averaged quantities and a
and Figs. 15e17). On the radial points, k-u and k-ε models have turbulent viscosity to produce the Reynolds Averaged Navier-
limited capability of predicting the correct air speed, the percent- Stokes (RANS) equations. In these time-averaged equations, the
age of agreement with the measurements being below 50%, and fluctuating quantities, called Reynolds turbulent stresses, are taken
there is no clear pattern that can explain why, at certain points, the into account by adding some extra terms in order to achieve
prediction is correct and in others it is not. When the RNG k-ε closure. Additional equations are then required to balance the
model is used, the agreement is excellent between r200 and r800 at added unknowns. The equations used to close the system and their
any given height (see Fig. 11), while it weakens farther away from numerical coefficients define the type of turbulence model [48].
the fan, where air speed is lower. The four turbulence models tested in this research are two-
equation eddy-viscosity models, which means that the Reynolds
4. Discussion stresses are assumed to be proportional to mean velocity gradients.
In the k-ε models, the two additional transport equations are used
4.1. The importance of the choice of the turbulence model to describe the kinetic energy k and the rate of dissipation of k per
unit mass ε, respectively:
The air flow generated by a ceiling fan is highly turbulent. As
 
shown in the previous section, it is characterised by elevated air vðrkÞ m
speed variations, and therefore by significant Reynolds number (Re) þ divðrkUÞ ¼ div t gradðkÞ þ 2mt Sij ,Sij  rε (1)
vt sk
variations. For this reason, choosing the most appropriate turbu-
lence model is essential in order to obtain accurate results.
 
The SST k-u turbulence model produced the most accurate and vðrεÞ m ε ε2
þ divðrεUÞ ¼ div t gradðεÞ þ C1ε 2mt Sij ,Sij  C2ε r
realistic results, being superior to the other three eddy-viscosity vt sε k k
turbulence model considering both qualitative and quantitative
(2)
results. This is in agreement with findings from previous research
which used the IESD-Fiala model to investigate the effect of air The equations contain five coefficients, whose value in the
movement generated by desktop fans on human thermal comfort standard k-ε model are constant [38]: Cm ¼ 0.09, sk ¼ 1.00, sε ¼ 1.30,
44 F. Babich et al. / Building and Environment 123 (2017) 37e49

Fig. 9. Air flow field generated by the ceiling fan at different time steps (SST k-u turbulence model).
 F. Babich et al. / Building and Environment 123 (2017) 37e49 45

Table 6 Wilcox [39e42], the turbulence frequency u ¼ ε/k is used as the

Comparison between measured and simulated air speed values. second variable, and the two additional equations are:
Location Number of points Agreement measurements - simulations
  
SST k-u k-u vðrkÞ m *
RNG k-ε k-ε
þ divðrkUÞ ¼ div m þ t gradðkÞ þ Pk  b rku (3)
Centre 3 100% 100% 100% 67%
vt sk
Perimeter 12 50% 50% 42% 42%
Radius 21 100% 71% 48% 48%    
vðruÞ m
þ divðruUÞ ¼ div m þ t gradðuÞ þ g1 2rSij ,Sij
Total 36 83% 67% 50% 47%
vt su
CPU time 6 h 35 m 6 h 00 m 3 h 34 m 2 h 10 m 
2 vUi
 ru dij  b1 ru2 (4)
3 vxi

C1ε ¼ 1.44, C2ε ¼ 1.92. In the RNG k-ε model [37], the C1ε is not fixed, The models constants are: sk ¼ 2.00, su ¼ 2.00, g1 ¼ 0.553,
and it represent a strain-dependent correction term introduced to b1 ¼ 0.075, b* ¼ 0.09. Subsequently, an improved k-u model was
improve the performance at low Re. In the k-u model developed by suggested, namely the SST k-u model [32e35], which is a combi-
nation of the standard keε model used in the fully turbulent region

Fig. 10. Measurements and CFD results comparison at increasing distance from the axis of the ceiling fan (SST k-u turbulence model).

Fig. 11. Measurements and CFD results comparison at increasing distance from the axis of the ceiling fan (RNG k-ε turbulence model).
46 F. Babich et al. / Building and Environment 123 (2017) 37e49

Fig. 12. Measurements and CFD results comparison at increasing distance from the axis of the ceiling fan (k-u turbulence model).

Fig. 13. Measurements and CFD results comparison at increasing distance from the axis of the ceiling fan (k-ε turbulence model).

far from the wall where Re is likely to be higher, and a keu model in Considering the three-dimensional transient ceiling fan model
the near-wall region. In the SST k-u model, the u equation changes proposed in this paper, it is now easier to understand why the SST
from equation (4) to: k-u model produced the best results followed by the RNG k-ε
model, while the other two turbulence models generated less ac-
    curate results. By combining the strengths of k-ε and k-u models,
vðruÞ mt
þ divðruUÞ ¼ div mþ gradðuÞ þ g2 2rSij ,Sij the SST k-u can deal with higher and lower Re, and it is able to
vt su;1
 accurately model the boundary layers and the flow separation
2 vU r vk vu under adverse pressure gradient conditions thanks to the extra
 ru i dij  b2 ru2 þ 2
3 vxi su;2 u vxk vxk cross-diffusion term on the right hand side of equation (5). The CFD
(5) results indeed demonstrate that this turbulence model is more
accurate in all considered locations, namely below the centre of the
There is an extra source term on the right hand side of the fan, on the perimeter, and on a radius at growing distance for any
equation, and the numerical coefficients are [36]: sk ¼ 1.00, given height and therefore air speed. On the other hand, the RNG k-
su,1 ¼ 2.00, su,2 ¼ 1.17, g2 ¼ 0.44, b2 ¼ 0.083, b* ¼ 0.09. Moreover, ε is as good as the SST k-u in the central and perimeter points,
blending functions are used to achieve a smooth transition be- where the flow is more turbulent and air speed is higher, while it
tween the two models, namely k-ε and keu.
 F. Babich et al. / Building and Environment 123 (2017) 37e49 47

Fig. 14. Measurements and CFD results comparison e perimeter points (SST k-u tur- Fig. 16. Measurements and CFD results comparison e perimeter points (k-u turbu-
bulence model). lence model).

Fig. 15. Measurements and CFD results comparison e perimeter points (RNG k-ε tur- Fig. 17. Measurements and CFD results comparison e perimeter points (k-ε turbulence
bulence model). model).

shows its main limitations by failing to predict the air speed farther
combining the two things, the simulation stopped due to the same
away from the fan.
overflow error after some iterations. Also a finer mesh was used,
In view of the fact that in none of these four turbulence models
but unsuccessfully. Both using a finer mesh and reducing the time
the agreement between measured and simulated values exceeded
step would also have the disadvantage of exponentially increasing
50% in the 12 perimeter points, two more advanced turbulence
the time required to reach a converged solution. Thus, evidence
models have been considered, namely the BSL [44] and SSG [43]
suggests that more advanced turbulence models such as Reynolds
Reynolds stress models. Theoretically, these models are more
stress models and also Large Eddy Simulations (LES) require a
suitable for complex flows, such as the swirling flow generated by
complete explicit model of the fan, the use of a moving mesh, and
a ceiling fan. Practically, their application proved to be fairly
therefore higher computational power. However, this would
difficult. Without changing any other input data or element of the
significantly limit the applicability of the model and increase the
model, the simulation immediately stopped due to fatal overflow
potential sources of uncertainties.
in linear solver using both models. Two relatively easy-to-test
changes were then tried: reducing the time step by setting a
limit for the maximum, rather than the RMS, Courant number
4.2. Required computational power
equal to 1, and using the best available results, namely those
calculated with the SST k-u model, as initial conditions. Even
In this study, simulations have been run on the High
48 F. Babich et al. / Building and Environment 123 (2017) 37e49

Performance Computer System at Loughborough University, a 5. Conclusions

2460-core 64-bit Intel Xeon cluster. Two nodes and 20 cores per
node, therefore 40 cores in total, have been used in these simula- This paper presents research findings on a transient three-
tions in order to make efficient use of the cluster architecture, and dimensional CFD implicit model of a ceiling fan. The research
reduce the computational time. In the pre-processing, the meshes question was whether a simple implicit model, that combines ac-
have been generated using a work station equipped with an Intel curacy with efficient computation, can be used for accurate thermal
Xeon E5520 CPU and 24 GB of RAM. All the other pre- and post- comfort studies. In order to validate the simulation results and
processing activities have been completed using a laptop with an therefore to fully address the research question, a comparison with
i5-3320M CPU and 8 GB of RAM. experimental data has been presented.
Due to extra terms in the two additional k and ε, or u, equations, The results confirm that the model developed is able to accu-
the CPU time required to achieve a converged solution using the rately, qualitatively and quantitatively, predict the air flow gener-
SST k-u and RNG k-ε models is twice the time required when the k- ated by a typical Indian ceiling fan in a room. When using the SST k-
u model is chosen, and three times higher than the time required u turbulence model, 83% of the simulated values are within the
when the standard k-ε model is selected (see Table 7). In detail, the error bars of the respective measured value, and both the swirling
SST k-u and RNG k-ε are both able to predict the meandering and downward air movements are effectively modelled.
behaviour of the plume over time (see Fig. 9), which means that Due to the small number of input parameters required, the
there are relevant differences between each time step and the relatively low CPU time required, and the use of an unstructured
following one, and therefore more iterations per time step are mesh, this model can be effectively used in any study in which it is
required to reach convergence. On the other hand, when the other important to accurately model the air movement generated by the
two simpler models (standard k-u and standard k- ε) are used, the fan, such as thermal comfort or air quality and contaminant dis-
differences between each time step and the following one become tribution research. The choice of the most appropriate turbulence
negligible after the initial time steps, and therefore only a few it- model is essential, and the best results have been achieved using
erations per time step are required thereafter. the SST k-u turbulence model, which is available in most of the CFD
Although there is no direct linear relationship between the packages.
number of cores and the CPU time due to the time required to split, Due to the use of a simplified implicit model and RANS turbu-
and then recombine, the results which increases with the number lence models, the major differences between measurements and
of cores, these total CPU times (see Table 7) could be reduced if simulated values occur in the narrow region below the perimeter of
more cores were available. Considering the CPU time and the ac- the fan, where there are rapid temporal and spatial variations in the
curacy of the results, the SST k-u model is still the best choice. flow field. Further work is required to investigate the possibility of
using more advanced Reynolds stress turbulence models and LES,
4.3. Usability of the developed model balancing the accuracy of the results and the computational effort
required to achieve a reliable converged solution.
Compared with the models previously developed [21e24], the More advanced measurement techniques such as particle image
model presented here is applicable to a wider range of fluid flow velocimetry could be used in the future to provide a better un-
problems due to the small number of input parameters required, derstanding of the flow field generated by the ceiling fan and better
the reasonably low CPU time required, and the use of an unstruc- validation data, including data for lower fan rotational speeds. This
tured mesh, which can accommodate complex geometries. This would also ease the use of more advanced turbulence models.
implicit model can therefore be used whenever a given piece of
research requires a focus on the flow field generated by a ceiling Acknowledgements
fan, not on the actual design of the fan.
When assessing the effect of air movement on human thermal This research was financially supported by the Engineering and
comfort using advanced thermal comfort models such as the IESD- Physical Sciences Research Council (EPSRC) via the London-
Fiala model [9,10], this CFD model is able to effectively predict the Loughborough Centre for Doctoral Research in Energy Demand
air flow generated by the fan at any distance from the fan and (LoLo) (grant number EP/H009612/1), and by the British Council
therefore it reliably estimates the air flow on a person. The limited under the Global Innovation Initiative, the latter involving an in-
agreement between measurements and simulated values at the ternational research collaboration between UC Berkeley (USA),
perimeter of the fan has little implication on the applicability of this CEPT University (India), Loughborough University and De Montfort
model in thermal comfort studies, since this perimeter region is University (UK). The authors express their gratitude for this
less than 10 cm wide, while any human body, regardless of age, support.
gender, height, and weight, is significantly bigger. Moreover, it is
almost always possible for a person to then move a short distance to References
increase or decrease the air movement on their body. Similarly, if
the spread of pollutants in a room is to be investigated, then any [1] S. Manu, Y. Shukla, R. Rawal, L.E. Thomas, R. de Dear, Field studies of thermal
comfort across multiple climate zones for the subcontinent: India Model for
error in the perimeter region does not affect the overall distribution Adaptive Comfort (IMAC), Build. Environ. 98 (2016) 55e70, http://dx.doi.org/
of the contaminants in the space, but only their concentration and 10.1016/j.buildenv.2015.12.019.
movement in that narrow perimeter region. [2] M. Indraganti, Thermal comfort in apartments in India: adaptive use of
environmental controls and hindrances, Renew. Energy 36 (2011)
1182e1189, http://dx.doi.org/10.1016/j.renene.2010.10.002.
[3] USAID PACE-D, HVAC Market Assessment and Transformation Approach for
Table 7 India, 2014. Washington DC.
Duration of the simulations and relative agreement with measurements. [4] N. Shah, W. Paul, P. Amol, Cooling the Planet: Opportunities for Deployment of
Superefficient Room Air Conditioners, 2013.
Simulation Total CPU time Agreement with measurements
[5] DSLDC, Hourly Demand Data from the State Load Dispatch Center, 2012.
SST k-u 6 h, 35 m 83% [6] P.O. Fanger, Thermal Comfort, McGraw-Hill, New York, USA, 1970.
RNG k-ε 6 h, 00 m 67% [7] G. Brager, H. Zhang, E. Arens, Evolving opportunities for providing thermal
Standard k-u 3 h, 34 m 50% comfort, Build. Res. Inf. 43 (2015) 274e287, http://dx.doi.org/10.1080/
09613218.2015.993536.
Standard k-ε 2 h, 10 m 47%
[8] ASHRAE. Standard 55-2013, Thermal Environmental Conditions for Human
 F. Babich et al. / Building and Environment 123 (2017) 37e49 49

Occupancy, 2013. [26] ANSYS. CFX 2016. http://www.ansys.com/products/Fluids/ANSYS-CFX

[9] D. Fiala, K. Lomas, M. Stohrer, A computer model of human thermoregulation (Accessed 23 November 2016).
for a wide range of environmental conditions: the passive system, J. Appl. [27] Q. Chen, Ventilation performance prediction for buildings: a method overview
Physiology 87 (1999) 1957e1972. and recent applications, Build. Environ. 44 (2009) 848e858, http://dx.doi.org/
[10] D. Fiala, K. Lomas, M. Stohrer, Computer prediction of human thermoregula- 10.1016/j.buildenv.2008.05.025.
tory and temperature responses to a wide range of environmental conditions, [28] BEE, Top Ten Ceiling Fans in India, 2016. https://www.bijlibachao.com/top-
Int. J. Biometeorology 45 (2001) 143e159. ten-appliances/best-ceiling-fan-khaitan-havells-usha-crompton-greaves-
[11] P. Cropper, T. Yang, M. Cook, D. Fiala, R. Yousaf, Coupling a model of human orient-in-india-electricity-consumption.html (Accessed 16 February 2017).
thermoregulation with computational fluid dynamics for predicting human- [29] Loveday D, Webb L, Verma P, Cook M, Rawal R, Valodaria, et al. The role of air
eenvironment interaction, J. Build. Perform. Simul. 3 (2010) 233e243, http:// motion for providing thermal comfort in residential/mixed mode buildings: a
dx.doi.org/10.1080/19401491003615669. multi-partner global innovation initiative (GII) project. 9th Windsor Confer-
[12] G. Brager, R. de Dear, A standard for natural ventilation, ASHRAE J. 42 (10) ence Making Comfort Relevant Network for Comfort and Energy Use in
(2000) 21. Buildings, Windsor, UK, 7th-10th April, 2016, p. 947e962.
[13] F. Nicol, M. Humphreys, Derivation of the adaptive equations for thermal [30] ANSYS. ICEM CFD 2016. http://resource.ansys.com/Products/OtherþProducts/
comfort in free-running buildings in European standard EN15251, Build. En- ANSYSþICEMþCFD. (Accessed 23 November 2016).
viron. 45 (2010) 11e17, http://dx.doi.org/10.1016/j.buildenv.2008.12.013. [31] A. Jain, R.R. Upadhyay, S. Chandra, M. Saini, S. Kale, Experimental investigation
[14] F. Babich, M. Cook, D. Loveday, P. Cropper, Numerical modelling of thermal of the flow field of a ceiling fan, in: ASME 2004 Heat Transfer/Fluids Engi-
comfort in non-uniform environments using real-time coupled simulation neering Summer Conference, American Society of Mechanical Engineers,
models, in: Proceedings of Building Simulation and Optimisation 2016: 3rd 2004, pp. 93e99.
IBPSA-England Conference, IBPSA, Newcastle upon Tyne, UK, 2016, pp. 4e11. [32] F.R. Menter, Performance of popular turbulence models for attached and
[15] C. Huizenga, Z. Hui, E. Arens, A model of human physiology and comfort for separated adverse pressure gradient flow, AIAA J. 30 (1992) 2066e2072.
assessing complex thermal environments, Build. Environ. 36 (2001) 691e699, [33] F.R. Menter, Improved Two-equation Keomega Turbulence Models for Aero-
http://dx.doi.org/10.1016/S0360-1323(00)00061-5. dynamic Flows - NASA Technical Memorandum TM-103975, 1992. Ames, CA.
[16] J. Stolwijk, A Mathematical Model of Physiological Temperature Regulation in [34] F.R. Menter, Two-equation eddy-viscosity turbulence model for engineering
Man, 1971. Washington DC. applications, AIAA J. 32 (1994) 1598e1605.
[17] J.A. Stolwijk, Mathematical Model of Thermoregulation. Physiological and [35] F.R. Menter, Eddy-viscosity transport equations and their relation to the
Behavioral Temperature Regulation, 1970, pp. 703e721. keepsilon model, J. Fluids Eng. -Transactions ASME 119 (1997) 876e884.
[18] N. Gao, J. Niu, H. Zhang, Coupling CFD and human body thermoregulation [36] F.R. Menter, M. Kuntz, R. Langtry, Ten Years of Industrial Experience with the
model for the assessment of personalized ventilation, HVAC&R Res. 12 (2006) SST Turbulence Model, in: Proceedings of the Fourth International Symposium
497e518, http://dx.doi.org/10.1080/10789669.2006.10391191. on Turbulence, Heat and Mass Transfer, Begell House, Redding, CT, 2003.
[19] N. Gao, H. Zhang, J. Niu, Investigating indoor air quality and thermal comfort [37] V.C. Patel, W. Rodi, G. Scheuerer, Turbulence models for near-wall and low
using a numerical thermal manikin, Indoor Built Environ. 16 (2007) 7e17, Reynolds number flows-a review, AIAA J. 23 (1985) 1308e1319.
http://dx.doi.org/10.1177/1420326X06074667. [38] B. Launder, D. Spalding, The numerical computation of turbulent flows,
[20] L. Schellen, M. Loomans, B.R.M. Kingma, M. De Wit, A.J.H. Frijns, W.D. van Comput. Methods Appl. Mech. Eng. 3 (1974) 269e289.
Marken Lichtenbelt, The use of a thermophysiological model in the built [39] D.C. Wilcox, Reassessment of the scale-determining equation for advanced
environment to predict thermal sensation: coupling with the indoor envi- turbulence models, AIAA J. 26 (1988) 1299e1310.
ronment and thermal sensation, Build. Environ. 59 (2013) 10e22, http:// [40] D.C. Wilcox, Comparison of two-equation turbulence models for boundary
dx.doi.org/10.1016/j.buildenv.2012.07.010. layers with pressure gradients, AIAA J. 31 (1993) 1414e1421.
[21] R. Bassiouny, N.S. Korah, Studying the features of air flow induced by a room [41] D.C. Wilcox, Turbulence Modelling for CFD, 1993. La Canada, CA.
ceiling-fan, Energy Build. 43 (2011) 1913e1918, http://dx.doi.org/10.1016/ [42] D.C. Wilcox, Simulating transition with a two-equation turbulence model,
j.enbuild.2011.03.034. AIAA J. 32 (1994) 247e255.
[22] Momoi Y, Sagara K, Yamanaka T, Kotani H. Modeling of Ceiling Fan Based on [43] S. Wallin, A. Johansson, An explicit algebraic Reynolds stress model for
Velocity Measurement for CFD Simulation of Airflow in Large Room. Pro- incompressible and compressible turbulent flows, J. Fluid Mech. 403 (2000)
ceedings of 9th international conference on air distribution in rooms, vol. 1, 89e132.
Coimbra, Portugal: 2004. [44] F.R. Menter, Zonal two equation k-turbulence models for aerodynamic flows,
[23] Momoi Y, Sagara K, Yamanaka T, Kotani H. Modeling of prescribed velocity AIAA J. (1993) 2906.
generated by ceiling fan based on velocity measurement for CFD simulation. [45] ANSYS. CFX, Solver Modeling Guide, 2015.
Proceedings of 10th international conference on air distribution in rooms, vol. [46] I. Celik, Procedure for estimation and reporting of uncertainty due to dis-
1, Helsinki, Finland: 2007. cretization in CFD applications, J. Fluids Eng. -Transactions ASME 130 (2008)
[24] E. Adeeb, A. Maqsood, A. Mushtaq, Effect of number of blades on performance na. doi:0.1115/1.2960953.
of ceiling fans, MATEC Web Conf. 28 (2015) 2002, http://dx.doi.org/10.1051/ [47] T. Chai, R.R. Draxler, Root mean square error (RMSE) or mean absolute error
matecconf/20152802002. (MAE)? -Arguments against avoiding RMSE in the literature, Geosci. Model
[25] S. Zhu, J. Srebric, S.N. Rudnick, R.L. Vincent, E.A. Nardell, Numerical modeling Dev. 7 (2014) 1247e1250, http://dx.doi.org/10.5194/gmd-7-1247-2014.
of indoor environment with a ceiling fan and an upper-room ultraviolet [48] H.K. Versteeg, W. Malalasekera, An Introduction to Computational Fluid Dy-
germicidal irradiation system, Build. Environ. 72 (2014) 116e124, http:// namics: the Finite Volume Method, Pearson Education, 2007.
dx.doi.org/10.1016/j.buildenv.2013.10.019.