Finite element thermal analysis of special building components

Finite element thermal analysis of special building components

ENE&G¥ AND ELSEVIER Energy and Buildings 22 (1995) 115-123 BUILDING Finite element thermal analysis of special building components L. Agnoletto, G...

658KB Sizes 2 Downloads 59 Views

ENE&G¥ AND

ELSEVIER

Energy and Buildings 22 (1995) 115-123

BUILDING

Finite element thermal analysis of special building components L. Agnoletto, G. Cortella, M. Manzan Universit~ degli Studi di Udine, Istituto di Fisica Tecnica e di Tecnologie lndustriali, Viale Ungheria 41, 33100 Udine, Italy Received 16 October 1993; revised 30 October 1994

Abstract

The finite element method, which is usually employed to study thermal conduction within solid components, has been utilized to solve coupled conduction, convection and long-wave radiation problems. As an example, the thermal behavior of special building envelope components under transient conditions is investigated, with the aim of establishing the rate of heat transfer due to the difference between the internal and external temperatures and to the absorption of the incoming solar radiation. The components considered in this paper consist of a two-pane glazing system, an open cavity with air circulation, and an opaque wall. Assuming a sinusoidal variation in the absorbed solar radiation, the temperature distribution and the heat fluxes are estimated. Comparisons with standard components are illustrated, and the efficiency of the system has been determined as a function of the climatic data. Keywords: Finite element method; Thermal analysis; Special building envelope components

I. I n t r o d u c t i o n

2. S y s t e m s i m u l a t i o n by the finite e l e m e n t m e t h o d

Interest in solar energy for the reduction of energy consumption for heating purposes began more than 20 years ago. Research, development and demonstration programs have been supported in many countries, and many theoretical and experimental studies have been carried out and are still proceeding [1--4]. A thorough knowledge of solar building components is needed to allow correct design and to take the best advantage of the available solar energy. In this paper we employ a new finite element approach based on the modeling of coupled conduction, convection and radiation heat transfer [5,6]. Using this method, most building components can be investigated accurately. As an example of its application, the paper deals with the evaluation of transient heat fluxes under winter conditions with reference to complex solar envelope components. Temperature fields in the building envelope are usually multidimensional and time dependent. Therefore, in general, they can be computed accurately by numerical procedures only. In the present study the common assumption of two-dimensional temperature fields was made. The resulting computer program has successfully been validated by comparison with reliable numerical and experimental data reported in the literature [7,8].

A schematic diagram of the building component investigated in this paper is presented in Fig. 1. As can be seen, we have a double-pane window (consisting of an external glazing pane, a closed cavity, and an internal glazing pane), a cavity connected with the inside and an opaque wall which can be insulated on the side facing the cavity. Inside and outside temperatures are assumed to be constant. The heat exchanges between the environment and the glass or wall surfaces are governed by combined convection and long-wave radiation. The heat transfer coefficients of convection and the linearized equivalent heat transfer coefficients of radiation are assumed to be constant. Two openings through the wall, at the lower and uigper levels, allow air circulation between the cavity and the room. A special device stops the air circulation when the temperature within the cavity is lower than the temperature in the room, thus preventing flow reversal. The heat transfer mechanisms are: - convective and radiative heat exchange between the outside and the double-pane glass; - conductive, convective and radiative heat exchange within the double-pane glass; - radiative heat exchange between adjacent walls in the cavity;

0378-7788/95/$09.50 © 1995 Elsevier Science S.A. All rights reserved SSDI 0378-7788(94)00908-3

116

L. Agnoletto et aL / Energy and Buildings 22 (1995) 115-123

work the fluid regions are assumed to always be narrow, and they can be schematized as channels, in which the fluid is flowing, or as enclosures, in which the flow rate and the average fluid velocity are equal to zero. In the solid region, the energy transport is governed by the heat conduction equation, which can be cast in the form

double pane window

V. (kVt) + g = O

air cavity

(1)

where ~ is a generalized source/sink term

thermal insulation

Ot

= 0 - pc O-O

opaque wall

(2)

expressing the difference, per unit volume, between the rates of energy generation and of energy storage in the material. Appropriate boundary conditions for Eq. (1) are t =tp

(3)

on the section of the boundary where the temperatures are prescribed, q"= - k V t . n

(4)

on the section of the boundary where the flux is prescribed, and q" = h c ( t - t f ) + h r ( t - t s ) = - k Vt.n Fig. 1. Schematic representation of the special building component investigated.

