---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
---------------------------------------------------------------------------------------------------------------------------------

Multiple Model Software for Airflow and Thermal
Building Simulation. A case study under tropical
humid climate, in Réunion Island.
H. BOYER
J. BRAU
J.C. GATINA

Université de la Réunion(*) / INSA de Lyon (**)
CETHIL, INSA de Lyon (**)
Université de la Réunion (*)

Abstract : The first purpose of our work has been to allow -as far as heat transfer modes,
airflow calculation and meteorological data reconstitution are concerned- the integration of
diverse interchangeable physical models in a single software tool for professional use,
CODYRUN. The designer's objectives, precision requested and calculation time consideration,
lead us to design a structure accepting selective use of models, taking into account multizone
description and airflow patterns. With a building case study in Reunion Island, we first analyse
the sensibility of the thermal model to diffuse radiation reconstitution on tilted surfaces. Then,
a realistic balance between precision required and calculation time leads us to select detailed
models for the zone of main interest, but to choose simplified models for the other zones.

I ) Display of the simulation tool :
Born from a joint research project
involving both the Université de la Réunion and
INSA de Lyon, this work aims at producing an
efficient building thermal simulation tool,
including some research and conception aspects,
and taking into consideration different types of
climates. More precisely, it is a multizone software
integrating both natural ventilation and moisture
transfers, called CODYRUN.

a) Architecture and original aspects
The three main parts are the building
description, the simulation program and the
operation of the results. As far as the description is
concerned, we have been brought up to break down
the building into three types of entities which are
the following ones : firstly, the Zones (from a

__________________________________
(*) Laboratoire de Génie Industriel, Université de
la Réunion - Faculté des Sciences, 15 rue Cassin
97489 Saint-Denis Cedex. Ile de la RéunionFRANCE
(**) Institut National des Sciences Appliquées.
Equipe Equipement de l'Habitat, CETHIL,
Bat. 307, 20 Avenue Albert Einstein.
69621 Villeurbanne Cedex, FRANCE

thermal point of view), the Inter-zones partitions
(zones separations, outdoor being considered as a
particular zone) and finally, the Components (i.e.
the walls, the glass partitions, air conditioning
system sets, so on and so forth). For a simulation,
the Project notion includes a weather input file, a
building and then, a result file. The following
treelike structure illustrates our organisation :
Meteorological file
Project

Zones

Components

Inter-Zones

Components

Building
Result file

Fig. 1 : Data organization
During a simulation, one of the most
interesting aspects is to offer the expert thermician
a wide range of choices between different heat
transfer models and meteorological reconstitution
parameter models. The aim of this simulation may
be, for a given climate, to realize studies of the
software sensitivity to these different models, in
order to choose those that should be integrated in a
suitable conception tool.
In the second part of this article, the first
application concerns the comparison of two models
of sky diffuse radiation. In the same way, for a

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------given climate, depending on the objective sought
The thermal model relies on INSA's
during the simulation (through tracking the
previous simulation code, CODYBA (BRAU,
temperatures or estimation of a yearly energy
1987). With the usual physical assumptions, we
consumption), it is interesting to have the liberty of
use the technique of nodal discretisation of the
selecting the models to be involved.
space variable by finite difference. In addition, the
In most existing simulation software, the
mass of air inside one zone is represented by a
choice of the models being already done, their
single thermal capacity. Thus, for a given zone, the
application is global to the whole building. If these
principle of energy conservation applied to each
models are complex, the multi-zone feature as well
concerned wall node, associated with the sensible
as airflow patterns leads quickly to calculation time
balance of the air volume, constitute a set of
which is not compatible with a conception tool. It
equations, that can be condensed in a matricial
has seemed interesting to us to allow, for some of
form :
the phenomena, a selective use of the models.
Thus, the choice of the indoor convection model is
dT
realized during the definition of a zone or else, the
(1)
C
 A T +B
