- Email: [email protected]
Solar heating system design: The part of simulation models
Energy and Buildings, 3 ( 1 9 8 1 ) 197 - 212 © Elsevier S e q u o i a S.A., L a u s a n n e - - P r i n t e d in T h e N e t h e r l a n d s
197
Solar Heating System Design: The Part of Simulation Models P. C H O U A R D
E.D~F. Les Renardidres, F-77250 Moret/Loing (France)
INTRODUCTION
Models, why ? All the branches of physics are deeply concerned with simulation: to simulate by speculative ways the behaviour of physical systems is a very powerful tool at the disposal of research workers and engineers in all the fields of their activities. The thermal evaluation of buildings is specially concerned with this tool because all the others methods in competition with simulation such as full-size field-experiments or laboratory tests are always very expensive and restricted: expensive because of high cost of the p r o t o t y p e s themselves and of the instrumentation needed for proper checking on performances, and restricted because it is impossible to scan all the parameters involved and to get a sufficient knowledge of them in their full range o f variation. Today, a lot of applications of solar energy systems for space heating, cooling, and domestic h o t water are available; this is fully consistent with these thermal simulation methods, because solar systems being basically driven by transient phenomena often require accurate methods for the prediction of their thermal performance or life-cycle cost analysis. As a matter of fact, the thermal models will be able to: a) get a better appreciation of the physical p h e n o m e n a involved and to evaluate their o u t c o m e from a practical point of view; b) predict the behaviour of materials, i.e. to estimate their reliability and duration; c) size the systems and their c o m p o n e n t s as a function of constraints and objectives. Moreover the research engineer also needs such a model when, having to conceive new techniques or new applications, he has to
know the sensibility of a given parameter from a technical and/or economical point of view.
Various categories of models In accordance with the numerous features listed above, a great number of solar models n o w exist; b u t the objectives, the accuracy, the costs, and the procedure of each model are very different one from another, and a classification of them seems necessary. A good criterion might be their relative complexity according to the c o m p u t e r size needed, or the expected cost of a typical testproblem. We shall in fact distinguish two categories of solar models: the "simplified" models and the "comprehensive" ones. The former are specially designed for consulting engineers at the stage of preliminary design, who have to get a first idea of a given system configuration quickly and at a low cost. At this stage, it is completely impossible to use big computer or comprehensive programs far t o o expensive to compete with the simplified ones when the components of the system or the configuration of the building are not yet defined or even when the "solar choice" has n o t yet been done by the customer! Several models belonging to this category already exist in Europe or in the U.S.A., but they cannot be compared with the comprehensive ones from an accuracy point of view. These comprehensive models are in fact used by research people because they are adapted to fundamental studies or to advanced projects, due to their high degree of accuracy. On the other hand, the calculations involved are so complex that a more or less sophisticated c o m p u t e r is always nec~.ssary in that case.
198 First, we shall describe in this paper the comprehensive models taking the computer program " C L I M A S O L " designed by Electricit~ de France as an example; in a second step, we shall describe the principles of the simplified simulations and take the U.S. f-chart m e t h o d and the French "Bilan therm i q u e " m e t h o d as examples.
COMPREHENSIVE MODELS
General features All the comprehensive models show c o m m o n features because there are often very few ways to reach the above-mentioned objectives. First of all, the computation procedure is always a step-by-step one, which means that a recurring m e t h o d is used: all the equations are solved at a given m o m e n t and these results are used as inputs for the next calculations at the next moment. As a matter of fact, building materials always exhibit a great inertia, therefore the lag between two calculation steps can be up to several hours, according to the materials and to the degree of accuracy needed. A one-hour period seems to be quite satisfactory in most cases. As a consequence of the objectives, specially of the transient phenomena requirements, it is necessary to go into all the details of the description of the building itself. In order to take all the heat transfers between rooms, or between inside and outside, at every m o m e n t into account one has to describe the location of the building, its environment in relation to other buildings, the room shapes, the physical properties of materials, and so on. A good knowledge of the meteorological environment is also necessary, and the accuracy of its description depends on the accuracy of the calculation procedure itself. As an example, if the calculations are repeated hourly, an hourly description of the climate will be necessary. If such climatic data are not available, an interpolation will be made. Internal heat gains are an essential parameter in m o d e m , well-insulated buildings. Therefore, the comprehensive models attach the greatest importance to them. Lighting (power dissipated, duration, type of lamps), appliances, and people (number, activity on an hourly basis) are terms of the thermal
balance of the building and appear in the equations, as well as the solar gains through the windows and even through the walls. Once all these data are known, the computation itself becomes possible. From a theoretical point of view, and except for the heat transfer through the walls, the algorithms are quite simple. The question is to take all the gains, all the transfers and all the above parameters into account and to write that all the instantaneous energy flows (positive and negative) are balanced at every moment. The difficulties come from the number of iterations (720 hours within a single m o n t h for instance) and from the solution of the transient equations for the heat transfer through walls and floors. Some methods will be described in the next section. For a large number of iterations a more powerful computer will be needed. A great flexibility does exist owing to the computer. The simulation can be done over long periods, often one year when a typical meteorological year is considered, but it could be still longer. Moreover, as all the variables are calculated (temperatures, energy flows, etc.), it is possible to get a good knowledge of the thermal response of the building in addition to its energy consumption. In fact, it is often a question of printed material because it may be difficult to choose the proper variable to be printed! A lot of additional possibilities still exist. Some models include financial computations, life cycle cost analysis, return on investment, and so on. In fact, everything is possible, and the biggest question looks like a dilemma in connection with the computer: if the computer program is too sophisticated, it cannot be handled by engineers or scientists, it becomes an informatical concept without any practical use. This might be a danger.
Heat transfer models Without giving too many details (it would require a few hundred pages), we are just going to set the problem and to describe its solution by two usual algorithms.
The problem o f the wall For instance, let us consider a wall between two mediums at different air temperatures: Tai (inside) and Tao (outside). The temperature T (x, t) in the wall is a solution for the
199
unidimensional heat transmission equation at every m o m e n t a
~2T
- ~x
-
2
8T ~t
where a: diffusivity of the material. The two limiting conditions of this equation result from the thermal equilibrium of the surfaces of the wall (inside and outside) inside: = hi[Tai -- T(1, t)] + ¢R outside: = ho[Tao -- T(0, t)] + ¢A
because we just need to know the inside temperature of the wall at every moment. There are t w o methods adapted to this problem: the finite difference method and the response factor method. The finite difference m e t h o d This is a computer based method. The above equation is discretized in both space and time, by a substitution of finite differences for derivates. The wall is divided into N elementary layers, and with the above equation, the temperature T (x, t + At) of each junction is calculated at time t + A t from the temperatures of the next junctions at time t and/or t + At. Mathematically, the equation ~2T ~T a -~x 2
where X:
conductivity of the wall (W/m °C) inside and outside convective transfer coefficient (W/m e °C) long wave radiation between the in¢R: side surface and its environment radiation absorbed by the outside CA: surface of the wall. Let us remark that if the wall is heterogeneous (n layers with different materials), the problem is identical. Instead of two, there are n + 1 limit conditions all looking like this one
is written like this
hi, ho :
~ki ~X
i=li
---- ~ki+ 1
} ~X/Xi+l
The initial time condition can correspond to a linear distribution of the temperatures in the wall, from T (0, t) up to T (1, t). The problem is to solve the equation 52T ~T
a
-~x 2
a
T(x + Ax, t ) - - 2 T ( x , t) + T ( x -- Ax, t)
=
(Ax) 2
T(x, t + A t) -- T ( x , t) At
which allows T (x, t + A t) to be evaluated as a function of T (x + Ax, t), T (x, t) and T (x -Ax, t); this is the so-called "explicit" method. The equation a
---- 0
~t
~2T
~T
8x 2
~t
could also be written like this a
T(x +Ax, t+At)--2T(x,t+
At)+T(x--AX,
(~x) 2 = T ( x , t + A t) -- T ( x , t)
- --
At
~t '
Outside
Wall
x
T=o ~.
Inside
i. OT
ax
0
Fig. 1. Heat transfer through a wall.
r.i
and thus T (x, t + A t) is evaluated as a function of T (x -- At, t + At), T ( x , t) and T ( x + Ax, t + At); this is the so-called "implicit" method. The former is quite satisfactory b u t the inequation a~t
1
(~x) 2
2
must be fulfilled for the convergence of the algorithm; therefore, the time interval At must be very small (a few minutes). The latter m e t h o d is always convergent, but it is time-
t+At
200
consuming and requires a bigger area in the computer core.
TCmpe,atu,e
\-/ The response factor method This m e t h o d is based on the thermal response of a wall to an elementary stimulus from a given signal. If we admit that all the equations of the systems are linear, and that the thermal properties of the materials are hot time-changing, we find that the response of the wall to any signal is the sum of the elementary responses we already know. Therefore, there are three steps in the procedure. We must: -- split the excitation function into elementary impulses (for instance into temperature steps or triangular functions); --calculate the response of the system (the wall) to each elementary stimulus; -- add up all these responses. As a matter of fact, every time-dependent function can be split into a series of numbers equal to the values of this function at the moments 0, At, 2 At, 3 At .... Thus, a series of triangular functions replaces the actual thermal function (see Fig. 2). Other possibilities exist, but the work of Mitalas [1] demonstrated t h a t the triangular function is best adapted to such a purpose. :;';:::r::o'0
The,real Triangutar functlons
.~r
Fig. 2. Splitting triangles.
2 ~t
3 .~t
4 .~t
of a c o n t i n u o u s
function
into
Therefore, it is only necessary to know the response of a given wall to a triangular stimulus from the outside air temperature. These responses at each time interval are the so-called response factors of the wall (see Fig. 3). Let us take as an example a variation of the outside temperature which generates a variation in the inner temperature distribution of the wall: T (x, t). The knowledge of this distribution leads to the expression for the energy flow crossing the inside surface of the wall
¢ (1,0 =-x\ x=l
4, /
..............
/
~
~ T,m~ ~t
? zt
3 ~t
4 ~r
5 at
6 ~t
Fig. 3. Definition of response factors.
The distribution itself can be evaluated by the finite difference method. It can also be measured if thermal data is available. From a theoretical point of view, the number of thermal response factors of a wall is infinite (as many as the time intervals n.At). As a matter of fact, the response factors become negligible after 20 hours in most practical cases.
Review o f the main existing computer programs It is impossible to describe all of the comprehensive models that exist all over the world. From a general point of view, these models belong either to the "response f a c t o r " category or to the "finite difference" category, but the i n p u t / o u t p u t procedures are often very different from one another according to the objectives of their designers. In fact, we would like only to demonstrate in this section t h a t a great number of models are operative and that very few comparative studies have been done. Therefore our review looks more or less like a list. U.S. models Among the numerous computer programs available in the United States, we shall just mention the most famous ones. The BLAST (Building Loads Analysis and System Thermodynamics) program is a comprehensive set of subprograms for predicting energy consumption in buildings. It has been designed at the U.S. Army Construction Engineering Research Laboratory (CERL). Each of the three major subprograms: --computes the space load on an hourly basis, using the ASHRAE methods; --simulates the air distribution system and calculates the h o t water, chilled water and electric energy demands; - - simulates the central plant (boilers, chillers, solar systems) and computes energy consumption and life cycle cost. The main
201
feature of BLAST is that a good simulation of the heating or cooling equipment is done on an hourly basis. TRNSYS has been designed at the University of Wisconsin by Klein et al. in 1973". It can simulate the dynamic behaviour of solar transient systems. It consists of components models (for collectors, controls, storage tanks, heat exchangers, furnaces, building loads) and an executive routine: the user selects both the components and their design parameters, and specifies the way in which the components are interconnected. The building is n o t simulated, the user only specifies a global loss coefficient (UA or GV in W/°C); on the contrary, all the details of the solar system are described and simulated. The DOE-1 (formerly Cal-ERDA) is sponsored by the U.S. Department of Energy. It is another program available to the public which simulates the thermal load of a building using the response factors methods. The authors** (Lokmanhekim et al.) simplified the calculation algorithms and as a result the program is very fast and its running cost rather low. European models Several models are available of which we will only mention the most detailed ones. In Belgium, the University of Li6ge has designed several computer programs called LPB1, LPB3, LPB4, each of them exhibiting an increasing degree of complexity in building simulation. In France, the COSTIC (Comit6 Scientifique et Technique des Industries du Chauffage) uses a detailed program which simulates the building by the response factor method and the solar system by modular subroutines. This program has been widely used in order to generate "rules of t h u m b " for domestic water heating. Unfortunately, this program is not publicly available and not much literature has been published.
*S.A. Klein e t al., TRNSYS: A transient system simulation program, Solar Energy Laboratory, University of Wisconsin at Madison, June 1979. **M. Lokmanhekim et al., DOE-2, Building Energy Analysis Group, Energy and Environment Division, Lawrence Berkeley Laboratory, Berkeley, CA, Feb. 1979.
The INSA (Institut National des Sciences Appliqu6es) high school in Lyon has done a very useful synthetic j o b in the field of comprehensive programs. Their model is based on the response factor method. We also mention here the CLIMASOL program, designed at Electricit~ de France, which will be described in the next section. E D F computer program: C L I M A S O L This program basically consists of two sections that can be described separately. The first section (CLIM) computes the instantaneous thermal heating and/or cooling loads of a given building taking the transient phenomena into account. The second section (MASOL) is a detailed model of a solar heating device which simulates the operation of a solar system in order to meet the loads previously calculated. As a result, MASOL gives the energy balance of the system (solar energy used and supplementary energy if any) on a daily, monthly or yearly basis. Thermal loads: CLIM program The general principles of the program are listed hereafter. First of all one describes the building, its climatic location and utilization. The building is divided into " m o d u l e s " which are (by definition) a group of rooms having the same thermal characteristics (inside and outside temperatures, orientation, utilization, etc.). From all this data, the program computes, step by step, the energy requirements in terms of power and power consumption in each module. Data acquisition. The computer reads all the useful data via punched cards of which the main ones are: TITRE (title): this card describes the geographic location of the building and the weather data to be considered. It also defines the category of the building (residential, commercial, etc.). MATI~RIAUX & PAROIS (materials and walls): these cards correspond to the frame of the buildings and give the composition of each wall, window, floor, ceiling, and so on. They describe the physical properties of the building materials (thermal conductivity, density and specific heat) as well as the possible composition of the walls (thickness of each layer, tilt angle from horizontal, surface
202 resistance, shading coefficient of the windowblinds, if any, etc.). R E N O U V E L L E M E N T D ' A I R (air change rate) : this is a card for the fresh air mass flow rate into each module, its humidity and temperature and the electric consumption of fans (if any). It is possible to define several air change ratios as a function of time, according for example to the number of people living in the building. APPORTS (free heat): it is possible to define hourly the rates which heat and moisture are separately given off by people, lighting, and/or appliances in each module. These heat gains play a great rSle in modern, wellinsulated buildings and this is why they are described in detail. PUISSANCES INSTALLI~ES (reference power of the heating or cooling equipment): this card represents the thermal maximum power of the available equipment (fundamental data as far as transient phenomena are concerned). MODULES: detailed description of the rooms (modules) themselves, dimensions, characteristics of the walls (from the MATl~RIAUX cards) and inside design temperatures (minimum and maximum); the contiguous modules are also identified by this card.
Weather data and solar radiation. Data on o u t d o o r weather conditions (dry bulb temperature, relative humidity and solar coefficients) have been stored in the computer m e m o r y for 26 stations in France for a typical year. As far as solar radiation is concerned, the solar position (altitude and azimuth from South) is c o m p u t e d at every m o m e n t taking into account such data as latitude, day of the year, hour of the day, spherical coordinates of the radiation receiver (tilt from horizontal, and deviation from South) and so on. These calculations together with the solar coefficients coming from the meteorological data allow an accurate estimation of the actual incident radiation to be considered for the energy balance of the building. This procedure which is based both on calculation and on meteorological data makes up for the shortage of relevant data (very few weather stations make measurements of radiation on tilted surfaces on an hourly basis). Computation o f thermal phenomena. As in all the models belonging to this category, the
problem is to calculate the air and wall temperature in each module (room) as a function of: the internal heat gains listed above; the design (inside and outside) conditions; and the temperature of the contiguous modules. In order to solve the heat transfer equations in each wall, a finite difference computation method has been chosen; the matrix representing the system is inversed by a ShermanMorrison [2] method. Simplifying assumptions have been made, the justification of which would be quite a lengthy paper itself. These assumptions concern: - - t h e inside temperature which is supposed to be homogeneous in the living space area; - - t h e heat flux which is supposed to be normal to the wall surfaces (heat losses through "thermal bridges" are calculated separately). - - t h e solar radiation transmitted by the window panes which is assumed to be evenly distributed on the floor surface (i.e., no shadow limit is calculated). From all these data, algorithms and assumptions, CLIM can c o m p u t e the thermal heating or cooling loads and the inside temperature in each module twice an hour (i.e., 24 X365 X 2 = 17 520 iteration steps are necessary to simulate a whole year). These loads will be interpreted as input data by the next subroutine: MASOL.
Solar system simulation: M A S O L program Most solar energy systems can be described as a group of three major interconnected subsystems namely collection of solar heat, storage of energy and distribution devices. Basically, MASOL consists of a detailed model of these three sub-systems: each subsystem includes its own inertia, heat losses and heat gains. General description. Figure 4 shows the block diagram of the system. Energy flows are represented by arrows and by the mass flow rates of the fluids. The distribution system includes a recirculated fluid loop (for space heating) and a non-recirculated loop (domestic hot water for example). These loops can be activated together or separately according to the actual system. The control subroutines are also copied closely from the actual control equipment:
203 FR:
d,s,,,b,~,on
collector heat removal efficiency factor UL: collector energy loss coefficient (W/ m 2 °C) On the other hand, the CLIM program computes the solar radiation transmitted by a given surface
Pv = r(i)PD COS i + 0.9 7(o)P d Col~ection
SIo,a~e
D,str,bution
Fig. 4. Schematic diagram of the solar system simulated by MASOL.
-- if the temperature at the outlet of the solar collectors exceeds the lowest temperature of the storage tank, solar heat is transferred from the collectors to the storage tank (mass flow rate at its typical value in the collector-tank loop) ; if the storage temperature exceeds a threshold value and if the rooms need that (positive load as c o m p u t e d by CLIM), stored heat is transferred to the distribution system (flow rate at its typical value in the distribution loop); -- if the load cannot be met by the solar distribution system, the supplementary energy system is activated. -
Thus, we can write (assuming a(i) = a, whatever i is)
Pu = FR (aPv -- ~A T) where Pv comes from CLIM, a is the plate absorptance, and ~ the collector overall energy loss coefficient ( W / m 2 °C); these numbers can best be determined by collector tests in which instantaneous efficiency is measured under specified conditions (e.g. [4] ). The general thermal balance of the collector will be written as dTcAPT
-
H e a t collection. The theory relating instantaneous collector performance to meteorological and operating conditions is well-known from the studies of Hottel et al. [3] (during the years 1942 to 1959). We can thus express the rate at which useful energy is collected as a function of solar radiation, of outside ambient temperature and of fluid temperature. The equation is
M1C1
dt
= ~Pv --/3(TcAPT -- TEXT)
--P 1Cidi(T~
-- Te)
where M i is the equivalent mass capacity of the collector, fluid + material (kg/m 2) C i the heat capacity of the fluid (J/kg °C) TCAPT the m e a n fluid temperature in the collector (°C) TEXT the ambient temperature (°C) p~ the density of the fluid (kg/m 3) dl the mass flow rate in the collector (m3/s) Ts, Te the fluid temperature at the outlet and at the inlet of the collector (°C). F r o m the previous equation, one can deduce that TCAPT(t) = (TcAPT(O) -- T~) exp --
+ Tc
Pu =
FR [ a ( i)7( i ) cos iP D + 0.9a (O)7(O)Pd-- UL A T . . . ]
where T~ is an equilibrium temperature and to a time constant:
where
s(~TExT + ~Pv) + 2 p i C i c l l T ~
c~(i):
solar absorptance of the collector plate r(i): solar transmittance of the transparent covers PD, Pd: solar radiation, direct and diffuse (W/ m 2)
AT:
difference between fluid inlet temperature and outside ambient temperature (°C)
Tc =
SfJ + 2 p l C l d i
where s is the collector area (m 2) To = .
sMi ci
s{3 + 2 p i C l d i
and as a result Pu = 2 p l c l d i ( T c A P T -- Te)
204 S t o r a g e tank. A sensible heat storage is simulated: the storage mass (usually water) is fed by the collected energy and feeds the space or water heating systems; moreover the storage leakage is simulated by a loss coefficient. If the storage temperature exceeds 95 °C (case of water), an extra leakage is created in order to keep the fluid temperature within security limits. As a matter of fact, the temperature of the storage tank is not homogeneous. The effects of natural convection lead to a vertical gradient in the tank which is never negligible. The best solution would be to solve the four equations of the phenomenon, but this is far too complicated from a practical, thermal point of view. A simpler model has been chosen, based on the works of Close [5] and Gutierrez [6]. The tank is divided into N isothermal layers (identical volumes V 1 = V 2 = . . . = V N and temperatures T1 ~< T2 <...~< TN). The fluid coming from the solar collectors (temperature T~, Fig. 4) feeds the upper part of the tank, and transfers solar energy to the layer of same density (layer k such as Tk+ 1 > T~/> Tk) and to the lower ones (temperatures Tl~ - 1,.-., T,), T1), the rate at which energy is transferred being of course proportional to the temperature differences and to the mass flow rates. As a result, there is no heat transfer in the layers k + 1, k + 2 , . . . N and in the other ones ( 1 , 2 , . . . , k -- 1, k) the exchange is proportional to T2 -- T1, T3 -- T2 .... T~ -- TI~ 1, T ~ - - T k.
With the above-mentioned notations, TN is also the temperature at the inlet of the space heating heat exchanger. Let T n be the temperature at the outlet of this heat exchanger ( " r e t u r n " temperature). Tn is a function of TN, of the load and of the heat exchanger efficiency. This " r e t u r n " fluid (Tn) transfers its energy back to the tank in the same way as the solar fluid (Ts) did. The thermal balance of each layer of the tank can be written by considering the heat gains as equal to the heat losses, including the losses to the environment and the supplementary losses. These gains and losses are listed hereafter formally. a) Gains from, and losses to, the collectors. According to the temperature and flow-rate condition, the layer i exchanges the following power Pi:
r,
=_
r,
Fig. 5. Storage
1
tank layers.
Pi = 0
if dl = 0 or T~ < T i
Pi = Y i P l C l d l (Ts -- Ti)
if
d 1
~ O and T i+1 > T~ ~> T i
Pi = Y i P l C l d l (T/+I --Ti)
if d I =~ 0 andT~ ~> Ti+ 1 where Yi and Yi are binary parameters (0 or 1) that identify the layer concerned with the heat exchange, then mathematically (~: Kronecker symbol) : yi = 1
if T i +1 > Ts >~ Ti, and 0 otherwise
N
Yi = ~
y~(1 --5i1)(1 --(~iN)
j=l
b) Space heating system. We can similarly write: Pi=0
if d2 = 0 0 r
TR > Ti
Pi = x i P 2 c 2 d 2 ( T R -- Ti)
if d 2 ¢ 0 and Ti ~> TR > T/_ 1 Pi = X i P 2 c 2 d 2 (Ti 1 -- Ti)
if d2 ¢ 0 a n d T a ~< Ti-1 where xi and X i are binary parameters: if T i ~ TR > T i - 1 ,
xi = l
and 0 otherwise i
Xi = E xj(1--~i])(1
-~1i)
]=1
c) Hot water. P/= 0
if
d 4 =
0
Pi = - - p 2 c 2 d 4 (Ti - - T i - l )
if d4 ¢ 0 where To is the cold water temperature.
205
d) Losses into the environment (temperature
Tac). P i = - - K S i ( T i - - Tac )
where K is a global coefficient (W/m 2 °C) and Si the loss area of the ith layer. The global thermal balance of the ith layer is:
dTi M2C2 --
= Y'(Pi)
dt where EPi is the sum of the various powers Pi listed above. Heating system simulation. From a practical point of view, stored heat is transferred to the rooms by emitters (convectors, heat exchangers, radiators or heating floor); the instantaneous power emitted can be written as a function of Tin t (inside room temperature) and Tc (fluid temperature in the emitters, Tc -~ (TN + TR)/2) P = LRa(Tc -- Tint)
where R a is a global heat exchange coefficient (W/m °C) and L a dimensional parameter. The heating system equation looks like the collector and storage equations:
dVc M3c2
~
= P2c2d2(Tn
-- Ta)--LR~(Tc
- - Tint)
Control system and auxiliary energy. When the energy in the storage tank is depleted, an auxiliary heating system is switched on. Such an on-off regulation poses a problem arising from the step by step procedure of the simulation model. Auxiliary heat is emitted for one step at least (half an hour for instance), and as a consequence, heating energy may exceed the actual need. The room temperature rises and the thermal loads of the next step become smaller. Mathematically, from the thermal balance of the r o o m (T cons = design temperature): Mmem
dTint - VmVm(Tcons dt
- - Tint) + A P
and the energy, AP exp(--GmVm/Mmcmt), will be substracted from the load at the next step. Domestic hot water. The demand for h o t water is simulated three times a day, according to the number of people living in the house: 8.8 1/day per person at 8 a.m. 6.4 1/day per person at midday 22.6 1/day per person at 8 p.m. In the computer program, each person is assumed to have been drawing water for a time step (or two steps in the evening) at the flow rate corresponding to the abovementioned values. If the storage temperature exceeds the design hot water temperature (60 °C), the flow rate is reduced accordingly. General computation method. All the differential equations in the previous sections are put together in order to form a complete set of N + 2 equations with the N + 2 variables being the temperatures of the N layers in the tank, and the temperatures of the solar collectors and the heat distribution equipment. This so-called matrix of coefficients is solved by a Runge-Kutta-Merson method, giving all the above temperatures at each step. Some details of the procedure are given in the following. At each step: - - T s (outlet of the collectors) is compared with each Ti (temperature of each layer), from which xi, Yi, Xi, Yi and the collector flow rate are calculated; if thermal needs (hot water and space heating) as predicted by CLIM are positive, then energy is transferred from the storage tank to the room if the tank temperature exceeds the design temperature; if not, the auxiliary system is activated; - - a l l the coefficients of the temperature matrix are calculated from the previous data, and the matrix equations are solved; the thermal p o w e r rates are calculated; results are recorded by the c o m p u t e r and the procedure starts again as the next step. --
--
- - t h e
We deduce the inside temperature Trot to be
"I
Tint = Tcon s + - 1--exp Gm Vm
( °mm)t t
imCm
" S I M P L I F I E D " MODELS
The excess of losses is only: G m V m (Tin t - - Tcons ) = A p
[
1 --exp
( GmWmt)] --Mmcm
General features As stated in the introduction, the objectives of these models are completely different
206 from the comprehensive methods. The accuracy is n o t so important, but cost, simplicity and public availability become essential. In fact, there is a great diversity among these models: the calculation procedures are often very different, as well as the thermal or meteorological data required. It is even impossible to describe such a m e t h o d without taking a definite one as an example. Nevertheless, before describing two of them, we can bring out some features c o m m o n to all of them. First of all, a big computer program is never necessary, a small computer may be useful in order to save time, but generally all the equations can be solved manually or with a slide-rule (this is why some people say " m a n u a l " models). An iterative m e t h o d is generally used as was the case for the comprehensive models, but the time-lag between two calculation steps is far bigger. This time interval is often one m o n t h because more than one calculation step per m o n t h exceeds the possibilities of a manual yearly based procedure. Therefore, the manual models never take transient phenomena into account, they are based on averaged parameters and on averaged data and the detailed description of the building or of the system is not necessary. This may be advantageous when the system or the building are not known but when orders of magnitude are already necessary. In fact, the best reliability criterion of these simplified models is based on experimental field measurements or on the results of a big computer model. Such a comparison is a reliable test and far better than speculative checkings because the assumptions included in such a manual procedure can never be justified from a theoretical point of view. T h e f-chart m e t h o d
This m e t h o d was developed in 1976-77 at the University of Wisconsin (U.S.A.) by Beckman, Klein and Duffie [7]. It is based on " c h a r t s " where the long term performance of liquid and air systems is correlated with two important dimensionless variables. These charts are presented in several publications in graphical and equation form, and they are referred to as "f-charts" because f is the symbol used by the authors to identify the frac-
tion of the m o n t h l y total heating load supplied by solar energy. The main problem is to identify these two dimensionless variables (two because a sheet of paper is a two-dimensional object). The authors empirically suggest that these variables be: X = A F R ' U L (Tre~ -- T a ) A t / L
Y = AFR'(V~ ) H T N / L where A is the area of the solar collectors (m 2) F a' the collector-heat exchanger efficiency factor, FR -- FR (heat removal efficiency factor) with a " g o o d " heat exchanger UL collector overall energy loss coefficient, (W/m 2 °C) /x t number of seconds in the m o n t h Tret reference temperature (100 °C) Ta m o n t h l y average ambient temperature
(°C) L m o n t h l y heating load (J) H w m o n t h l y average dally radiation incident on the collector surface (J/m 2) N number of days in the m o n t h (r-d) m o n t h l y average transmittanceabsorptance product of the collector. X is the ratio of the collector energy loss at a reference temperature to the heating load of the m o n t h , and Y is the ratio of the energy absorbed on the collector plate surface to the heating load. All of these variables are obtained from collector test results or from m o n t h l y meteorological data, except for HT which must be calculated from the m o n t h l y average daily radiation data on a horizontal surface which are the only radiation data available. This calculation is slightly complicated by the diffuse radiation ratio which is barely known. The correlation between X, Y and f is given in Fig. 6 for liquid-based space heating systems and in Fig. 7 for air systems. This correlation has been developed from the results of a detailed computer simulation (TRNSYS program). To determine f for a given system and for a given month, the first thing to do is to calculate X and Y for the collector and heating load in question. Then f is determined at the intersection of X and Y on the chart, by interpolation if necessary. This is done for each m o n t h of the year. The annual solar contribu-
207
f
j
ature and average ambient temperature. From this chart, a graphic m e t h o d leads to the monthly solar contribution.
f=0
wl<
Creation o f the charts The well-known Hottel equation of a flat plate collector gives the useful gain Qu from the collector (W/m2):
f=0.4
Q3 !laJ
gl "r I
f=O.5
Qu =FR 0
4 8 12 v % REF. COLLECTORLOSS A HEATING LOAD
. f= 09
< .-J a/E ,,,
ii
..,c. -r
f=0.8
2 -
~ ~
f
=
/
/
O
~
f .
:
6
o
z f=0.5 f=0.4
I
:o 0
i
0
t
i
i
i
i
i
I
4
i
i,
i,
ira)]
where G is the global incident radiation (W/m 2) TF the inlet temperature of the fluid (°C) Ta the ambient temperature (°C) (at) the effective absorptance-transmittance product of the transparent cover F a the heat removal efficiency factor UL the heat loss coefficient (W/m 2 °C) provided G exceeds a threshold value (G ~> [UL/(OZT)](T F - - T a ) , and Qu = 0 of G < [UL/ (ar)] (TF -- Ta ). From these two equations, one can derive the average daily energy collected and used (Wh/m 2) :
16
Fig. 6. f-chart for liquid systems.
e e O
[(C~T)G-- UL(T F --
P
Gu = J Qu dt P
L
I
i , , , , ,
,
i
,
,
8 12 X '~ REF. COLLECTORLOSS HEATING LOAD
,
,
,
,
,
~
,
,
,
16
Fig. 7. f-chart for air systems.
tion is the sum of the m o n t h l y energy contribution Q T calculated by the equation Q T =
fL. A great advantage of this method is that the f-charts that are provided by the authors are universal (only for all the U.S. climatic conditions). The great disadvantage is the non-flexibility {other systems, solar assisted heat pumps, etc.) and the preliminary calculation procedure which may be long. The EDF graph method: "Bilan thermique d 'une maison solaire" This m e t h o d has been designed at Electricit~ de France in 1976-76 by Chouard, Michel and Simon [8]. It is based on charts giving the useful solar energy available in a given system configuration as a function of a parameter A T, the difference between fluid temper-
I
where I is the time interval(s) during which G exceeds the threshold value. Let us remark that for a given collector, this value depends only on our parameter A T. Consequently, for a given collector configuration, a given location and a given period (e.g., one month) the energy Gu depends only on the parameter A T. The charts are derived on this basis (A T = T e - - T.). In order to evaluate Gu as a function of AT, let us consider two extreme cases: a sunny day and a cloudy day. In the first case (no clouds all day long), Gu reaches a maximum: Gu = Gu,max = FR f(~r)G dt- UL A TI I
The site latitude and collector tilt being given, the radiation G can be obtained as a function of time t for the typical day of each month. Then the time interval I is easily obtained from the equation: UL G= AT
(~r)
The integral of the above formula is com-
208 p u t e d with the trapezoidal rule: it is the h a t c h e d area in Fig. 8.
2600
2400 Energy
22OO
((~T) -~T
& 2
2000
13-
1800
threshold
E
1600
Solar t i m e
Fig. 8. C o m p u t a t i o n of Gu, max.
In the s e c o n d case ( c l o u d y day), the procedure is identical, e x c e p t t h a t it is n o longer based on G (global radiation) b u t on D (diffuse radiation). T h e useful e n e r g y c o l l e c t e d is G ..... in. T h e real energy Gu is calculated f r o m G ........ G~.mi~ and n/N, the ratio of the average n u m b e r of h o u r s o f bright sunshine to the m a x i m u m possible n u m b e r o f hours o f bright sunshine f o r the same p e r i o d assuming a c o m p l e t e absence o f c l o u d c o v e r (e.g., one m o n t h ) . T h e f o r m u l a is:
v
3
4000/1400
>. oJ
3800/1200
2 3600/1000
3400/800
3200/600
3000/400
2800/200
G,
=~Gu,ma
x +
1-
, Gu,mi n 0
This m e t h o d basically supposes t h a t the n hours o f bright sunshine are u n i f o r m l y dist r i b u t e d over the day which m e a n s t h a t the day is m a d e o f a l t e r n a t e sunshine periods lasting (n/N) dt and c l o u d y periods lasting (1 -n/N) dt. This a s s u m p t i o n is o n l y true if the time period is long e n o u g h (greater t h a n one week is a c c e p t a b l e b u t one m o n t h is preferable). T h e energy G~ is c o m p u t e d r e p e a t e d l y for several values o f the p a r a m e t e r A T which allows curves f o r G~ as a f u n c t i o n o f A T for each m o n t h to be d r a w n for each l o c a t i o n and each s y s t e m c o n f i g u r a t i o n (e.g. Fig. 9). Of course, these charts have been m a d e o n c e and f o r all, and have been published for a lot o f various climatic c o n d i t i o n s . T h e only thing the designer has to do is to m a k e use o f t h e m as described in the n e x t paragraph.
Solar project Let us consider t h r e e main f u n c t i o n s which are integrated in a solar h e a t e d h o m e n a m e l y solar e n e r g y collection, storage and heat distribution. T h e p a r a m e t e r which these functions have in c o m m o n is the t e m p e r a t u r e o f the water (or the air) flowing into the solar
10
20
30
40
50
60
70
(TF - Ta)(°C)
Fig. 9. Single p a n e c o l l e c t o r : m o n t h l y p e r f o r m a n c e profiles for z o n e E ( o r i e n t a t i o n 0 ° , i n c l i n a t i o n 45 ° ).
collectors when the h e a t e x c h a n g e r t e m p e r ature gradient in the storage t a n k is negligible. As a c o n s e q u e n c e , in a p r o j e c t , o n e has t o det e r m i n e u n d e r c o n s t a n t average c o n d i t i o n s the t e m p e r a t u r e o f the w a t e r (or air) and f r o m its k n o w l e d g e the m o n t h l y a m o u n t o f energy really saved due to the use o f solar heat. One can n o t e t h a t in this a p p r o a c h the size of the storage t a n k o n l y influences the transient b e h a v i o u r o f the system (transient means in a building, 1 - 10 h o u r s d u r a t i o n ) the practical sizes o f storage being of the order o f 4 - 5 m 3 o f w a t e r per h o u s e or apartment.
Average storage temperature in a given month. If T~., is the t e m p e r a t u r e o f the water flowing into the solar collectors and Ts the t e m p e r a t u r e o f the storage t a n k t h e n u n d e r c o n s t a n t c o n d i t i o n s TF = T~. T h e useful energy Gu p r e s e n t e d as a f u n c t i o n o f A T = (Tv -- Ta) can be evaluated as a f u n c t i o n o f
209
T~ because T, given by meteorological data can be averaged monthly for space heating evaluation in building. Let us draw for a given month on the same graph as a function of Ts the following characteristics: useful solar energy Eu (translated from Gu) - - s p a c e heating energy load, (kWh)Ec -
I
2000 = W
1600
\
--
Ta)
+
-
~ \ _E~_ _~
\
1400 Febru,~ry
EH) N
\
1200
where Go is the heat losses coefficient in W/m 3 °C V the living space volume in m 3 TNc the internal house temperature below which space heating is needed (depends on G and is equal to 15 °C for G = lW/ m 3 °C) Ta the external average air temperature in °C N the number of days of the month considered EH the domestic hot water energy demand (in France I0 kWh/day on average) solar energy distribution system fitted into the house:
E, ~J
1800
-
"t 2 4 G ° y (TNc E~ = \ iooo
I
2200
Eo\~i~i ,il [
~ 1000 April
,,,
\
800
60o
I: I
Ill',
4oo
t!',l l '
Ill
0
E~-A
,
•,
~
I
ii'
200
~
~-
"
~ ~
~
70
80
hot water dem~
,, i I
10
20TsTS°30 Tsl,40 April
50
60
90
Ts (°C)
Fig. 10. Calculation chart for February, April and August (Nimes house example).
-
heating
EDI = 24Nk I (T, -- Ti) kWh/month
water (at 60 °C)
ED2 ,=
k2N
Ts--Tw 50
where T i is the design temperature
kWh/ month
of the house
( 2 0 °C)
T~ kl k2
the cold water temperature (10 °C) the exchange coefficient (kW/°C) the maximal h o t water demand (10 kWh/day) The functions Eu = f(T~) and ED1 + ED2 have opposite slopes and will intersect at a point at which the useful solar energy (Eo) is greater or less than Ec (see Fig. 10). First case: Eo > Ec. The coordinates Eo and To of the above-defined intersection represent respectively the energy saved in heating demand o f the house for the month considered, and the average temperature of the water in the storage (T~m) and in the solar heat distribution system (Tdm). In this case To = Tam = T~m. Second case: Eo < Ec. The supplementary solar energy heats the storage tank up to an equilibrium average temperature equal to the intersection temperature T~I of the functions Eu = f ( T s ) and Ec = constant. The regulation
of the heat distribution system adapts the temperature of the water (either by an "on and o f f " or a progressive opening of a three way value )to the heating demand (Ts2). As a consequence the average water temperature Tdm of the heating system is different from the one of the storage and is equal to Ts 2. If Tsm = Tsl exceeds 60 °C, solar energy will be able to cover the whole h o t water demand of the house. On the contrary, if T~m = T~I is less than 60 °C, solar energy, though sufficient for the whole energy demand of the house (recall that we are studying the case E o > Ec ) will not be able to heat hot water up to 60 °C, and supplementary energy will be necessary. This supplementary energy is equal to the difference between the h o t water heating requirement (k2N) and the solar energy which heats the water from Tw (10 °C) to T~m (ordinate of the point of the water heating characteristics ED2 = f(Ts) for which the abscissa is Tsm). Examples
In this section, we are going to compare the t w o manual methods we have just been describing. For this purpose, the same project has been studied by the two methods one after another. The problem is to determine the
210
Solarradiation Outside temperature (°C)
Degreedays (basis 16,5°C)
R~f~rence Gj GT (90°,south) (90° , south) kWh/m2.day ~Wh/m2month
Fundamental
Total
AT
load
(°C)
(kWh/month)
Solar
parameters
f
X
Y
energy kWh/month
January
6.5
282
2,25
69,7
93,5
3028
2,93
0,5
0,28
848
February
9,3
192
2,8?
79
90,7
2100
3,71
0,81
0,47
987
March
11,0
170
3,14
97,3
89
1900
4,45
1,12
6,63
1197
April
12,8
122
3,15
94,5
87,2
1410
5,69
1,46
0,73
1028
May
15.2
70
2,75
85,3
84,8
892
9,03
2,08
0,83
740
June
23,4
3
2,58
77,4
76,6
210
33,5
8
1
210
July
23,4
0
2,85
88,3
76,6
186
39,1
10,3
1
186
August
22,1
0
3,18
98,6
77,9
186
39,8
11,4
1
186
September
17,5
20
3,38
101,4
82,5
382
19,8
6,78
1
381
October
13,5
110
3,06
94,9
86,5
1294
6,35
1,59
0,77
996
November
8,9
215
2,31
69,3
91,1
2347
3,57
0,64
0,35
821
December
6,8
274
1.94
60,1
93,2
2948
3,00
0,44
0,24
707
0,491
8287
Year
1458
1016
16883
j
Fig. 11. Typical project by the f-chart method.
Outside
Degree-
Temperature
days
(°C)
I(basis 16,5°C)
January
6,5
282
February
9,3
March
Specific useful
Distribution
Storage
Storage
temperature (Oc)
efficiency
solar
temperature (Oc)
2,791
0,980
25,5
25,5
192
2,145
1,260
26,5
11.0
170
1,751
1,290
April
12,8
122
1,342
May
15,2
70
June
23,4
July
energy (kWh/m 2. day)
Global useful
Solar
energy (kWh/month) load
solar
Contribution (%)
0,760
3028
808
26,7
26,5
0,753
2100
930
44,2
27,0
27,0
0,749
1900
1049
55,2
1,150
27,0
27,0
0,749
1410
905
64.2
0,822
0,709
24,5
26,5
0,753
892
579
64,9
3
0,200
0,200
21,5
56,0
0,538
210
113
53,8
23,4
0
0,171
0,171
--
64,6
0,476
186
89
47,6
August
22,1
0
0,171
0,171
-
81,0
0,356
186
66
35,6
Septembar
17,5
20
0,363
0,363
21,5
84,0
0,336
382
128
33,5
October
13,5
110
1,193
1,135
26,5
38,0
0,669
1294
824
63,6
November
8,9
215
2,235
1,030
25,5
25,5
0,760
2347
822
35,0
December
6,8
274
2,717
0,840
24,5
24,5
0,767
2948
699
23,7
16883
7012
41,5
Months
Year
load
1458
Fig. 12. Typical project by the EDF method.
211 solar c o n t r i b u t i o n o f the f o l l o w i n g s y s t e m o n a y e a r l y basis (space h e a t i n g + d o m e s t i c h o t water). Building: o n e - f a m i l y h o u s e , l o c a t e d near Avignon ( F r a n c e ) . Latitude: 4 3 ° N . Living area: 140 m 2 Loads: t h e r m a l loss c o e f f i c i e n t : G = 1.2 W/m 3 °C, degree-days: as listed in Figs. 11 and 12, d o m e s t i c h o t water: 150 1/day ( t e m p e r a t u r e s 14 - 50 °C) Collectors: area: 35 m2; s o u t h facing; tilt angle: 90 ° (vertical); d o u b l e glazing: F R U L = 4 W/m 2 °C, F ~ (sT) -- 0.65 Storage tank: c a p a c i t y : 3 0 0 0 1 (85.7 1/m 2) D i s t r i b u t i o n heating floor, emission coeffiequipment: c i e n t 300 W/°C Auxiliary electric h e a t i n g b y individual energy: convectors T h e results are s h o w n in Figs. 11 and 12.
CONCLUSIONS It was n o t possible in such a p a p e r to s h o w all t h e details o f the solar models, n o r even to describe t h e i r t h e o r e t i c a l and practical procem o n s t r a t e t h a t solar m o d e l s play a great rSle
dures in an extensive way. We p r e f e r r e d t o dein all solar research programs t h a n k s t o their flexibility, and t h a t u n d e r n o c i r c u m s t a n c e s m u s t c o m p r e h e n s i v e and simplified m o d e l s be c o m p a r e d . Each o f t h e m has its o w n scope of applications.
REFERENCES 1 G. P. Mitalas and D. G. Stephenson, Cooling load calculations by thermal response factor method, A S H R A E Trans., 73 (1967) part I, paper 2018. 2 R. S. Varga, Matrice iterative analysis, Series in Automatic Computation, Prentice-Hall, 1962. 3 H. C. Hottel and A. Whillier, Evaluation of flatplate collector performance, in Trans. of the Conference on the use of Solar Energy, Vol. II, Thermal Processes, Univ. of Arizona (1955). 4 E. Aranovitch, The joint solar collector testing programme of the European community, in Proc. ISES (UK Section) Meeting (April 1977). 5 D. J. Close, A design approach for solar processes, Solar energy, 11 (2)(1967). 6 G. Gutierrez et al., Simulation of forced circulation water heaters: effects of auxiliary energy supply, load type and storage capacity, Solar Energy, 15 (1974) 287 - 298. 7 W. A. Beckman, S. A. Klein and J. A. Duffle, Solar Heating Design by the f-Chart Method, Wiley, NewYork, 1977. 8 Ph. Chouard, H. Michel and M. F. Simon, Bilan thermique d'une maison solaire -- Mdthode de calcul rapide, Eyrolles, Paris, 1977.
197
Solar Heating System Design: The Part of Simulation Models P. C H O U A R D
E.D~F. Les Renardidres, F-77250 Moret/Loing (France)
INTRODUCTION
Models, why ? All the branches of physics are deeply concerned with simulation: to simulate by speculative ways the behaviour of physical systems is a very powerful tool at the disposal of research workers and engineers in all the fields of their activities. The thermal evaluation of buildings is specially concerned with this tool because all the others methods in competition with simulation such as full-size field-experiments or laboratory tests are always very expensive and restricted: expensive because of high cost of the p r o t o t y p e s themselves and of the instrumentation needed for proper checking on performances, and restricted because it is impossible to scan all the parameters involved and to get a sufficient knowledge of them in their full range o f variation. Today, a lot of applications of solar energy systems for space heating, cooling, and domestic h o t water are available; this is fully consistent with these thermal simulation methods, because solar systems being basically driven by transient phenomena often require accurate methods for the prediction of their thermal performance or life-cycle cost analysis. As a matter of fact, the thermal models will be able to: a) get a better appreciation of the physical p h e n o m e n a involved and to evaluate their o u t c o m e from a practical point of view; b) predict the behaviour of materials, i.e. to estimate their reliability and duration; c) size the systems and their c o m p o n e n t s as a function of constraints and objectives. Moreover the research engineer also needs such a model when, having to conceive new techniques or new applications, he has to
know the sensibility of a given parameter from a technical and/or economical point of view.
Various categories of models In accordance with the numerous features listed above, a great number of solar models n o w exist; b u t the objectives, the accuracy, the costs, and the procedure of each model are very different one from another, and a classification of them seems necessary. A good criterion might be their relative complexity according to the c o m p u t e r size needed, or the expected cost of a typical testproblem. We shall in fact distinguish two categories of solar models: the "simplified" models and the "comprehensive" ones. The former are specially designed for consulting engineers at the stage of preliminary design, who have to get a first idea of a given system configuration quickly and at a low cost. At this stage, it is completely impossible to use big computer or comprehensive programs far t o o expensive to compete with the simplified ones when the components of the system or the configuration of the building are not yet defined or even when the "solar choice" has n o t yet been done by the customer! Several models belonging to this category already exist in Europe or in the U.S.A., but they cannot be compared with the comprehensive ones from an accuracy point of view. These comprehensive models are in fact used by research people because they are adapted to fundamental studies or to advanced projects, due to their high degree of accuracy. On the other hand, the calculations involved are so complex that a more or less sophisticated c o m p u t e r is always nec~.ssary in that case.
198 First, we shall describe in this paper the comprehensive models taking the computer program " C L I M A S O L " designed by Electricit~ de France as an example; in a second step, we shall describe the principles of the simplified simulations and take the U.S. f-chart m e t h o d and the French "Bilan therm i q u e " m e t h o d as examples.
COMPREHENSIVE MODELS
General features All the comprehensive models show c o m m o n features because there are often very few ways to reach the above-mentioned objectives. First of all, the computation procedure is always a step-by-step one, which means that a recurring m e t h o d is used: all the equations are solved at a given m o m e n t and these results are used as inputs for the next calculations at the next moment. As a matter of fact, building materials always exhibit a great inertia, therefore the lag between two calculation steps can be up to several hours, according to the materials and to the degree of accuracy needed. A one-hour period seems to be quite satisfactory in most cases. As a consequence of the objectives, specially of the transient phenomena requirements, it is necessary to go into all the details of the description of the building itself. In order to take all the heat transfers between rooms, or between inside and outside, at every m o m e n t into account one has to describe the location of the building, its environment in relation to other buildings, the room shapes, the physical properties of materials, and so on. A good knowledge of the meteorological environment is also necessary, and the accuracy of its description depends on the accuracy of the calculation procedure itself. As an example, if the calculations are repeated hourly, an hourly description of the climate will be necessary. If such climatic data are not available, an interpolation will be made. Internal heat gains are an essential parameter in m o d e m , well-insulated buildings. Therefore, the comprehensive models attach the greatest importance to them. Lighting (power dissipated, duration, type of lamps), appliances, and people (number, activity on an hourly basis) are terms of the thermal
balance of the building and appear in the equations, as well as the solar gains through the windows and even through the walls. Once all these data are known, the computation itself becomes possible. From a theoretical point of view, and except for the heat transfer through the walls, the algorithms are quite simple. The question is to take all the gains, all the transfers and all the above parameters into account and to write that all the instantaneous energy flows (positive and negative) are balanced at every moment. The difficulties come from the number of iterations (720 hours within a single m o n t h for instance) and from the solution of the transient equations for the heat transfer through walls and floors. Some methods will be described in the next section. For a large number of iterations a more powerful computer will be needed. A great flexibility does exist owing to the computer. The simulation can be done over long periods, often one year when a typical meteorological year is considered, but it could be still longer. Moreover, as all the variables are calculated (temperatures, energy flows, etc.), it is possible to get a good knowledge of the thermal response of the building in addition to its energy consumption. In fact, it is often a question of printed material because it may be difficult to choose the proper variable to be printed! A lot of additional possibilities still exist. Some models include financial computations, life cycle cost analysis, return on investment, and so on. In fact, everything is possible, and the biggest question looks like a dilemma in connection with the computer: if the computer program is too sophisticated, it cannot be handled by engineers or scientists, it becomes an informatical concept without any practical use. This might be a danger.
Heat transfer models Without giving too many details (it would require a few hundred pages), we are just going to set the problem and to describe its solution by two usual algorithms.
The problem o f the wall For instance, let us consider a wall between two mediums at different air temperatures: Tai (inside) and Tao (outside). The temperature T (x, t) in the wall is a solution for the
199
unidimensional heat transmission equation at every m o m e n t a
~2T
- ~x
-
2
8T ~t
where a: diffusivity of the material. The two limiting conditions of this equation result from the thermal equilibrium of the surfaces of the wall (inside and outside) inside: = hi[Tai -- T(1, t)] + ¢R outside: = ho[Tao -- T(0, t)] + ¢A
because we just need to know the inside temperature of the wall at every moment. There are t w o methods adapted to this problem: the finite difference method and the response factor method. The finite difference m e t h o d This is a computer based method. The above equation is discretized in both space and time, by a substitution of finite differences for derivates. The wall is divided into N elementary layers, and with the above equation, the temperature T (x, t + At) of each junction is calculated at time t + A t from the temperatures of the next junctions at time t and/or t + At. Mathematically, the equation ~2T ~T a -~x 2
where X:
conductivity of the wall (W/m °C) inside and outside convective transfer coefficient (W/m e °C) long wave radiation between the in¢R: side surface and its environment radiation absorbed by the outside CA: surface of the wall. Let us remark that if the wall is heterogeneous (n layers with different materials), the problem is identical. Instead of two, there are n + 1 limit conditions all looking like this one
is written like this
hi, ho :
~ki ~X
i=li
---- ~ki+ 1
} ~X/Xi+l
The initial time condition can correspond to a linear distribution of the temperatures in the wall, from T (0, t) up to T (1, t). The problem is to solve the equation 52T ~T
a
-~x 2
a
T(x + Ax, t ) - - 2 T ( x , t) + T ( x -- Ax, t)
=
(Ax) 2
T(x, t + A t) -- T ( x , t) At
which allows T (x, t + A t) to be evaluated as a function of T (x + Ax, t), T (x, t) and T (x -Ax, t); this is the so-called "explicit" method. The equation a
---- 0
~t
~2T
~T
8x 2
~t
could also be written like this a
T(x +Ax, t+At)--2T(x,t+
At)+T(x--AX,
(~x) 2 = T ( x , t + A t) -- T ( x , t)
- --
At
~t '
Outside
Wall
x
T=o ~.
Inside
i. OT
ax
0
Fig. 1. Heat transfer through a wall.
r.i
and thus T (x, t + A t) is evaluated as a function of T (x -- At, t + At), T ( x , t) and T ( x + Ax, t + At); this is the so-called "implicit" method. The former is quite satisfactory b u t the inequation a~t
1
(~x) 2
2
must be fulfilled for the convergence of the algorithm; therefore, the time interval At must be very small (a few minutes). The latter m e t h o d is always convergent, but it is time-
t+At
200
consuming and requires a bigger area in the computer core.
TCmpe,atu,e
\-/ The response factor method This m e t h o d is based on the thermal response of a wall to an elementary stimulus from a given signal. If we admit that all the equations of the systems are linear, and that the thermal properties of the materials are hot time-changing, we find that the response of the wall to any signal is the sum of the elementary responses we already know. Therefore, there are three steps in the procedure. We must: -- split the excitation function into elementary impulses (for instance into temperature steps or triangular functions); --calculate the response of the system (the wall) to each elementary stimulus; -- add up all these responses. As a matter of fact, every time-dependent function can be split into a series of numbers equal to the values of this function at the moments 0, At, 2 At, 3 At .... Thus, a series of triangular functions replaces the actual thermal function (see Fig. 2). Other possibilities exist, but the work of Mitalas [1] demonstrated t h a t the triangular function is best adapted to such a purpose. :;';:::r::o'0
The,real Triangutar functlons
.~r
Fig. 2. Splitting triangles.
2 ~t
3 .~t
4 .~t
of a c o n t i n u o u s
function
into
Therefore, it is only necessary to know the response of a given wall to a triangular stimulus from the outside air temperature. These responses at each time interval are the so-called response factors of the wall (see Fig. 3). Let us take as an example a variation of the outside temperature which generates a variation in the inner temperature distribution of the wall: T (x, t). The knowledge of this distribution leads to the expression for the energy flow crossing the inside surface of the wall
¢ (1,0 =-x\ x=l
4, /
..............
/
~
~ T,m~ ~t
? zt
3 ~t
4 ~r
5 at
6 ~t
Fig. 3. Definition of response factors.
The distribution itself can be evaluated by the finite difference method. It can also be measured if thermal data is available. From a theoretical point of view, the number of thermal response factors of a wall is infinite (as many as the time intervals n.At). As a matter of fact, the response factors become negligible after 20 hours in most practical cases.
Review o f the main existing computer programs It is impossible to describe all of the comprehensive models that exist all over the world. From a general point of view, these models belong either to the "response f a c t o r " category or to the "finite difference" category, but the i n p u t / o u t p u t procedures are often very different from one another according to the objectives of their designers. In fact, we would like only to demonstrate in this section t h a t a great number of models are operative and that very few comparative studies have been done. Therefore our review looks more or less like a list. U.S. models Among the numerous computer programs available in the United States, we shall just mention the most famous ones. The BLAST (Building Loads Analysis and System Thermodynamics) program is a comprehensive set of subprograms for predicting energy consumption in buildings. It has been designed at the U.S. Army Construction Engineering Research Laboratory (CERL). Each of the three major subprograms: --computes the space load on an hourly basis, using the ASHRAE methods; --simulates the air distribution system and calculates the h o t water, chilled water and electric energy demands; - - simulates the central plant (boilers, chillers, solar systems) and computes energy consumption and life cycle cost. The main
201
feature of BLAST is that a good simulation of the heating or cooling equipment is done on an hourly basis. TRNSYS has been designed at the University of Wisconsin by Klein et al. in 1973". It can simulate the dynamic behaviour of solar transient systems. It consists of components models (for collectors, controls, storage tanks, heat exchangers, furnaces, building loads) and an executive routine: the user selects both the components and their design parameters, and specifies the way in which the components are interconnected. The building is n o t simulated, the user only specifies a global loss coefficient (UA or GV in W/°C); on the contrary, all the details of the solar system are described and simulated. The DOE-1 (formerly Cal-ERDA) is sponsored by the U.S. Department of Energy. It is another program available to the public which simulates the thermal load of a building using the response factors methods. The authors** (Lokmanhekim et al.) simplified the calculation algorithms and as a result the program is very fast and its running cost rather low. European models Several models are available of which we will only mention the most detailed ones. In Belgium, the University of Li6ge has designed several computer programs called LPB1, LPB3, LPB4, each of them exhibiting an increasing degree of complexity in building simulation. In France, the COSTIC (Comit6 Scientifique et Technique des Industries du Chauffage) uses a detailed program which simulates the building by the response factor method and the solar system by modular subroutines. This program has been widely used in order to generate "rules of t h u m b " for domestic water heating. Unfortunately, this program is not publicly available and not much literature has been published.
*S.A. Klein e t al., TRNSYS: A transient system simulation program, Solar Energy Laboratory, University of Wisconsin at Madison, June 1979. **M. Lokmanhekim et al., DOE-2, Building Energy Analysis Group, Energy and Environment Division, Lawrence Berkeley Laboratory, Berkeley, CA, Feb. 1979.
The INSA (Institut National des Sciences Appliqu6es) high school in Lyon has done a very useful synthetic j o b in the field of comprehensive programs. Their model is based on the response factor method. We also mention here the CLIMASOL program, designed at Electricit~ de France, which will be described in the next section. E D F computer program: C L I M A S O L This program basically consists of two sections that can be described separately. The first section (CLIM) computes the instantaneous thermal heating and/or cooling loads of a given building taking the transient phenomena into account. The second section (MASOL) is a detailed model of a solar heating device which simulates the operation of a solar system in order to meet the loads previously calculated. As a result, MASOL gives the energy balance of the system (solar energy used and supplementary energy if any) on a daily, monthly or yearly basis. Thermal loads: CLIM program The general principles of the program are listed hereafter. First of all one describes the building, its climatic location and utilization. The building is divided into " m o d u l e s " which are (by definition) a group of rooms having the same thermal characteristics (inside and outside temperatures, orientation, utilization, etc.). From all this data, the program computes, step by step, the energy requirements in terms of power and power consumption in each module. Data acquisition. The computer reads all the useful data via punched cards of which the main ones are: TITRE (title): this card describes the geographic location of the building and the weather data to be considered. It also defines the category of the building (residential, commercial, etc.). MATI~RIAUX & PAROIS (materials and walls): these cards correspond to the frame of the buildings and give the composition of each wall, window, floor, ceiling, and so on. They describe the physical properties of the building materials (thermal conductivity, density and specific heat) as well as the possible composition of the walls (thickness of each layer, tilt angle from horizontal, surface
202 resistance, shading coefficient of the windowblinds, if any, etc.). R E N O U V E L L E M E N T D ' A I R (air change rate) : this is a card for the fresh air mass flow rate into each module, its humidity and temperature and the electric consumption of fans (if any). It is possible to define several air change ratios as a function of time, according for example to the number of people living in the building. APPORTS (free heat): it is possible to define hourly the rates which heat and moisture are separately given off by people, lighting, and/or appliances in each module. These heat gains play a great rSle in modern, wellinsulated buildings and this is why they are described in detail. PUISSANCES INSTALLI~ES (reference power of the heating or cooling equipment): this card represents the thermal maximum power of the available equipment (fundamental data as far as transient phenomena are concerned). MODULES: detailed description of the rooms (modules) themselves, dimensions, characteristics of the walls (from the MATl~RIAUX cards) and inside design temperatures (minimum and maximum); the contiguous modules are also identified by this card.
Weather data and solar radiation. Data on o u t d o o r weather conditions (dry bulb temperature, relative humidity and solar coefficients) have been stored in the computer m e m o r y for 26 stations in France for a typical year. As far as solar radiation is concerned, the solar position (altitude and azimuth from South) is c o m p u t e d at every m o m e n t taking into account such data as latitude, day of the year, hour of the day, spherical coordinates of the radiation receiver (tilt from horizontal, and deviation from South) and so on. These calculations together with the solar coefficients coming from the meteorological data allow an accurate estimation of the actual incident radiation to be considered for the energy balance of the building. This procedure which is based both on calculation and on meteorological data makes up for the shortage of relevant data (very few weather stations make measurements of radiation on tilted surfaces on an hourly basis). Computation o f thermal phenomena. As in all the models belonging to this category, the
problem is to calculate the air and wall temperature in each module (room) as a function of: the internal heat gains listed above; the design (inside and outside) conditions; and the temperature of the contiguous modules. In order to solve the heat transfer equations in each wall, a finite difference computation method has been chosen; the matrix representing the system is inversed by a ShermanMorrison [2] method. Simplifying assumptions have been made, the justification of which would be quite a lengthy paper itself. These assumptions concern: - - t h e inside temperature which is supposed to be homogeneous in the living space area; - - t h e heat flux which is supposed to be normal to the wall surfaces (heat losses through "thermal bridges" are calculated separately). - - t h e solar radiation transmitted by the window panes which is assumed to be evenly distributed on the floor surface (i.e., no shadow limit is calculated). From all these data, algorithms and assumptions, CLIM can c o m p u t e the thermal heating or cooling loads and the inside temperature in each module twice an hour (i.e., 24 X365 X 2 = 17 520 iteration steps are necessary to simulate a whole year). These loads will be interpreted as input data by the next subroutine: MASOL.
Solar system simulation: M A S O L program Most solar energy systems can be described as a group of three major interconnected subsystems namely collection of solar heat, storage of energy and distribution devices. Basically, MASOL consists of a detailed model of these three sub-systems: each subsystem includes its own inertia, heat losses and heat gains. General description. Figure 4 shows the block diagram of the system. Energy flows are represented by arrows and by the mass flow rates of the fluids. The distribution system includes a recirculated fluid loop (for space heating) and a non-recirculated loop (domestic hot water for example). These loops can be activated together or separately according to the actual system. The control subroutines are also copied closely from the actual control equipment:
203 FR:
d,s,,,b,~,on
collector heat removal efficiency factor UL: collector energy loss coefficient (W/ m 2 °C) On the other hand, the CLIM program computes the solar radiation transmitted by a given surface
Pv = r(i)PD COS i + 0.9 7(o)P d Col~ection
SIo,a~e
D,str,bution
Fig. 4. Schematic diagram of the solar system simulated by MASOL.
-- if the temperature at the outlet of the solar collectors exceeds the lowest temperature of the storage tank, solar heat is transferred from the collectors to the storage tank (mass flow rate at its typical value in the collector-tank loop) ; if the storage temperature exceeds a threshold value and if the rooms need that (positive load as c o m p u t e d by CLIM), stored heat is transferred to the distribution system (flow rate at its typical value in the distribution loop); -- if the load cannot be met by the solar distribution system, the supplementary energy system is activated. -
Thus, we can write (assuming a(i) = a, whatever i is)
Pu = FR (aPv -- ~A T) where Pv comes from CLIM, a is the plate absorptance, and ~ the collector overall energy loss coefficient ( W / m 2 °C); these numbers can best be determined by collector tests in which instantaneous efficiency is measured under specified conditions (e.g. [4] ). The general thermal balance of the collector will be written as dTcAPT
-
H e a t collection. The theory relating instantaneous collector performance to meteorological and operating conditions is well-known from the studies of Hottel et al. [3] (during the years 1942 to 1959). We can thus express the rate at which useful energy is collected as a function of solar radiation, of outside ambient temperature and of fluid temperature. The equation is
M1C1
dt
= ~Pv --/3(TcAPT -- TEXT)
--P 1Cidi(T~
-- Te)
where M i is the equivalent mass capacity of the collector, fluid + material (kg/m 2) C i the heat capacity of the fluid (J/kg °C) TCAPT the m e a n fluid temperature in the collector (°C) TEXT the ambient temperature (°C) p~ the density of the fluid (kg/m 3) dl the mass flow rate in the collector (m3/s) Ts, Te the fluid temperature at the outlet and at the inlet of the collector (°C). F r o m the previous equation, one can deduce that TCAPT(t) = (TcAPT(O) -- T~) exp --
+ Tc
Pu =
FR [ a ( i)7( i ) cos iP D + 0.9a (O)7(O)Pd-- UL A T . . . ]
where T~ is an equilibrium temperature and to a time constant:
where
s(~TExT + ~Pv) + 2 p i C i c l l T ~
c~(i):
solar absorptance of the collector plate r(i): solar transmittance of the transparent covers PD, Pd: solar radiation, direct and diffuse (W/ m 2)
AT:
difference between fluid inlet temperature and outside ambient temperature (°C)
Tc =
SfJ + 2 p l C l d i
where s is the collector area (m 2) To = .
sMi ci
s{3 + 2 p i C l d i
and as a result Pu = 2 p l c l d i ( T c A P T -- Te)
204 S t o r a g e tank. A sensible heat storage is simulated: the storage mass (usually water) is fed by the collected energy and feeds the space or water heating systems; moreover the storage leakage is simulated by a loss coefficient. If the storage temperature exceeds 95 °C (case of water), an extra leakage is created in order to keep the fluid temperature within security limits. As a matter of fact, the temperature of the storage tank is not homogeneous. The effects of natural convection lead to a vertical gradient in the tank which is never negligible. The best solution would be to solve the four equations of the phenomenon, but this is far too complicated from a practical, thermal point of view. A simpler model has been chosen, based on the works of Close [5] and Gutierrez [6]. The tank is divided into N isothermal layers (identical volumes V 1 = V 2 = . . . = V N and temperatures T1 ~< T2 <...~< TN). The fluid coming from the solar collectors (temperature T~, Fig. 4) feeds the upper part of the tank, and transfers solar energy to the layer of same density (layer k such as Tk+ 1 > T~/> Tk) and to the lower ones (temperatures Tl~ - 1,.-., T,), T1), the rate at which energy is transferred being of course proportional to the temperature differences and to the mass flow rates. As a result, there is no heat transfer in the layers k + 1, k + 2 , . . . N and in the other ones ( 1 , 2 , . . . , k -- 1, k) the exchange is proportional to T2 -- T1, T3 -- T2 .... T~ -- TI~ 1, T ~ - - T k.
With the above-mentioned notations, TN is also the temperature at the inlet of the space heating heat exchanger. Let T n be the temperature at the outlet of this heat exchanger ( " r e t u r n " temperature). Tn is a function of TN, of the load and of the heat exchanger efficiency. This " r e t u r n " fluid (Tn) transfers its energy back to the tank in the same way as the solar fluid (Ts) did. The thermal balance of each layer of the tank can be written by considering the heat gains as equal to the heat losses, including the losses to the environment and the supplementary losses. These gains and losses are listed hereafter formally. a) Gains from, and losses to, the collectors. According to the temperature and flow-rate condition, the layer i exchanges the following power Pi:
r,
=_
r,
Fig. 5. Storage
1
tank layers.
Pi = 0
if dl = 0 or T~ < T i
Pi = Y i P l C l d l (Ts -- Ti)
if
d 1
~ O and T i+1 > T~ ~> T i
Pi = Y i P l C l d l (T/+I --Ti)
if d I =~ 0 andT~ ~> Ti+ 1 where Yi and Yi are binary parameters (0 or 1) that identify the layer concerned with the heat exchange, then mathematically (~: Kronecker symbol) : yi = 1
if T i +1 > Ts >~ Ti, and 0 otherwise
N
Yi = ~
y~(1 --5i1)(1 --(~iN)
j=l
b) Space heating system. We can similarly write: Pi=0
if d2 = 0 0 r
TR > Ti
Pi = x i P 2 c 2 d 2 ( T R -- Ti)
if d 2 ¢ 0 and Ti ~> TR > T/_ 1 Pi = X i P 2 c 2 d 2 (Ti 1 -- Ti)
if d2 ¢ 0 a n d T a ~< Ti-1 where xi and X i are binary parameters: if T i ~ TR > T i - 1 ,
xi = l
and 0 otherwise i
Xi = E xj(1--~i])(1
-~1i)
]=1
c) Hot water. P/= 0
if
d 4 =
0
Pi = - - p 2 c 2 d 4 (Ti - - T i - l )
if d4 ¢ 0 where To is the cold water temperature.
205
d) Losses into the environment (temperature
Tac). P i = - - K S i ( T i - - Tac )
where K is a global coefficient (W/m 2 °C) and Si the loss area of the ith layer. The global thermal balance of the ith layer is:
dTi M2C2 --
= Y'(Pi)
dt where EPi is the sum of the various powers Pi listed above. Heating system simulation. From a practical point of view, stored heat is transferred to the rooms by emitters (convectors, heat exchangers, radiators or heating floor); the instantaneous power emitted can be written as a function of Tin t (inside room temperature) and Tc (fluid temperature in the emitters, Tc -~ (TN + TR)/2) P = LRa(Tc -- Tint)
where R a is a global heat exchange coefficient (W/m °C) and L a dimensional parameter. The heating system equation looks like the collector and storage equations:
dVc M3c2
~
= P2c2d2(Tn
-- Ta)--LR~(Tc
- - Tint)
Control system and auxiliary energy. When the energy in the storage tank is depleted, an auxiliary heating system is switched on. Such an on-off regulation poses a problem arising from the step by step procedure of the simulation model. Auxiliary heat is emitted for one step at least (half an hour for instance), and as a consequence, heating energy may exceed the actual need. The room temperature rises and the thermal loads of the next step become smaller. Mathematically, from the thermal balance of the r o o m (T cons = design temperature): Mmem
dTint - VmVm(Tcons dt
- - Tint) + A P
and the energy, AP exp(--GmVm/Mmcmt), will be substracted from the load at the next step. Domestic hot water. The demand for h o t water is simulated three times a day, according to the number of people living in the house: 8.8 1/day per person at 8 a.m. 6.4 1/day per person at midday 22.6 1/day per person at 8 p.m. In the computer program, each person is assumed to have been drawing water for a time step (or two steps in the evening) at the flow rate corresponding to the abovementioned values. If the storage temperature exceeds the design hot water temperature (60 °C), the flow rate is reduced accordingly. General computation method. All the differential equations in the previous sections are put together in order to form a complete set of N + 2 equations with the N + 2 variables being the temperatures of the N layers in the tank, and the temperatures of the solar collectors and the heat distribution equipment. This so-called matrix of coefficients is solved by a Runge-Kutta-Merson method, giving all the above temperatures at each step. Some details of the procedure are given in the following. At each step: - - T s (outlet of the collectors) is compared with each Ti (temperature of each layer), from which xi, Yi, Xi, Yi and the collector flow rate are calculated; if thermal needs (hot water and space heating) as predicted by CLIM are positive, then energy is transferred from the storage tank to the room if the tank temperature exceeds the design temperature; if not, the auxiliary system is activated; - - a l l the coefficients of the temperature matrix are calculated from the previous data, and the matrix equations are solved; the thermal p o w e r rates are calculated; results are recorded by the c o m p u t e r and the procedure starts again as the next step. --
--
- - t h e
We deduce the inside temperature Trot to be
"I
Tint = Tcon s + - 1--exp Gm Vm
( °mm)t t
imCm
" S I M P L I F I E D " MODELS
The excess of losses is only: G m V m (Tin t - - Tcons ) = A p
[
1 --exp
( GmWmt)] --Mmcm
General features As stated in the introduction, the objectives of these models are completely different
206 from the comprehensive methods. The accuracy is n o t so important, but cost, simplicity and public availability become essential. In fact, there is a great diversity among these models: the calculation procedures are often very different, as well as the thermal or meteorological data required. It is even impossible to describe such a m e t h o d without taking a definite one as an example. Nevertheless, before describing two of them, we can bring out some features c o m m o n to all of them. First of all, a big computer program is never necessary, a small computer may be useful in order to save time, but generally all the equations can be solved manually or with a slide-rule (this is why some people say " m a n u a l " models). An iterative m e t h o d is generally used as was the case for the comprehensive models, but the time-lag between two calculation steps is far bigger. This time interval is often one m o n t h because more than one calculation step per m o n t h exceeds the possibilities of a manual yearly based procedure. Therefore, the manual models never take transient phenomena into account, they are based on averaged parameters and on averaged data and the detailed description of the building or of the system is not necessary. This may be advantageous when the system or the building are not known but when orders of magnitude are already necessary. In fact, the best reliability criterion of these simplified models is based on experimental field measurements or on the results of a big computer model. Such a comparison is a reliable test and far better than speculative checkings because the assumptions included in such a manual procedure can never be justified from a theoretical point of view. T h e f-chart m e t h o d
This m e t h o d was developed in 1976-77 at the University of Wisconsin (U.S.A.) by Beckman, Klein and Duffie [7]. It is based on " c h a r t s " where the long term performance of liquid and air systems is correlated with two important dimensionless variables. These charts are presented in several publications in graphical and equation form, and they are referred to as "f-charts" because f is the symbol used by the authors to identify the frac-
tion of the m o n t h l y total heating load supplied by solar energy. The main problem is to identify these two dimensionless variables (two because a sheet of paper is a two-dimensional object). The authors empirically suggest that these variables be: X = A F R ' U L (Tre~ -- T a ) A t / L
Y = AFR'(V~ ) H T N / L where A is the area of the solar collectors (m 2) F a' the collector-heat exchanger efficiency factor, FR -- FR (heat removal efficiency factor) with a " g o o d " heat exchanger UL collector overall energy loss coefficient, (W/m 2 °C) /x t number of seconds in the m o n t h Tret reference temperature (100 °C) Ta m o n t h l y average ambient temperature
(°C) L m o n t h l y heating load (J) H w m o n t h l y average dally radiation incident on the collector surface (J/m 2) N number of days in the m o n t h (r-d) m o n t h l y average transmittanceabsorptance product of the collector. X is the ratio of the collector energy loss at a reference temperature to the heating load of the m o n t h , and Y is the ratio of the energy absorbed on the collector plate surface to the heating load. All of these variables are obtained from collector test results or from m o n t h l y meteorological data, except for HT which must be calculated from the m o n t h l y average daily radiation data on a horizontal surface which are the only radiation data available. This calculation is slightly complicated by the diffuse radiation ratio which is barely known. The correlation between X, Y and f is given in Fig. 6 for liquid-based space heating systems and in Fig. 7 for air systems. This correlation has been developed from the results of a detailed computer simulation (TRNSYS program). To determine f for a given system and for a given month, the first thing to do is to calculate X and Y for the collector and heating load in question. Then f is determined at the intersection of X and Y on the chart, by interpolation if necessary. This is done for each m o n t h of the year. The annual solar contribu-
207
f
j
ature and average ambient temperature. From this chart, a graphic m e t h o d leads to the monthly solar contribution.
f=0
wl<
Creation o f the charts The well-known Hottel equation of a flat plate collector gives the useful gain Qu from the collector (W/m2):
f=0.4
Q3 !laJ
gl "r I
f=O.5
Qu =FR 0
4 8 12 v % REF. COLLECTORLOSS A HEATING LOAD
. f= 09
< .-J a/E ,,,
ii
..,c. -r
f=0.8
2 -
~ ~
f
=
/
/
O
~
f .
:
6
o
z f=0.5 f=0.4
I
:o 0
i
0
t
i
i
i
i
i
I
4
i
i,
i,
ira)]
where G is the global incident radiation (W/m 2) TF the inlet temperature of the fluid (°C) Ta the ambient temperature (°C) (at) the effective absorptance-transmittance product of the transparent cover F a the heat removal efficiency factor UL the heat loss coefficient (W/m 2 °C) provided G exceeds a threshold value (G ~> [UL/(OZT)](T F - - T a ) , and Qu = 0 of G < [UL/ (ar)] (TF -- Ta ). From these two equations, one can derive the average daily energy collected and used (Wh/m 2) :
16
Fig. 6. f-chart for liquid systems.
e e O
[(C~T)G-- UL(T F --
P
Gu = J Qu dt P
L
I
i , , , , ,
,
i
,
,
8 12 X '~ REF. COLLECTORLOSS HEATING LOAD
,
,
,
,
,
~
,
,
,
16
Fig. 7. f-chart for air systems.
tion is the sum of the m o n t h l y energy contribution Q T calculated by the equation Q T =
fL. A great advantage of this method is that the f-charts that are provided by the authors are universal (only for all the U.S. climatic conditions). The great disadvantage is the non-flexibility {other systems, solar assisted heat pumps, etc.) and the preliminary calculation procedure which may be long. The EDF graph method: "Bilan thermique d 'une maison solaire" This m e t h o d has been designed at Electricit~ de France in 1976-76 by Chouard, Michel and Simon [8]. It is based on charts giving the useful solar energy available in a given system configuration as a function of a parameter A T, the difference between fluid temper-
I
where I is the time interval(s) during which G exceeds the threshold value. Let us remark that for a given collector, this value depends only on our parameter A T. Consequently, for a given collector configuration, a given location and a given period (e.g., one month) the energy Gu depends only on the parameter A T. The charts are derived on this basis (A T = T e - - T.). In order to evaluate Gu as a function of AT, let us consider two extreme cases: a sunny day and a cloudy day. In the first case (no clouds all day long), Gu reaches a maximum: Gu = Gu,max = FR f(~r)G dt- UL A TI I
The site latitude and collector tilt being given, the radiation G can be obtained as a function of time t for the typical day of each month. Then the time interval I is easily obtained from the equation: UL G= AT
(~r)
The integral of the above formula is com-
208 p u t e d with the trapezoidal rule: it is the h a t c h e d area in Fig. 8.
2600
2400 Energy
22OO
((~T) -~T
& 2
2000
13-
1800
threshold
E
1600
Solar t i m e
Fig. 8. C o m p u t a t i o n of Gu, max.
In the s e c o n d case ( c l o u d y day), the procedure is identical, e x c e p t t h a t it is n o longer based on G (global radiation) b u t on D (diffuse radiation). T h e useful e n e r g y c o l l e c t e d is G ..... in. T h e real energy Gu is calculated f r o m G ........ G~.mi~ and n/N, the ratio of the average n u m b e r of h o u r s o f bright sunshine to the m a x i m u m possible n u m b e r o f hours o f bright sunshine f o r the same p e r i o d assuming a c o m p l e t e absence o f c l o u d c o v e r (e.g., one m o n t h ) . T h e f o r m u l a is:
v
3
4000/1400
>. oJ
3800/1200
2 3600/1000
3400/800
3200/600
3000/400
2800/200
G,
=~Gu,ma
x +
1-
, Gu,mi n 0
This m e t h o d basically supposes t h a t the n hours o f bright sunshine are u n i f o r m l y dist r i b u t e d over the day which m e a n s t h a t the day is m a d e o f a l t e r n a t e sunshine periods lasting (n/N) dt and c l o u d y periods lasting (1 -n/N) dt. This a s s u m p t i o n is o n l y true if the time period is long e n o u g h (greater t h a n one week is a c c e p t a b l e b u t one m o n t h is preferable). T h e energy G~ is c o m p u t e d r e p e a t e d l y for several values o f the p a r a m e t e r A T which allows curves f o r G~ as a f u n c t i o n o f A T for each m o n t h to be d r a w n for each l o c a t i o n and each s y s t e m c o n f i g u r a t i o n (e.g. Fig. 9). Of course, these charts have been m a d e o n c e and f o r all, and have been published for a lot o f various climatic c o n d i t i o n s . T h e only thing the designer has to do is to m a k e use o f t h e m as described in the n e x t paragraph.
Solar project Let us consider t h r e e main f u n c t i o n s which are integrated in a solar h e a t e d h o m e n a m e l y solar e n e r g y collection, storage and heat distribution. T h e p a r a m e t e r which these functions have in c o m m o n is the t e m p e r a t u r e o f the water (or the air) flowing into the solar
10
20
30
40
50
60
70
(TF - Ta)(°C)
Fig. 9. Single p a n e c o l l e c t o r : m o n t h l y p e r f o r m a n c e profiles for z o n e E ( o r i e n t a t i o n 0 ° , i n c l i n a t i o n 45 ° ).
collectors when the h e a t e x c h a n g e r t e m p e r ature gradient in the storage t a n k is negligible. As a c o n s e q u e n c e , in a p r o j e c t , o n e has t o det e r m i n e u n d e r c o n s t a n t average c o n d i t i o n s the t e m p e r a t u r e o f the w a t e r (or air) and f r o m its k n o w l e d g e the m o n t h l y a m o u n t o f energy really saved due to the use o f solar heat. One can n o t e t h a t in this a p p r o a c h the size of the storage t a n k o n l y influences the transient b e h a v i o u r o f the system (transient means in a building, 1 - 10 h o u r s d u r a t i o n ) the practical sizes o f storage being of the order o f 4 - 5 m 3 o f w a t e r per h o u s e or apartment.
Average storage temperature in a given month. If T~., is the t e m p e r a t u r e o f the water flowing into the solar collectors and Ts the t e m p e r a t u r e o f the storage t a n k t h e n u n d e r c o n s t a n t c o n d i t i o n s TF = T~. T h e useful energy Gu p r e s e n t e d as a f u n c t i o n o f A T = (Tv -- Ta) can be evaluated as a f u n c t i o n o f
209
T~ because T, given by meteorological data can be averaged monthly for space heating evaluation in building. Let us draw for a given month on the same graph as a function of Ts the following characteristics: useful solar energy Eu (translated from Gu) - - s p a c e heating energy load, (kWh)Ec -
I
2000 = W
1600
\
--
Ta)
+
-
~ \ _E~_ _~
\
1400 Febru,~ry
EH) N
\
1200
where Go is the heat losses coefficient in W/m 3 °C V the living space volume in m 3 TNc the internal house temperature below which space heating is needed (depends on G and is equal to 15 °C for G = lW/ m 3 °C) Ta the external average air temperature in °C N the number of days of the month considered EH the domestic hot water energy demand (in France I0 kWh/day on average) solar energy distribution system fitted into the house:
E, ~J
1800
-
"t 2 4 G ° y (TNc E~ = \ iooo
I
2200
Eo\~i~i ,il [
~ 1000 April
,,,
\
800
60o
I: I
Ill',
4oo
t!',l l '
Ill
0
E~-A
,
•,
~
I
ii'
200
~
~-
"
~ ~
~
70
80
hot water dem~
,, i I
10
20TsTS°30 Tsl,40 April
50
60
90
Ts (°C)
Fig. 10. Calculation chart for February, April and August (Nimes house example).
-
heating
EDI = 24Nk I (T, -- Ti) kWh/month
water (at 60 °C)
ED2 ,=
k2N
Ts--Tw 50
where T i is the design temperature
kWh/ month
of the house
( 2 0 °C)
T~ kl k2
the cold water temperature (10 °C) the exchange coefficient (kW/°C) the maximal h o t water demand (10 kWh/day) The functions Eu = f(T~) and ED1 + ED2 have opposite slopes and will intersect at a point at which the useful solar energy (Eo) is greater or less than Ec (see Fig. 10). First case: Eo > Ec. The coordinates Eo and To of the above-defined intersection represent respectively the energy saved in heating demand o f the house for the month considered, and the average temperature of the water in the storage (T~m) and in the solar heat distribution system (Tdm). In this case To = Tam = T~m. Second case: Eo < Ec. The supplementary solar energy heats the storage tank up to an equilibrium average temperature equal to the intersection temperature T~I of the functions Eu = f ( T s ) and Ec = constant. The regulation
of the heat distribution system adapts the temperature of the water (either by an "on and o f f " or a progressive opening of a three way value )to the heating demand (Ts2). As a consequence the average water temperature Tdm of the heating system is different from the one of the storage and is equal to Ts 2. If Tsm = Tsl exceeds 60 °C, solar energy will be able to cover the whole h o t water demand of the house. On the contrary, if T~m = T~I is less than 60 °C, solar energy, though sufficient for the whole energy demand of the house (recall that we are studying the case E o > Ec ) will not be able to heat hot water up to 60 °C, and supplementary energy will be necessary. This supplementary energy is equal to the difference between the h o t water heating requirement (k2N) and the solar energy which heats the water from Tw (10 °C) to T~m (ordinate of the point of the water heating characteristics ED2 = f(Ts) for which the abscissa is Tsm). Examples
In this section, we are going to compare the t w o manual methods we have just been describing. For this purpose, the same project has been studied by the two methods one after another. The problem is to determine the
210
Solarradiation Outside temperature (°C)
Degreedays (basis 16,5°C)
R~f~rence Gj GT (90°,south) (90° , south) kWh/m2.day ~Wh/m2month
Fundamental
Total
AT
load
(°C)
(kWh/month)
Solar
parameters
f
X
Y
energy kWh/month
January
6.5
282
2,25
69,7
93,5
3028
2,93
0,5
0,28
848
February
9,3
192
2,8?
79
90,7
2100
3,71
0,81
0,47
987
March
11,0
170
3,14
97,3
89
1900
4,45
1,12
6,63
1197
April
12,8
122
3,15
94,5
87,2
1410
5,69
1,46
0,73
1028
May
15.2
70
2,75
85,3
84,8
892
9,03
2,08
0,83
740
June
23,4
3
2,58
77,4
76,6
210
33,5
8
1
210
July
23,4
0
2,85
88,3
76,6
186
39,1
10,3
1
186
August
22,1
0
3,18
98,6
77,9
186
39,8
11,4
1
186
September
17,5
20
3,38
101,4
82,5
382
19,8
6,78
1
381
October
13,5
110
3,06
94,9
86,5
1294
6,35
1,59
0,77
996
November
8,9
215
2,31
69,3
91,1
2347
3,57
0,64
0,35
821
December
6,8
274
1.94
60,1
93,2
2948
3,00
0,44
0,24
707
0,491
8287
Year
1458
1016
16883
j
Fig. 11. Typical project by the f-chart method.
Outside
Degree-
Temperature
days
(°C)
I(basis 16,5°C)
January
6,5
282
February
9,3
March
Specific useful
Distribution
Storage
Storage
temperature (Oc)
efficiency
solar
temperature (Oc)
2,791
0,980
25,5
25,5
192
2,145
1,260
26,5
11.0
170
1,751
1,290
April
12,8
122
1,342
May
15,2
70
June
23,4
July
energy (kWh/m 2. day)
Global useful
Solar
energy (kWh/month) load
solar
Contribution (%)
0,760
3028
808
26,7
26,5
0,753
2100
930
44,2
27,0
27,0
0,749
1900
1049
55,2
1,150
27,0
27,0
0,749
1410
905
64.2
0,822
0,709
24,5
26,5
0,753
892
579
64,9
3
0,200
0,200
21,5
56,0
0,538
210
113
53,8
23,4
0
0,171
0,171
--
64,6
0,476
186
89
47,6
August
22,1
0
0,171
0,171
-
81,0
0,356
186
66
35,6
Septembar
17,5
20
0,363
0,363
21,5
84,0
0,336
382
128
33,5
October
13,5
110
1,193
1,135
26,5
38,0
0,669
1294
824
63,6
November
8,9
215
2,235
1,030
25,5
25,5
0,760
2347
822
35,0
December
6,8
274
2,717
0,840
24,5
24,5
0,767
2948
699
23,7
16883
7012
41,5
Months
Year
load
1458
Fig. 12. Typical project by the EDF method.
211 solar c o n t r i b u t i o n o f the f o l l o w i n g s y s t e m o n a y e a r l y basis (space h e a t i n g + d o m e s t i c h o t water). Building: o n e - f a m i l y h o u s e , l o c a t e d near Avignon ( F r a n c e ) . Latitude: 4 3 ° N . Living area: 140 m 2 Loads: t h e r m a l loss c o e f f i c i e n t : G = 1.2 W/m 3 °C, degree-days: as listed in Figs. 11 and 12, d o m e s t i c h o t water: 150 1/day ( t e m p e r a t u r e s 14 - 50 °C) Collectors: area: 35 m2; s o u t h facing; tilt angle: 90 ° (vertical); d o u b l e glazing: F R U L = 4 W/m 2 °C, F ~ (sT) -- 0.65 Storage tank: c a p a c i t y : 3 0 0 0 1 (85.7 1/m 2) D i s t r i b u t i o n heating floor, emission coeffiequipment: c i e n t 300 W/°C Auxiliary electric h e a t i n g b y individual energy: convectors T h e results are s h o w n in Figs. 11 and 12.
CONCLUSIONS It was n o t possible in such a p a p e r to s h o w all t h e details o f the solar models, n o r even to describe t h e i r t h e o r e t i c a l and practical procem o n s t r a t e t h a t solar m o d e l s play a great rSle
dures in an extensive way. We p r e f e r r e d t o dein all solar research programs t h a n k s t o their flexibility, and t h a t u n d e r n o c i r c u m s t a n c e s m u s t c o m p r e h e n s i v e and simplified m o d e l s be c o m p a r e d . Each o f t h e m has its o w n scope of applications.
REFERENCES 1 G. P. Mitalas and D. G. Stephenson, Cooling load calculations by thermal response factor method, A S H R A E Trans., 73 (1967) part I, paper 2018. 2 R. S. Varga, Matrice iterative analysis, Series in Automatic Computation, Prentice-Hall, 1962. 3 H. C. Hottel and A. Whillier, Evaluation of flatplate collector performance, in Trans. of the Conference on the use of Solar Energy, Vol. II, Thermal Processes, Univ. of Arizona (1955). 4 E. Aranovitch, The joint solar collector testing programme of the European community, in Proc. ISES (UK Section) Meeting (April 1977). 5 D. J. Close, A design approach for solar processes, Solar energy, 11 (2)(1967). 6 G. Gutierrez et al., Simulation of forced circulation water heaters: effects of auxiliary energy supply, load type and storage capacity, Solar Energy, 15 (1974) 287 - 298. 7 W. A. Beckman, S. A. Klein and J. A. Duffle, Solar Heating Design by the f-Chart Method, Wiley, NewYork, 1977. 8 Ph. Chouard, H. Michel and M. F. Simon, Bilan thermique d'une maison solaire -- Mdthode de calcul rapide, Eyrolles, Paris, 1977.