i b

_

/,4

m Fig. 2. Domain of definition for coupled conduction and convection problems with (a) solid regions; (b) channels where a fluid flows; (c) enclosures where the flow rate is zero.

- convective heat exchange between the cavity walls and the circulating air; absorption of the incoming solar radiation; - heat conduction through the opaque wall; - convective and radiative heat exchange between the wall and the room. In order to study all these phenomena we use the finite element technique [6]. A typical domain of definition for coupled conduction and convection problems is represented in Fig. 2, where a solid and a fluid region are interconnected in various ways [5]. In this

(5)

on the section of the boundary where heat transfer occurs by convection and radiation. Convection occurs with a fluid at a specified temperature h, while radiation occurs with surfaces at a specified temperature ts. In Eqs. (4) and (5), a positive value of q implies flow of heat out of the domain, in the direction of the outward normal n. In coupled convection and conduction problems we must also take into account the contributions from the fluid regions, i.e., from the channel and the enclosure within the solid. For surface convection in a narrow enclosure (d/L << 1), the energy exchanges at any point can be expressed as dq = (hc(ts,i - If) + hr(ts.i- ts0)) d F

(6)

Only the bulk fluid temperature variations are of interest for a fluid within the channel. Thus, for a duct of high aspect ratio (d/L << 1), the most general form of the energy equation is

,fc,S,

+ u,-d ) = g

(7)

In Eq. (7) the surface coefficients are first evaluated separately and then specified as input data, while the mass flow rate is evaluated taking into account the temperature variation along the channel. This results in a stack effect that is similar to an internal rate of pressure generation. On a channel F the pressure

L. Agnoletto et al. /Energy and Buildings 22 (1995) 115-123

117

difference between the inlet and the outlet becomes