dt
choice of the conduction model is done during the
description of a wall. The aim is then to link the
The automation of the setting up of the
level of complexity of the models concerning one
previous equation requires a decomposition of A
entity to the interest beared for this same entity.
into twelve elementary matrices, each one having a
And that is how we will show in the second part of
particular feature.
this article the importance in choosing a detailed
model for the indoor convection in the main
dT
interest zone and simplified models for the other
C

zones.
dt
Acond + A cvi _lin+ ...+ A connex T 
Bint_load +...+B connex

(2)

Airflow Model

Step

Thermal Model

Humidity Model

Fig. 2 : General flowchart
In relation with calculation program, the
software executes at each step the reckowing of
airflow patterns, temperature field and specific
humidity of each zone.
The most simplified airflow model
considers as known the airflow rates between all
zones, whereas the more detailed model calculates
with a pressure model the airflow through each of
the openings, which can be large ones. The
building is represented as a network of pressure
nodes, connected by non-linear equations giving
the flows as a function of the pressure difference.
This detailed airflow calculation goes through the
iterative solution of the system of non linear
equations made up with air mass conservation
inside each zone.

Thus, Acond gathers features linked to heat
conduction through the walls whereas Acvi_lin
gathers the features depending on indoor and
linearised convective exchanges. In the same way,
for each step, the filling up of B vector is made
easier by its decomposition into fifteen elementary
vectors. The physical coupling of the considered
zone to the other ones is realized through an
iterative connection process, through the filling up
of a matrice Aconnex and a vector Bconnex.
At each step, the resolution of equation (1)
uses an implicit finite difference procedure and the
coupling iterations between the different zones
make it possible to calculate the evolution of
temperatures as well as those of sensible powers
needed in case of air conditioning. Having in mind
a compromise precision/calculation time, it is to be
noticed that the thermician using the software
keeps the control of the solving methods, of the
iterations (mainly concerning connection process
and airflow model) and of the different
convergence criterions.

b) WINDOWS Front-End :
Developed on PC micro computer with
Microsoft WINDOWS interface, implemented in
C-language, the software beneficits of all the user-

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------friendliness required for
a conception tool
The chosen building is a cubic shape with a
(windowing, mouse, ...).
side of 6 m, in contact with the ground. We
consider a thermal split up in three zones : the
A more technical aspect linked to this
ground floor, the eastern first floor and the western
system is very interesting for simulation tools
first floor. The following sketch illustrates our
based on PC, such as our software : memory
description :
managing supplied by WINDOWS allows to
allocate much more than the classical 640 Ko
limit, which is most necessary as regards the
capacity and the variety of the matrices and vectors
that interfere in the simulation process.
The user's interface proposes to begin the
description of a building with the following
window :

East
South
Fig 4 : Building sketch, Southern face

Fig. 3 : Building description window
The push buttons zone, interzone and
components make possible the access to this
previously defined entities. The access to the
models linked to the entity "building" is also
possible from this window through the push
buttons of the screen part called "models". It is
possible in this way to select the chosen models for
outside heat convection transfer, reconstitution of
meteorological parameters (diffuse and sky
temperature) as well as the airflow model.

The main characteristics are the following
ones : let us suppose all the walls are made up with
of dense concrete of 12 cm, the slab on grade being
of 30 cm in the same material. On the eastern and
western sides are display bay windows (simple
glass) of an elementary surface of four square
meters and of usual optical and conductive quality.
Later, the conductive model kept for each wall and
glass is a simple two capacitors model (R2C).
Besides, we suppose the presence of an air
conditioning system in the western zone, upstairs,
but we shall come back on this component in a
following paragraph. The building is described
with three zones, sixteen inter-zones and twenty
two components.

II ) Case study :
a) Description :
Réunion Island is located in the Indian
Ocean by 55°1 longitude East and 21°5 latitude
South. The climatic conditions being those of a
tropical and humid climate, we have reconstituted
an hourly file gathering all the meteorological
parameters necessary for our simulations, i.e.
mainly solar radiation data, outdoor dry
temperature and moisture rate, and wind. With
CODYRUN, for some periods of the considered
year, we propose to reproduce a building indoor
thermal conditions.