Ap =g[3pf1(t(r ) -

tref) d F ' e

(8)

F

where r is a curvilinear coordinate, and £ is the unit vector parallel to the gravity field and pointing upward. This pressure difference induces a mass flow rate governed by the relationship

ap (9)

pfufSf- g

In the present problem, the fluid is flowing in a narrow channel where fluid convection, surface convection and surface radiation occur simultaneously. To take all these contributions into account, two kinds of elements are introduced: surface convection elements and fluid convection elements. Then, as shown in Fig. 3, we can utilize fluid convection elements, such as the element V, between surface convection elements, such as elements III and IV. In this way the surface convection

--7 / / / I

/

/

/

ds

/

" m

5,I

t

d

6

s=0

(a)

Is 1

2

4 4

1

3

4

1 5 6

21 4

-0--(

1'3

Uf

×

-1

2

- -

J

3

(b)

Fig. 3. Fluid and surface convection in narrow channels. (a) A fluid is flowing in a narrow channel where there is both fluid and surface convection; (b) the surface convection elements III and IV take into account the energy exchanges between the fluid and surfaces I and II, respectively, while the fluid convection element V takes into account the additional contributions to the energy balance due to the bulk motion of the fluid and element VI takes into account the long-wave radiative exchange between the opposed walls of the channel.

Fig. 4. Space-discretized mesh of the most general building component investigated•

elements take into account the energy exchanges between the solid surfaces and the fluid, while the fluid convection elements take into account the additional contributions to the energy balance due to the bulk motion of the fluid. A surface convection element, such as element VI, is added to take into account the longwave radiative heat exchange between the opposite walls of the channel. With this model, no special input data are required for the mesh definition, even if additional parameters, such as fluid velocities, convective and radiative heat transfer coefficients, are obviously needed for the characterization of the coupled problem• Following the standard finite element approach, the quantities resulting from the element energy balances are transferred to the nodes. In this way, after assembling the contributions of all the elements involved, the usual space-discretized forms are obtained [5]. The resulting space-discretized equations can be integrated in time by means of one of the many popular

L. Agnoletto et al. / Energy and Buildings 22 (1995) 115-123

118

50

14

(a)

(b)

45

12 40

6"

10

® 8

r

35

o

~ 25

16 E

~. eo ~15

4

/

10

2

5

0 [ 0

.

.

.

.

.

.

.

6

.

.

.

12

. 18

i

0

.

0

24

6

12

40 :

60 f

-

(c)

18

24

hour

hour

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

.

i (d)

55

35

50 30 ~

~

[

O" 25 ,

-.

~. 45

/.

~ 4o ~ as

e 15 -~ 25 10 " 20

°i 0

15 •

i 1

10 - 6

12

18

24

6

0

12

hour

18

24

hour

Fig. 5. Component (a): (a) temperature of the external glazing pane for A"= 500 W m 2; (b) temperature of the internal glazing pane for A"=500 W m-Z; (c) temperature of the air in the cavity for A"=500 W m-Z; (d) temperature of the external wall surface for A"=500 W m-2.

200

1.0 . . . .

/"~\

i!

0.9 160 !

0.8

~"\ " \

// /

0.7

120 ~

A"

~0.6

A"/

] \ \,,

/

~ 0.5

~ o.4

/

/

//

'~ 0.3

40 i

"\..

A'

// / /

0.1 0.0 6

./"

I

/,

0.2

o k, '.

/

.

. 12

.

7,, 18

Ay

_

........

J/

"J

i

-40 L

24

0

6

Fig. 6. Component (a): air velocity in the cavity for A'=200 and A"=500 W m -z. finite difference s c h e m e s available. In this paper w e favor the implicit m e t h o d because it is unconditionally stable, and does not lead to oscillations in the c o m p u t e d results.

3. E x a m p l e of an application The building e n v e l o p e c o m p o n e n t s c h e m a t i z e d in Fig. 1 was considered. D e p e n d i n g on the wall, the following cases w e r e investigated:

12

18

24

hour

hour Fig.

7. C o m p o n e n t

A"=500

W

(a):

h e a t flux by c o n d u c t i o n

for A ' = 2 0 0

and

m -2.

component (a): the opaque wall is c o m p o s e d of a single layer (non-insulated wall) of thickness s = 2 4 0 m m , thermal conductivity k = 1.2 W m -1 K -~, and density p = 1800 kg m - 3 ; component (b): the opaque wall (insulated wall) is c o m p o s e d o f two layers: an external insulating layer of thickness s = 40 m m , thermal conductivity k = 0.035 W m - 1 K - 1, and density p = 40 kg m 3; and an internal brick layer of thickness s = 200 m m , thermal conductivity k = l . 2 W m -1 K -1, and density p = 1 8 0 0 kg m -3.

L. Agnoletto et al. / Energy and Buildings 22 (1995) 115-123 1.0

The outside and inside temperatures are: outside to = 0 °C; inside t~= 20 °C. The temperature of the air leaving the lower level of the room and entering the open cavity is assumed to be t=18 °C. On the external surface of the wall, the absorbed solar radiation is assumed to change according to a sinusoidal law with a 24 h period:

0.9 0.8

119

A"

0.7

~0.6

~, o.s 0.4

"~ 0.3 0.2 0.1

q~=A sin[~2(0- 6)]

0.0 0

6

12

18

24

hour

Fig. 8. C o m p o n e n t (b): air velocity in the cavity for A ' = 2 0 0 A " = 5 0 0 W m -2.

and

200 160

g

120 8o 40 0 A' -40 6

12

18

24

hour

Fig. 9. C o m p o n e n t (b): h e a t flux by c o n d u c t i o n for A ' = 2 0 0 A " = 5 0 0 W m -2.

and

The glazing system consists of a double-pane window whose global conductance C is assumed to be C = 6 W m - Z K -1

(10)

The convective heat transfer coefficient to the air in the cavity, from both the glass and the wall, is given by the relationship hc=2.8+3v W m - 2 K -1

(12)

(11)

where v is the air velocity (m s-l). The long-wave radiative heat transfer coefficient, between the double-pane glass and the wall, is linearized and assumed to be constant. Its average value is hr = 5.5 W m - 2 K -1 The geometry of the channel in which the air flow takes place is characterized by the following parameters: width of the cavity d = 40 mm; height of the wall L = 3 m.

The convective and long-wave radiative heat exchanges between the window pane and the outside, and between the wall and the inside, are assumed to be characterized by the global average coefficients: outside ho=14.5 W m -2 K - l ; inside hi=7.0 W m -2 K -1.

Obviously, during the night period (6:00 p.m. to 6:00 a.m.) the absorbed solar radiation is set identically equal to zero. Two different peak values of the absorbed solar radiation have been considered: A = A ' = 2 0 0 W m -2 (average cloudy sky conditions), and A =A"=500 W m -2 (average clear sky conditions). Using the finite element code, the following parameters have been determined: the temperature distribution at each relevant surface; the air temperature and air velocity distribution in the cavity; the hourly values of the heat flux transmitted by conduction through the wall; the hourly values of the rate of energy supplied by the circulating air; the daily amount of total energy transferred to the interhal environment by the circulating air and by conduction through the wall; the fraction of the incoming solar radiation that can be effectively utilized. In order to compare the performance of the components investigated with that of standard building components, the daily energy gain was determined also for the following reference cases: component (c): the same component as (b), with the double-pane glazing glass replaced by an opaque shield and no air circulation between inside and outside. component (d): a multilayer insulated wall like that in component (b), but without any shielding or glazing element. In Fig. 4 we report the space-discretized mesh corresponding to the general building component represented in Fig. 1.

3.1. Non-insulated wall with external double-pane glass: component (a) In Fig. 5, for a peak value of the absorbed solar radiation A ' = 500 W m - z , w e report the hourly values of the temperature of the external and internal glazing panes, of the air in the cavity, and of the external surface of the wall, at the heights x/L = 0.2, 0.4, 0.6, 0.8 and 1.0. In Fig. 6 we report the hourly values of the air velocity in the open cavity for A ' = 200 and A"= 500 W m -2.

120

L. A g n o l e t t o

et aL / Energy and Buildings 22 (1995) 115-123

40

5O

(a)

45

(b)

35

40

//

30

35

0

25

--30

.Q

em 25

Z. 2o E ~ 15

[

10

10 5 ~ .

.

.

.

.

A xa-

J

12

18

0

24

6

12

30

90

(c)

80

18

24

18

24

hour

hour

70 6"60

, :

5O

F! ( d )

26 6" • ~ 22

R

E

[

20

18

10 14

0 6

12

18

0

24

6

12

hour

hour

Fig. 10. C o m p o n e n t (b): (a) temperature of the internal glazing pane for A " = 500 W m - Z ; (b) t e m p e r a t u r e o f the air in the cavity for A " = 500 W m - 2 ; (c) temperature of the insulating layer surface for A " = 500 W m - z ; (d) temperature of the insulating layer-wall interface for A " = 500 W m -2

450

900

400

800

350

l

A}/

'

700 I 600

! 300

'\

~" 2so ~2°° I 150 ~

50 ~ ,/

/

\

A'

/,/

//

100 ..j

\

~

i

/' ~ / /

,

, / / ~ / /'

\'. I

\

//

300

//

100 5

0

~ 400

/'

/\ \

A'

\

0 12

18

24

hour

__ ~\,,

0

6

12

'\ 18

24

hour

Fig. 11. C o m p o n e n t (a): h e a t flux due to air circulation for . 4 ' = 200 and A " = 5 0 0 W m -2.

Fig. 12. C o m p o n e n t (b): heat flux due to air circulation for A ' = 2 0 0 and A " = 5 0 0 W m -2.

In Fig. 7 we report the behavior of the conduction heat flux for both clear sky and cloudy weather conditions. For clear sky conditions (A"= 500 W m-2), air circulation occurs for more than 20 h per day, as shown in Fig. 6. The temperature of the external surface of the wall is almost always higher than the temperature of the internal surface. Consequently, the temperature of the air leaving the cavity is usually higher than that of the incoming air (assumed to be equal to 18 °C). We can also notice that, as shown in Fig. 7, the heat

flux by conduction through the wall is always entering the room. At each instant the solar radiation transmitted inwards is higher than the heat loss due to the temperature difference between the inside and the outside. For cloudy weather conditions (.4'= 200 W m-2), air circulation occurs only during a few sunny hours, as shown in Fig. 6. In this case, for many hours the temperature of the external surface of the wall is lower than the temperature of the internal surface, and consequently, as shown in Fig. 7, the conduction heat flux is directed outwards. In fact, the heat loss due to the

L. Agnoletto et al. / Energy and Buildings 22 (1995) 115-123

121

Table 1 Daily energy transferred through the building components investigated A (W m -2)

Daily energy (kJ) C o m p o n e n t (a)

100 150 200 250 300 350 400 500 600

C o m p o n e n t (b)

- 3 000 - 800 1 800 4 700 7 600 10 700 13 800 20 100 26 600

-600 1 800 4 600 7 400 10 500 13 600 16 500 23 100 29 900

Table 2 Overall transmittance of the wall as a function of air velocity in the cavity Air velocity (m s -1)

0.00 0.05 0.10 0.15 0.20 0.25 0.30

U (W m -2 K - ' ) C o m p o n e n t (a)

C o m p o n e n t (b)

1.41 1.54 1.65 1.74 1.81 1.85 1.89

0.55 0.57 0.58 0.59 0.60 0.61 0.61

temperature difference is higher than the heat gain due to the absorbed solar radiation.

3.2. Insulated wall with external double-pane glass: component (b) A thermal insulation layer is applied on the external surface of the wall in order to reduce the heat loss by conduction. As a consequence, an increase has to be expected in the rate of energy supplied by the circulating air.

C ompone nt (c)

Component (d)

- 2 900 - 2 900 - 2 800 - 2 700 - 2 700 - 2 600 - 2 600 - 2 500 - 2 400

- 3 000 - 2 800 - 2 600 - 2 400 - 2 300 - 2 100 - 1 900 - 1 500 - 1 200

In Fig. 8 the hourly values of the air velocity within the cavity are presented for both A' =200 W m -2 and A"=500 W m -2. By comparing Fig. 8 with Fig. 6, referred to the non-insulated wall, we can clearly detect shorter air circulation periods. In fact air circulation, with insulated walls, can be expected only in the periods of direct exposure to solar radiation. In Fig. 9 we report the hourly values of the heat flow by conduction through the wall for A' =200 W m -2 and A"=500 W m -2. For clear sky conditions (A"=500 W m-2), at any instant the heat flow by conduction is directed towards the inside, but its value is lower than that obtained with the non-insulated wall. For cloudy weather conditions (A'---200 W m-2), the heat flux by conduction is directed from inside to outside for several hours. The temperatures of the internal glazing pane, of the air in the cavity, of the external surface of the insulating material and of the interface between the insulating material and the wall are shown in Fig. 10 for A"=500 W m -2. As can be seen, during the sunny period both the temperature of the external surface of the insulating layer and the temperature difference across the insulating material are very high.

Table 3 Solar efficiency rt of the solar and standard components investigated A (W m -~)

100 150 200 250 300 350 400 500 600

Solar efficiency "O C o m p o n e n t (a)

C o m p o n e n t (b)

Component (c)

C omp o n en t (d)

0.52 0.53 0.56 0.58 0.60 0.62 0.64 0.67 0.68

0.30 0.38 0.45 0.50 0.54 0.57 0.59 0.63 0.66

0.000 0.000 0.003 0.005 0.007 0.008 0.008 0.009 0.010

0.000 0.003 0.013 0.020 0.024 0.027 0.031 0.032 0.034

122 0,7

L. Agnoletto et aL / Energy and Buildings 22 (1995) 115-123 .

.

.

.

.

.

.

.

.

In Eq. (13), the daily energy loss QL is defined as

.

QL = U( O~- 0o)(86.4)

5~r 0,S

I 0,4

15

q Ol3

C~

0,2 0,I 0,05

(14)

Iil

i 0,10

O,15

0,20

0,25

0,30

0,35

Q ~/Qs

Fig. 13. Solar efficiency of the building components investigated as a function of the climate parameter QtJQs: ( ~ ) component (a); ([151) component (b); (©) component (c); (©) component (d).

3.3. Transient thermal behavior of the different components

where U is the overall transmittance of the wall. In Table 2, the values of the overall transmittance coefficients of insulated and non-insulated solar components are reported as a function of air velocity within the cavity. The values of U at v = 0 represent the thermal transmittance between the room and the internal surface of the cavity. Thus, from the data given in Section 3, we have: U m i n = 1.41 W m -2 K -j for the non-insulated solar wall (a); Umi,=0.55 W m -2 K -~ for the insulated solar wall (b); and, similarly, we obtain Um~n=0.55 W m 2 K-1 for components (c) and (d), which present the same structure up to the cavity. The solar efficiency is therefore defined as

QL+Q A comparison of the transient thermal behavior of components (a) and (b) leads to interesting considerations. In both cases the effective heat capacity of the wall is significant, and consequently the energy associated with solar radiation is transmitted inside with some hours of delay. The mechanisms by which the heat is supplied to the inside are very different in components (a) and (b). In component (a) the heat is transferred mainly by conduction through the wall, while in component (b) the presence of the insulating layer increases the energy contribution from air circulation. The rates of energy supplied by the circulating air are represented in Figs. I1 and I2, referring to components (a) and (b), respectively. In Table 1 the values of the daily energy Q transmitted to the room through the building components (a) and (b) and through the reference components (c) and (d) are reported as a function of the absorbed solar radiation. As we can see, the daily energy transferred to the inside through the solar components (a) and (b) is positive up to very low values of A, typical of cloudy days. The daily energy transfer is higher when wall (b) is used. However, the highest contribution is due to the circulation of air and is limited to sunny hours. This occurrence may cause overheating in the room, particularly in presence of windows. Therefore, the solar energy collected by the circulating air can be advantageously supplied to rooms that are not adjacent to the solar component (for example, rooms with northern exposure). In order to characterize the daily energy gain of the components, the solar efficiency r/ can be introduced. The net daily energy transmitted through a solar component can be written as Q = r/Q~- QL

(13)

Os

(15)

where QL is the daily energy loss, Q is the daily energy transferred through the component, and Qs is the daily absorbed solar radiation. In Table 3 the solar efficiency -q for the various cases is reported as a function of A. As can be expected, the efficiency of the non-insulated solar wall is higher than that of the insulated solar wall. In components (c) and (d) the efficiency is equal to zero for the lowest values of A, for which Q = - Qz~. in Fig. 13 the relationship between the efficiency and the climate parameter (QL/Qs) is reported for each component. It can be noted that for the insulated wall (d) the efficiency is lower than 5%. This means that 95% of the absorbed solar radiation is rejected in the environment. The efficiency increases to 68% for the non-insulated solar wall and to 66% for the insulated solar wall. The results reported in Fig. 13 can help building designers to determine the net energy transmitted through standard and solar components as a function of simple parameters. This kind of investigation can easily be extended to long period analyses in order to also evaluate the monthly efficiencies. 4. Conclusions

Coupled conduction, convection and radiation have successfully been studied with the finite element technique utilizing a simple new description of the domain. This procedure has been applied to the investigation of transient multidimensional heat transfer in some special building components. Two walls with an open air cavity and two standard walls have been considered. Their behavior in terms of temperatures reached and heat transferred to the inside has been investigated. Comparison of the results obtained confirms the high

L. Agnoletto et al. / Energy and Buildings 22 (1995) 115-123

efficiency of some building components in collecting solar radiation. The present method can be extended to many other different structures, in order to evaluate their thermal performance and give useful design guidelines. 5. Nomenclature

A c C d g h k L n p q q" 0 Q r,z R s S t U v

peak value of absorbed solar radiation specific heat thermal conductance width acceleration due to gravity heat transfer coefficient thermal conductivity height of the wall outwards normal to the boundary surface perimeter of flow passage heat flux heat flow rate rate of internal heat generation energy radial and axial coordinates pressure resistance linear coordinate cross section temperature overall transmittance velocity

Greek letters fl coefficient of volumic expansion F boundary ~7 solar efficiency

0 p

123

time density

Subscripts c f i o p r s

convection fluid inside outside prescribed temperature radiation surface

References [1]

[2]

[3]

[4]

[5]

[61

[7]

[8]

F. Trombe and A. Le Phat Vinh, Etude sur le chauffage des habitations par utilisation du rayonnement solaire, Rev. Gen. Therm., 48 (1965) 223-241. H. Akbari and T.R. Borgers, Free convective laminar flow within the Trombe wall channel, Sol. Energy, 22 (2) (1979) 165-174. T.R. Borgers and H. Akbari, Free convective turbulent flow within the Trombe wall channel, Sol. Energy, 33 (3/4) (1984) 253-264. W. Smolec and A. Thomas, Theoretical and experimental investigations of heat transfer in a Trombe wall, Energy Convers. Management, 34 (5) (1993) 385--400. G. Comini, O. Saro and M. IVianzan, A physical approach to finite element modeling of coupled conduction and convection, Numer. Heat Transfer, Part B, 24 (1993) 243-261. R.W. Lewis, W.K. Sze and H.C. Huang, Some novel techniques for the finite element analysis of heat and mass transfer problems, Int. J. Numer. Methods Eng., 25 (1988) 611-624. A. Fanchiotti et al., Passive solar systems, Final report of Progetto Finalizzato Energetica (in Italian), C.N.R., Rome, Italy, Sept. 1983. Passive Solar Design Handbook, US Department of Energy, Washington, DC, Jan. 1980.