b) Diffuse radiation model choice :
Under tropical-humid climate, the readings
made on the short wave diffuse radiation show the
importance of this kind of inputs. A quick analysis
of the graphs on the site of Réunion Island shows
an important diffuse radiation when the direct
beam is low. Tropical humid climate leads to take
into account the total diffuse radiation (sky diffuse
and ground reflection), when designing the
protection devices, i.e. screens, shadow masks, ...
(Cabirol, 1984).
Most of solar energy calculations consider
this diffuse radiation as being isotropic. As given

d ( s,  ) =

(4)

1  cos s
Max(cos i, 0)
( F . C ( s ).
 ( 1  F ).
) . dh
2
sin h
2

6

with C( s )  1.00115  3.54.10 s  2.46.10 s
Gh
dh
F  1
(1 
)
GhExt
Gh
and

i

mean angle of incidence of beam radiation
with respect to the surface normal
h
solar altitude
 Gh
site horizontal global radiation
 GhExt extraterrestrial global radiation
Considering the sharing (between direct
beam and diffuse part) of the solar heat gains
entering a zone (through the windows), the model
considers that the direct beam is incident on the
floor and that the diffuse part is split up depending
on the proportion of the surfaces. Through the
optical characteristics of the materials in which the
surfaces of the walls in consideration are made up,
it is necessary to solve a linear set of equations to
obtain the radiation absorbed by the surface nodes
of the walls. These radiations are constituting one
of the elementary vector, Bswi.
For the study, we have taken from the yearly
file a sequence of two extreme consecutive days,
one cloudy and the other sunny. The following
graphs show the evolution of the horizontal direct
beam and diffuse radiation for the simulation
period.
1200

Horizontal radiation (W/m²)

Dh
1000

dh

800
600
400
200
0
1

6

11

16

21

26

Time step

31

36

41

46

Indoor diffuse radiation / 100
(W)

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------dh the diffuse and horizontal radiation, on a tilted
Fig. 5 : Horizontal radiation
plane (azimuth , inclination s), this diffuse
Through the glasses, one part of the energy
radiation is :
load is due to the beam part of the transmitted
radiation and the other part to its diffuse part. For
1  cos s
(3)
d ( s,  )  (
). dh (W/m²)
the first day, the following graphic plots show the
2
evolution of the total diffuse incident radiation in
Meanwhile, anisotropic models having been
the eastern floor zone, cases isotropic and
validated (Gueymard, 1987), we have integrated
anisotropic (Willmott).
the Willmott model (Lebru 1983), besides the
isotropic model for reconstitution of incident
diffuse radiation on a tilted plane. With the
9
following notations, the proposed expression is :
8
Willmott

7

Isotropic

6
5
4
3
2
1
0
1

6

11

16

21

Time step (Hour)

Fig. 6: Eastern floor diffuse radiation
To integrate a non isotropic diffuse radition
pattern is the same as to consider as non isotropic a
part of the diffuse radiation preliminary calculated
with the isotropic model. As a consequence, the
diffuse radiation incident in a zone will be lesser
with the anisotropic model, as shown by the
previous graph. For the incident direct beam, the
live would show the opposite order.
The following figure allows comparison
between graphs of dry temperatures evolution
concerning the eastern first floor, and that for the
two cases previously mentioned.
Willmott

30

Isotropic

29
28
27
26
25
24
23
1

6

11

16

21

26

31

36

41

46

Time step (Hour)

Fig. 7 : Eastern floor dry temperature
Considering our building and the mentioned
period, the graphs are practically merge. This
simulation shows in our particular case that the
whole building model is little sensitive to the sky
diffuse radiation model. A simple interpretation of
this phenomenon can be delivered. In the
anisotropic case, the diffuse radiation deriving
from the isotropic model is divided into two parts,
one of them is isotropic, the other one directing. As
far as the energetic behaviour of a building is

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------concerned, the coupling of the external short wave
The first one being a linear model with a
radiation and the zone in consideration is realized
coefficient of constant convective exchanges (but
through the glasses. For simple type of glass such
that can be modified by the conceptor), the second
as the one in which our windows are made, the
one, linear with a coefficient depending on the type
bibliographic expressions of the transmittance for
of wall (i.e. floor, ceiling, vertical wall) and the
direct or diffuse radiation are but little different.
last model is a non linear one. In this case, the
Moreover, the diffuse transmittance is not
exchange coefficient depends in a non linear way
dependent of the incidence angle. Thus, the
on the temperature difference between air and the
difference in dividing the energy load according to
given surface. The integration of such a model in
the isotropic criterion will have but little incidence
our system goes then through the filling up of a
on the zone's thermal behaviour.
vector dealing with non linear and convective
flows, Bcvi_lin and through a process of iterative
resolution with a convergence criterion that can be
c ) Choice of the indoor convection
modulated for example as regards the indoor dry
model
temperature. This non-linear model needs several
resolutions of the equation (1), at each step. With a
c -1) HVAC System Component :
convergence criterion of 10-3 °C on the air
temperature of the zone in consideration, the
The western zone being air conditioned, for
number of iterations is three or four. The use of
a better understanding of the following graphs, we
this model tends then to penalize the tool as far as
first display and comment the data window capture
the calculation time is concerned, but this model
of the Air Conditioning System component.
can be considered as the most reliable.

c -3) Exploitation
Considering the conditions of simulation
of the previous a) paragraph, our concern is now
the evolution of the dry temperature inside the
western zone floor as well as the sensible power
needed. Firstly, we consider the A case in which
the convection integrated with the help of a linear
model with a constant coefficient (5 W/m².K), in
each zone. The following graphs are obtained :
23

22

21

Fig 8 : Component description window
Thus, for an air conditioning system,
information relative to the sensible and latent loads
are to be entered. For our study, we'll consider the
sensible part. We must then define the threshold
temperature values (high and low), the hourly
schedules and the available heating and cooling
powers. Most of the time, simulation softwares
consider the heating power as convective.
However, to allow a better integration of systems
such as heating floors, it is possible to choose for
the involved power convective and radiative ratios.
In particular case of the air conditioning system
sizing, the software considers that the available
power is infinite.

c -2) The indoor convection
The choice of the indoor convection model
is made during the zone description. We have set
up three models as regard the indoor convection.

20

19
1

12

23

Fig. 9 : Air Temperature, Western floor
1

6

11

16

21

0
-0,5
-1
-1,5
-2
-2,5
-3

Fig. 10: AC Power, Western floor

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------The graph of the power in dotted
with graphs very close to the case B, both in
corresponds to equipment sizing, in which case we
temperature and in power. In parallel, the
consider the air conditioning power as infinite.
comparison of the simulation times (outside
This graph enables us to determine the cooling
building initialization period) for a given period of
power to be set up, about 3kW, in order to respect
the day is given in the following board.
at any time the specified 20°C.
If the cooling equipment power is inferior,
(2 kW in our case), the result is an increase of the
indoor temperature during the second half of the
day (Fig. 9). Quantitatively, the two previous
graphs allow a designer to realize the overheating
(and its length) due to the under-sizing of the air
conditioning system.
In case B, we integrate the indoor
convection as non linear, for the three zones. This
simulation is very costly as far as calculation time
is concerned. In the last case, C, we'll deal with the
non linear convection only in the zone in which we
are most interested, i.e. the western floor, and we'll
use the linear model with a constant coefficient for
the two other zones.
The superposition of the graphs relative to
the three cases produces the following results.

Indoor Temperature (°C)

23
case A

22

Case
B
A
C

Temp. Err. (°C)
0 (reference)
1.2
0.1

Power Err. (kW)
0 (reference)
0.5
0.1

Time (mn"s)
2"53
0"54
1"35

Thus, there is approximativaly a one to a
half ratio between the simulation time in case B
and C. This ratio can be found again by other
simple considerations. Indeed, it is possible in first
approximation to suppose as constant the time
necessary to set up and solve the state equation of
one zone, i.e. t seconds. If Nb_zones is the number
of zones, in the case A, the resolution for one step
requires tA = (Nb_zones . t) seconds. Let us
suppose constant and equal to i the number of
iterations in a zone, interations introduced by the
non linear model of the indoor convection. In case
B, all the different zones being concerned with the
non linearity, the calculation time equals tB =
(Nb_zones . i . t) seconds. In case C, this time
becomes :
tC = (i + (Nb_zones - 1) . t ) seconds.

case B
case C

With values corresponding to our case
(Nb_zones = 3, i = 3), the ratio (t C/tB) is quite
close to 0.5.

21

20

19
1

6

11

16

21

Fig. 11: Air Temperature, Western floor
1

6

11

16

21

AC Sensible power (kW)

0
case A

-0,5

case B

-1

case C

-1,5
-2
-2,5

Fig. 12: AC Power, Western floor
If we consider as right the values proposed
for case C, the inaccuracy of the model with a
constant coefficient applied to all the zones is
clearly visible. The mistake on the temperature
reaches more than one degree and these on the
power needed is of 0,5 kW. On the contrary, the
case C shows the accuracy of the selective
application of the non linear convection model,

III) Conclusion
Through the previous simulations, we have
introduced CODYRUN, focusing on a few points.
Many aspects developed to this day in this software
have been purposely put aside, in particular in
regard to airflow simulation, which is a dominant
heat transfer mode under tropical climate. As far as
heat conduction transfer is concerned, the different
implemented models haven't been reviewed. In the
same way, as regard this transfer mode, tools such
as physical aggregation of walls or modal
reduction (Roux 1988), allowing a notable
decreasing of calculation time, will be the objects
of future improvements. So, we join one of the
preoccupations of his study, i.e. the balance
between precision required and calculation time.

REFERENCES
Auffret, P.; Gatina, J.C.; Hervé ,P. 1984. "Habitat
et climat à la Réunion. Construire en pays tropical

---------------------------------------------------------------------------------------------------------------------------------IBPSA’93, International Building Performance and Simulation Association, Adelaïde, Australie
--------------------------------------------------------------------------------------------------------------------------------humide." Documents et recherches n°11.
Université de la Réunion.
Brau, J.; Roux, J.J.; Depecker., P.; Croizer, J.;
Gaignou, A.; Finas, R. 1987."Micro-informatique
et comportement thermique des bâtiments en
régime dynamique." Génie Climatique, 1987,
n°11 (Oct-Nov), 15-23
Cabirol, T. ; 1984; "Habitat bioclimatique :
l'incidence de l'ensoleillement et du vent". Afrique
Expansion, 1984, n° 6 (Nov. ), 45-48
Lebru, A. 1983. "Estimation des irradiations
solaires horaires dans un plan quelconque à partir
de la donnée de l'irradiation horaire globale (et
éventuellement diffuse) horizontale." Research
Report N° 239 cahier 1847. CSTB, Sophia
Antipolis FRANCE.
Roux, J.J.; Depecker; P., Krauss, G.; 1988.
"Pertinence and performance of a thermal model
adapted to CAD context." In Proceedings of the
sixth International PLEA Conference (Porto,
Portugal, July 27-31),749-754
Geymard, C.; 1987; "An anisotropic solar
irradiance model for tilted surface and its
comparison with selected engineering algorithms."
Solar Energy, vol. 38, N° 5, 367- 386
Walton, G.N.; 1984; National Bureau of Standards,
Washington, DC. "A computer algorithm for
predicting Infiltration and Interroom Airflows".
ASHRAE Transactions 84-11, N°3

