- Email: [email protected]
A two-port envelope model for building heat transfer
\ PERGAMON
Building and Environment 23 "0888# 08Ð29
A two!port envelope model for building heat transfer C[ Lombard \ E[ H[ Mathews Department of Mechanical Engineering\ University of Pretoria\ South Africa\ 9991 Received 08 October 0884 ^ revised 08 August 0885 ^ accepted 8 October 0886
Abstract A zone envelope thermal model is derived from _rst principles\ based on the exact two!port solution to the heat conduction equation[ A single set of matrices\ one for each frequency of interest\ is used to describe the thermal behaviour of the total envelope[ The set of matrices must be computed once for a particular zone and cooling loads for various conditions\ or the passive response\ can then be obtained[ The method is based on the same principles underlying the CTF method[ However\ the transfer function is calculated for the whole envelope by combining the two!port matrices of the individual walls\ after simpli_cation of the room interior heat transfer network[ The conduction transfer coe.cients of the individual walls forming the zone are not required[ Except for the simpli_cation of interior radiation\ the method is completely exact and easy to implement[ It is also numerically e.cient[ For load prediction\ an implementation in the frequency domain\ based on the prime factor FFT algorithm is convenient[ For time domain simulation with active control\ the envelope conduction coe.cients are easily obtained by transforming to the Z!domain[ This represents a modern solution of the envelope conduction problem\ which is convenient to implement in modern object oriented computer programming languages[ Þ 0887 Elsevier Science Ltd[ All rights reserved[
Nomenclature A area of wall ðm1Ł A\ B\ C\ D elements of the transmission matrix C zone e}ective heat storage capacitance ðJ:KŁ Ca capacitance of interior air mass ðJ:KŁ Cm lumped capacitance of interior furnishings ðJ:KŁ Ea stored thermal energy in interior air mass ðJŁ Em stored thermal energy in furnishings ðJŁ F state transition matrix G state input matrix h convection coe.cient ðW:m1 = KŁ l thickness of layer ðmŁ k thermal conductance ðW:m = KŁ q heat ~ow ðWŁ Qc interior convective load ðWŁ Qcr system load for speci_ed interior temperature Tir ðWŁ Qir radiative loads on interior surfaces ðWŁ Qr radiative load ðWŁ Ra interior surface _lm resistance ðK:WŁ Rf e}ective resistance of interior forcing temperature ðK:WŁ Ri envelope e}ective interior resistance ðK:WŁ Rm resistance between interior air and furnishings ðK:WŁ Corresponding author[ Tel[ ] 99!16!01!319!320 ^ fax ] 99!16!01! 251420 ^ e!mail ] sto}el[LombardÝeng[up[ac[za
Ro envelope e}ective exterior resistance ðK:WŁ Rv ventilation resistance ðK:WŁ s Laplace domain variable ðrad:sŁ Tm transmission matrix Tmw wall transmission matrix Ti zone e}ective interior temperature ðKŁ Tir speci_ed interior temperature for system load Qcr ðKŁ To bulk temperature of outside air ðKŁ Tsa sol!air temperature of exterior surface ðKŁ Tf e}ective interior forcing temperature ðKŁ U state input vector X state vector Y thermal admittance matrix y00\ y01\ y10\ y11 elements of the admittance matrix Zo characteristic thermal impedance of a layer ðm1K:WŁ[ Greek symbols a thermal di}usivity ðm1:sŁ v radian frequency ðrad:sŁ u z"v:1a#"0¦i#l argument of two!port matrix[ Subscripts e external wall surface conditions i interior wall surface conditions s surface conditions t refers to the combined substructures[
S9259Ð0212:88:,*see front matter Þ 0887 Elsevier Science Ltd[ All rights reserved PII ] S 9 2 5 9 Ð 0 2 1 2 " 8 6 # 9 9 9 5 4 Ð 5
19
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Superscripts 0\ 1 indicating di}erent walls
0[ Introduction A very important part of a thermal simulation package for buildings is the zone thermal model[ It is also a very problematic part because the heat ~ow in a building is a complicated process involving conduction through the walls\ interior convection and radiation\ external insu! lation and convection\ penetrating radiation through fen! estration\ heat from equipment as well as interior occupant behaviour[ The graphical depiction of the vari! ous heat transfer paths by Clarke as well as his excellent discussion of the various issues are most instructive ð6Ł[ Clarke treats the development of four generations of thermal models[ These methods di}er mostly in the way they model the conduction of heat through the envelope[ In this study an extension of the highly simpli_ed model of Mathews\ Richards and Lombard ð12Ł is presented which is based on a rigorous treatment of the conduction process[ The method is based on representing the total zone envelope with one distributed heat ~ow element[ A thermal model for building designers can attempt to treat the heat ~ows in a building from _rst principles\ in a fairly rigorous\ exact manner[ This requires detailed knowledge of the _nal construction and environment ^ and also a good deal of heat transfer knowledge from the user of the model[ Another approach is to assume that few details are available and to model only the most important aspects[ This latter approach is viable if the answers so obtained are su.ciently valid to _x major passive design features\ such as window size and orien! tation[ However\ for prediction of the dynamic behaviour of HVAC systems\ the dynamics of the zone thermal response are important[ The simplest methods solve the conduction problem by lumping all the thermal mass in a single node ð0\ 04\ 7\ 19\ 21\ 11Ł\ so that a single! or second!order heat ~ow network is the result[ Solution of these models only require the integration of one or two di}erent equations\ instead of solving the partial di}er! ential equation for conduction[ The more accurate methods avoid lumping and tries to solve the conduction equation[ This article shows how a highly simpli_ed method can be extended to become a rigorous method\ without sacri_cing the bene_ts of simplicity\ by hiding the details of the conduction process[ The model of Mathews\ Richards and Lombard ð12Ł was derived with one overriding objective ^ to use a single node for all of the thermal mass[ This objective was motivated by the need to provide designers of passive buildings with a simple tool\ easy to understand\ and useful for evaluating design choices at the very beginning of the design stage when no details are available[ The validity of this highly simpli_ed method was established
Fig[ 0[ The single node model avoids the conduction equation and attempts to represent the building zone with a single heat storage capacitance[
by comparing predictions with actual measurements in 21 buildings[ The model represents the envelope with the Ro\ C and Ri elements in the network of Fig[ 0[ The other elements in the _gure are Tsa the outside e}ective sol! air temperature ðKŁ\ Qir radiative loads distributed over interior surfaces ðWŁ\ Ra interior surfaces e}ective _lm resistance ðK:WŁ\ Qc interior convective loads ðWŁ\ Rv ventilation resistance ðK:WŁ and To temperature of ven! tilating air ðKŁ[ This simple model is adequate for passive design pur! poses but it is not su.ciently detailed for simulation of the energy requirements of buildings with active climate control[ In an actively controlled building the system input may vary rapidly\ and only the innermost layer of the mass of the envelope often participates in the heat exchange[ Also\ HVAC system designers require detailed knowledge of the contribution of each surface to the total load and about the distribution of radiant energy across the interior surfaces[ For these purposes the single node model\ with no distinction between surfaces and all the mass lumped together\ is clearly inadequate[ 1[ Thermal models for buildings The prediction of the behaviour of building!environ! ment systems is discussed from _rst principles by Lord and Wilson ð10Ł[ The underlying principles are those of mass and energy conservation\ together with rate equa! tions from transport theory[ Building analysis models are necessarily simpli_cations of these fundamental prin! ciples[ Widely di}ering implementation and numerical strategies are possible[ A plethora of models and software tools exist[ Clarke discusses the state of the art in 0878 ð6Ł[ For a comparison of some North American work refer to Winkelmann ð24Ł and Tuddenham ð22Ł[ See also ð4\ 25\ 1Ł and the user|s perspective of Newton\ James and Bartholomew ð14Ł[ All these methods are attempts to solve the following two fundamental problems of thermal analysis of buildings ] 0[ The heat transmission through the building envelope must be calculated[ This requires the solution of the conduction equation\ which is a second order partial
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
10
di}erential equation in space and time coordinates ^ with appropriate boundary conditions[ 1[ The interior heat balance must be found[ The heat transfer between the interior air and interior surfaces\ the heat carried by ventilating air and air!conditioned air\ also by convective interior loads\ radiative interior loads as well as radiation penetrating through windows\ must be in balance with the energy stored in the bulk of the interior air[ The various methods di}er mostly with respect to the _rst of the two problems namely the solution of the heat conduction equation[ For the second problem only two approaches are clearly distinguishable ] the _rst is the heat balance method which solves the net of interior heat ~ows more or less exactly[ The other is to simplify the network by introducing the concept of an e}ective radi! ant or e}ective room temperature ð8\ 01\ 00\ 09Ł[ The _rst problem\ solving the conduction equation\ can be attacked by all of the very large body of analytical\ graphical and numerical methods available for solving partial di}erential equations[ Consequently\ the various thermal models for buildings are largely distinguishable by the methods and simpli_cations they employ for solv! ing the heat conduction equation[ It is probably fair to say that all the general thermal analysis programs for buildings simplify the conduction problem by assuming one!dimensional heat ~ow through walls[ This considerably reduces the problem as it allows surface temperature distributions to be approximated by a single representative temperature[ The assumption obviously fails for heat ~ow underneath buildings\ envel! ope corners and edges\ complicated envelope geometries and envelope structures with thermal bridges and special programs are available for dealing with these particular problems[ Nevertheless\ the simpli_cation to one!dimen! sional heat ~ow has been well justi_ed by experience over the past few decades[ Transverse heat transfer\ e[g[ from one part of a wall to another part of the same wall\ apparently\ plays a very small part[ If this possibility may be ignored\ it is always possible to divide a surface into smaller independent areas[ For instance\ a part of a wall which receives solar radiation through a window can be regarded as distinct from the rest of the wall which receives no radiation[ A feature of the one!dimensional conduction problem\ when applied to building walls\ is the existence of a num! ber of layers of material with di}erent thermal charac! teristics[ For strictly numerical methods this presents no complication at all ^ for analytical methods it implies the solution of the partial di}erential equation for every layer and then matching the boundary conditions[ Many approximate methods for solving the one!dimensional conduction equation change the partial di}erential equa! tion into a set of ordinary di}erential equations\ by lump! ing the distributed capacitance and resistance of the layer
Fig[ 1[ The capacitance of a massive element can be lumped to obtain an RC!ladder representation[
together at nodes\ as in Fig[ 1[ With these methods the degree of approximation is related to the number of nodes*which may be _xed or variable[ In a lumped method the temperatures between the nodes are unde! _ned[ The wall is represented with a discrete set of elements and a _nite number of temperatures[ It can be shown that most of the popular thermal models\ analytical and numerical\ make assumptions which reduce to lumping[ This is true of simpli_ed methods which explicitly lump\ and also of methods which use conduction transfer factors e[g[ the well!known ASHRAE!CTF method ð2Ł[ To construct the latter method the conduction equation was solved analytically\ and from this solution was determined a digital approxi! mation\ which is tabulated for use by designers ð08\ 20\ 02\ 03Ł[ The distributed element ideally requires an in_nite number of transfer factors[ In practice\ a _nite set is used which approximates the response of the in_nite set\ in similar manner to the approximation of the distributed wall element by a _nite number of RC!sections[ A popu! lar numerical method is to represent the partial deriva! tives of the conduction equation with _nite di}erences\ a process synonymous with lumping[ A numerical method which avoids lumping is the _nite element method[ With this method the temperature distributions of small elements are estimated by curve _tting while also ensuring continuity at the boundaries of the elements[ The _nite element method is very often used for three!dimensional heat ~ow in structures ^ but\ for computational reasons and because it is di.cult to use\ it is not favoured in building heat transfer[ For a general comparison of these methods refer to ð13Ł[ The one!dimensional heat conduction equation is an equation with space x and time t as independent vari! ables[ Lumping is a method for simplifying the spatial dimensions by concentrating the heat stored in the struc! ture at a few points in space[ A di}erent approach is to simplify the time coordinate instead[ The conduction equation is a linear equation\ hence complicated solutions can be constructed by superimposing simpler solutions[ Fourier methods are based on Fourier|s theorem that all
11
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
periodic solutions can be constructed from sinusoids[ Fourier|s theorem allows the construction of arbitrary periodic responses by adding the responses of the Fourier components[ Since no simpli_cation of the spatial coor! dinate is assumed the procedure is exact for periodic functions[ The method calculates the response under the assumption that the time coordinate is periodic[ This assumption is well suited to the quasi!periodic time dependencies experienced in practice[ Non!periodic responses can be calculated to any degree of accuracy by assuming a su.ciently long period[ However\ for very long periods\ and for short sudden chances\ it is less attractive from a computational point of view ^ requiring a large number of frequencies to be taken into account[ The suitability of this frequency domain approach has been most vigorously argued by Haghighat and Athi! enithis ð05Ł and long before them by Rao and Chandra ð16Ł[
Fig[ 3[ The two!port representation of a sandwich structure is simply the product of the two!port matrices for each layer[
2[ The matrix method The matrix method uses the solution for the con! duction equation with sinusoidal forcing functions[ The construction of the solution is well described in the litera! ture ð07\ p[ 49Ł\ ð08\ 20Ł[ It represents conduction through a layer with a 1×1 matrix as in Fig[ 2[ Mathematically the solution to the conduction equation is written in the form of a matrix operator which transform conditions at one surface into the corresponding conditions at the other surface[ With T0 the temperature and q0 the heat ~ow at surface 0 and T1\ q1 the corresponding quantities on the other the two!port model are T1
$ % $ q1
A B C
D
T0
%$ % ×
q0
"0#
where the two!port cascade "transmission# matrix is A cosh u B Zo sinh u C
0 sinh u D cosh u[ Zo
Fig[ 2[ Representing a massive layer with a two!port matrix[
"1#
In this equation u z"v:1a#"0¦i#l and Zo l:Aku are parameters of the layer which depend on the frequency v and physical properties ] A is the wall area\ l is the layer thickness and k is the thermal conductivity of the layer[ The matrix solution is very convenient for a layered structure[ The matrix for a combination of layers is sim! ply the product of the matrices for the individual layers as in Fig[ 3[ The matrix solution to the conduction equation characterises the heat transfer with four complex elements[ Of this four only three are independent and the fourth element is _xed by the condition that the determinant of the matrix be unity[ The CTF method\ as described in ASHRAE Fundamentals ð2Ł\ uses two of the elements of the matrix to calculate conduction from the external surface to the interior of the wall and to the wall interior surface[ To model the ~ux of heat from the interior surface to the interior mass of the wall the other two matrix elements are required[ Here the CTF method is inconsistent ] it combines the walls and uses a room transfer function "RTF# to take into account the storage of penetrating radiation in the walls\ while determining the sensible heat gain by conduction through the walls individually for each wall[ Besides being convenient for solving heat ~ow through a layered structure\ the matrices can also be used for adding the contributions from di}erent walls in a con! sistent manner[ If the model is consistent it will allow the decomposition of the total response into the con! tributions of the individual walls[ With this method the matrix for the whole zone envelope is found once and the matrix is then used as a complete representation of the envelope for all further calculations[ In the next section it is shown that a single 1×1 frequency dependent matrix characterises the whole building envelope and that it pro!
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
vides an exact solution for the heat ~ow through the envelope[ Most of the more accurate building thermal simulation programs are based on the matrix method[ The ASH! RAE!CTF method and programs like DOE1e\ BLAST\ HVACSIM¦ use the matrix method to calculate the heat transfer through a wall and from this exact solution an approximate di}erence equation is extracted[ This di}erence equation can be determined by approximating either the frequency response or the impulse response of the matrices ð20Ł[ The di}erence equation steps from the frequency domain back into the time domain[ For a speci_ed wall the CTF|s are determined beforehand and can then be used to determine the response for any\ known beforehand or not\ forcing functions[ ASHRAE Fundamentals tabulates the sets of CTF values of bn and dn for various walls and the user may look them up and apply them directly[ The CTF method treats each wall on its own and combines the heat!~ows through the walls with other loads to _nd a heat load to the interior air\ the conduction heat gain "CHG#[ The cooling includes the conduction heat paths to the interior air\ direct con! vective loads and also radiation penetrating through the windows and striking the interior surfaces[ This heat is initially absorbed by the structure and slowly released to the interior air[ The CTF procedure takes this into account by the use of another di}erence equation\ the room transfer function "RTF# ð17Ł[ The room transfer function represents the combined e}ect of all the walls of the zone[ The sets of RTF values of v and w are also tabulated in ASHRAE for a small number of rep! resentative zones[ An argument against the Fourier series method is that it cannot be applied to systems simulation[ In simulation a controller intervenes to regulate the system[ The inter! vention is based on the instantaneous interior tem! perature so that the response of the system is not known beforehand[ For Fourier analysis\ the forcing functions must be known beforehand to be able to calculate the spectra[0 But\ for passive analysis\ and to _nd the required system response that will maintain a speci_ed interior temperature\ the frequency domain method is accurate and easy to implement[ To simulate ~oating interior tem! peratures it is better to transform to the time domain and to characterise the envelope with a di}erence equation[ The ASHRAE!CTF method\ as presented in the Fun! damentals ð2\ Chap[ 15Ł\ is also limited in this regard[ The tabulated coe.cients do not include the full set required 0 This argument is not entirely valid[ For some control laws it is possible to solve the set of equations in the frequency domain and to determine both the response of the building and that of the system simultaneously[ This is sometimes done by control engineers to _nd a control law which optimises the system response ^ e[g[ minimum energy controller[ However\ solving the equations for actual systems is not often feasible because of the strong non!linearities in practical control laws[
12
Fig[ 4[ Equivalent building zone thermal model with walls represented by a two!port matrix[
for a time dependent interior temperature\ only the summed c series "Sci# is given[ This means the interior temperature must be constant[ The simulation program HVACSIM¦ calculates the full set of CTF values via the Laplace transform method[ This means the solution for each wall is obtained in the frequency domain and then transformed to the time domain\ for solving the interior energy balance[ Davies ð03Ł also presents a purely time domain method for obtaining the full set[
2[0[ A two!port zone model The building thermal models used by HVACSIM¦ and in the ASHRAEÐCTF method are complicated networks of energy ~ow through the walls and between the interior surfaces of the walls[ It is much easier to explain building heat transfer in terms of a single envelope through which heat enters or leaves the zone[ This approach emphasises the importance of designing a building shell with the correct overall thermal properties[ To obtain this simple representation of building heat ~ow it is necessary to _nd a single element description of the envelope and to combine the forcing functions into external and internal e}ective functions[ This can be accomplished without sacri_cing any accuracy in the manner explained below[ Figure 4 is an extension of Fig[ 0[ In the _gure a single two!port matrix characterises the envelope[ The forcing functions on the exterior surfaces of the walls are com! bined in a single e}ective forcing function Tsa[ The matrix includes both exterior and interior _lm coe.cients[ Ti represents the e}ective interior temperature of the zone with a single node[ The other elements are ] Ca which represents the lumped mass of the interior air\ Cm the lumped mass of furniture and other interior structures\ Rm the thermal resistance between the interior e}ective temperature and the interior mass\ Tf and Rf represent direct heat ~ow to the interior air caused by ventilation\ air conditioning and convective heating as explained in Fig[ 5 Tf
RvRsQc¦RsTo¦RvTs Rs¦Rv
"2#
13
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Fig[ 5[ De_nition of direct heat gain to interior air source and resistance[
Rf
RvRs [ Rv¦Rs
"3#
Qr represents an e}ective radiative heat ~ow to the interior surface of the envelope[ Ra is the e}ective _lm resistance between the interior envelope surface and the e}ective interior node[ The de_nitions of Ra and Ti are based on simpli_cation of the interior heat transfer[ Dav! ies ð09Ł has published an extensive analysis of the issues involved[ He presents a model which characterises the room interior with two temperatures\ the rad!air model\ which is logically founded and allows the calculation of interior comfort temperature[ Furthermore\ it enables additional simpli_cation by allowing all the interior radi! ation sources to be connected only to the two interior nodes[ For the purposes of this paper it is su.cient to assume that it is possible to use this representation with some degree of accuracy[ In our implementation we have used the approximation Ra 0:"7A#\ where A is the area of the wall\ and we assume that Ti will closely correspond with the interior air temperature[ The inclusion of a more re_ned model awaits further investigation[ In Fig[ 4 the resistance Ra is shown absorbed in the two!port for the envelope[ The model in Fig[ 4 is a direct extension of the single node model of Mathews\ Richards and Lombard ð12Ł as in Fig[ 0[ It represents an attempt to retain the simplicity and ease of understanding of the single node model but to extend its accuracy and applicability by re_ning the model of the envelope and the inclusion of an interior air mass node and a separate node for interior massive structures[ Although it is convenient to have a single matrix rep! resentation for a zone envelope it is not a requirement[ It is possible to retain the matrix description for individual surfaces and with these matrices construct a heat ~ow graph which can be solved with standard network methods[ This is the approach described by Athienitis ð3Ł[ His method is a full frequency domain method which represents the zone with a network of admittances[ The solution procedure involves _nding an admittance matrix which combines the heat ~ows from the various walls with other loads to determine the _nal response[ This admittance matrix characterises the whole zone and not just the envelope[ In Figure 4 the simpli_cations are ]
0[ Massive elements in the interior of the zone are rep! resented with a lumped capacitance and an air contact resistance[ This approximation is not essential[ Anal! ogous to the two!port representation of the envelope a one!port representation of interior mass is possible[ 1[ The radiant energy exchange between interior surfaces is simpli_ed[ A single node is used which is involved in all energy exchanges\ both convective and radiative[ This simpli_cation is often employed since it allows the simpli_cation of the network of interior heat ~ows as in Fig[ 6 to a single node[ For a full discussion of its bene_ts and problems see Davies ð01\ 09Ł[ The same approximation is used to DOE1e ð29Ł to calculate RTF coe.cients and in the ASHRAE!CTF method[ In this model no lumping of envelope elements takes place[ The envelope is represented by the two!port matrix[ The outside forcing function is also completely exact[ The envelope can always be decomposed again into individual walls and the contributions of the walls will correctly sum to the total response[ However\ in the numerical procedures for solving the model of Fig[ 4 inaccuracies may arise[ These inaccuracies have no direct relation with the heat ~ow model[ For calculating the passive response to given outside con! ditions and speci_ed internal loads\ and for _nding the required cooling load for a speci_ed interior temperature\ a full frequency domain method is convenient[ The full frequency domain method _nds a Fourier series rep! resentation of the forcing functions\ calculates a matrix for each frequency and combine the responses for each frequency to obtain the _nal response[ For simulation and intelligent control a time domain approach is more convenient[ This requires transforming the matrix to a set of di}erence equations which can be used in a similar way as in the CTF method[ However\ in contrast with the CTF method\ a consistent approach is followed[ The various elements of the same matrix are used both for calculating the conduction through the walls and the response to interior radiative loads[ In mathematical terms the di}erence equations represent a transform of the matrix operator\ which represents the envelope\ from the Fourier domain to the Z!domain\ and all the tech! niques for accomplishing this are applicable[
Fig[ 6[ Transforming the interior heat ~ow network to the star rep! resentation yields a single node involved in all heat transfers[
14
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
2[1[ The envelope matrix
Ti qiRa¦Ts
In Fig[ 4 the envelope is represented with a single matrix[ It is convenient if this matrix includes both sur! face _lm resistances[ The radiative sources\ which work directly on the surfaces and not through the _lm resist! ances\ are then connected to points which are interior to the matrix[ Also\ the forcing functions on the surfaces of the walls are di}erent for each surface[ The single matrix representation requires that the forcing functions be reduced to two e}ective forcing functions ^ one for the outside and another for the inside surface of the envelope[ This reduction of the number of forcing functions is a crucial aspect of the simple representation[ On the outside surface the forcing functions can be combined in a single temperature forcing function by using the well established Sol!Air concept[ The external solar!\ long wave!radiation and convection with outside air is combined by _nding the Thevenin equivalent circuit of these forcing functions[ The Thevenin source is then the Sol!Air temperature and the Thevenin resistance the e}ective _lm resistance which can be included in the matrix[ For the interior surfaces\ which are not termin! ated\ but will be connected to other surfaces and the interior air node it is not so simple[ In Fig[ 7 the model for a single wall is given[ In this _gure the exterior e}ective "Thevenin# temperature is rep! resented by Te and the exterior surface _lm resistance has been absorbed into the two!port for the wall\
and
Tmw
$
Aw Bw Cw
%
Dw
Ts
$% $ qs
Aw Bw Cw
Ti
×
"qs−Qr#Ra¦Ts
$% $ qi
Te qe
0 Ra
9
×
0
Ts
Ra
−Qr
qs
0
[
"7#
But with eqn "4# conditions on the interior surface can be written in terms of exterior conditions ] Ti
0 Ra
$% $ %$ qi
9
×
0
Aw B w
Te
Ra
%$%
%$Cw
×
Dw
qe
−Qr
"8#
0
showing that despite the presence of the source Qr\ the matrix for the wall and interior _lm coe.cient is still the product of the two individual matrices and both _lm resistances can be included in the de_nition of the wall! matrix Tm
$
A B C
0 Ra
% $ %$
D
9
0
×
Aw Bw Cw
Dw
%
"09#
[
The _nal solution for the wall is
A B C
D
Te
Ra
%$% $% ×
qe
−Qr
0
"00#
[
To determine the heat ~ows through the surfaces of the wall eqn "00# must be solved for qe and qi in terms of the surface temperatures[ This gives qe
−ATe¦Ti¦RaQr 0 B "BC−AD#Te¦DTi¦"DRa−B#Qr
$% $ qi
%
"01#
which can conveniently be written qe
"4#
% $ %$%
%$qs−Qr
¦
qi
%
%$Dw
Hence
$% $
The heat ~ow across the exterior surface is qe[ Ts is the temperature on the inside surface of the wall and qs the heat ~ow across the inside surface[ Ra is the _lm resistance between the interior surface and the interior air\ which is at temperature Ti[ The heat ~ow to the interior surface from the interior air is qi[ Qr is radiation heat ~ow to the interior surface[ The matrix solution for the interior surface of the wall gives ]
"6#
qi qs−Qr[
Ti
[
"5#
−A 0 0 B BC−AD D
$% $ qi
%$
Te Ti¦QrRa
9
% $ % ¦Qr
−0
[ "02#
and from the network
In eqn "02# we see that the e}ective interior wall surface forcing temperature is Ti¦QrRa[ However\ it is more convenient to keep the forcing functions separate and to write qe
−A 0 0 B BC−AD D
$% $ qi
Fig[ 7[ The thermal model for a wall consists of a two!port matrix and exterior and interior forcing functions as de_ned in the text[
Te
Ra Qr [ B DRa−B
%$ % $ Ti
¦
%
"03#
In eqn "03# the two!port matrix for the wall is in the admittance or Y format[ If the B element of the Tm matrix is non!zero it can always be transformed to the
15
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Y
$
y00 y01
qet
%
qit
format with y00 −
A 0 y01 B B "04#
y00 0 B y01 y01 y00y11 y11 D [ y01 y01
"05#
Qre
Ra y10Ra Qr Qr B DRa−B y11Ra−0
$ % $ Qri
% $
%
"06#
eqn "03# can be written compactly "07#
q YT¦Qr with the further de_nitions q
qe
$% qi
qi1
"19#
and T
Te
$% Ti
[
"10#
where subscript t is used to indicate variables which are representative of the combined structure[ This important equation gives the contribution of the various forcing functions on the individual walls to the total heat ~ow across the envelope[ To _nd a single two!port representation for the envel! ope eqn "10# must agree with qt YtTt¦Qrt[
"11#
This implies that
The T! and Y!representation of the two!port are math! ematically equivalent[ The T!representation is convenient when systems are in cascade while the Y!representation is used for systems in parallel where the ~ows add[ If a radiative source vector is de_ned by Qr
qe1
qt Y0T0¦Qr0¦Y1T1¦Qr1
If y01 is non!zero the inverse transform is
C y10−
qi0
¦
or
AD D y11 [ y10 C− B B
A −
qe0
$% $ %$ %
y10 y11
"08#
The building envelope consists of a number of surfaces which enclose the interior space[ The envelope heat ~ow is the sum of the heat ~ows through all the surfaces[ In Fig[ 8 the envelope\ for simplicity\ has only two walls[ The total heat ~ow through the envelope is
Yt Y0¦Y1
"12#
Yt Y−0 t "Y0T0¦Y1T1#
"13#
Qrt Qr0¦Qr1
"14#
which de_nes the e}ective envelope Y!matrix and forcing functions[ The equation for the combined heat ~ow of two walls is in exactly the same form as that for a single wall[ This makes it possible to extent the model from two walls to any number of walls by accumulating the walls one by one into a total representation for the whole envelope[ The Y!matrices of the walls are added together and the e}ective forcing functions are calculated from a weighted average where the weights are also obtained from the elements of the Y!matrices[ 2[2[ Solution With the envelope matrix known the circuit of Fig[ 4 can be solved with the envelope regarded as a black box described by eqn "11#[ For passive analysis the interior temperature is the unknown variable[ For calculating the HVAC system load the heat ~ow to the interior air for a speci_ed interior air temperature is the unknown[ To solve the network of Fig[ 4 standard nodal analysis ð23Ł is used[ For state variables we choose the thermal energy in the interior air mass Ea TiCa and the thermal energy in the massive interior structures Em TmCm[ At node i the heat ~ows must add to zero\ hence qi¦qa¦qm¦qf 9[
"15#
But qa
Fig[ 8[ Thermal model for an envelope consisting of only two surfaces[
d Ea dt
"16#
d Ti−Tm Em dt Rm
"17#
qm
16
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
qf
Ti−Tf Rf
"18#
and qi is supplied by eqn "11# "29#
qi y10Tsa¦y11Ti¦Qri
where the variables now refer to the totals\ summed for all walls forming the envelope[ It is convenient to de_ne qr y10Tsa¦Qri qi−y10Ti\
"20#
a quantity which is known as it involves only the known forcing functions[ If these are substituted in eqn "15# we have
0
qr¦y11
1
0
1
0 0 d d 0 E a ¦ Ea ¦ Em ¦ Ea− Tf 9[ Ca dt dt CaRf Rf "21#
At node m the heat balance is ] qm
d Ti−Tm 0 0 Ea− Em[ Em dt Rm CaRm CmRm
G
0
0
Ea
0
1
1
0 0 Em− Tf 9[ CmRf Rf
X
$ % Em
and U
$
Tf:Rf−qr 9
%
"39#
"25#
"30#
Tf−To −Qc Rf
Qcr
"31#
where "32#
and qf −"qi¦qa¦qm# qi qr¦y11Tir
d X FX¦GU dt
"26#
where
1
0 0 K −0 H Ca y11¦ Rm ¦ Rf F H H 0 k CaRm and
1
[
are required from the envelope so that only two matrix elements are calculated and stored[ In the above equa! tions Tf\ Ea\ qr\ Qa\ Tsa\ y11 and y10 are complex numbers giving the amplitude and phase at frequency v[ It is also possible to _nd the required extra convective load Qcr which will ensure a given interior temperature Tir[
Tf Tir−Rfqf
the state equation in standard vector form is
0
0
0 Rf 0¦Rmy11¦sRmCa− ¦Rm sRmCm¦0
"24#
Equations "24# and "22# are the state equations for the circuit[ If the state vector X\ and the source vector U\ are de_ned by Ea
CaRm"Tf−Rfqr#
qr y10 ( Tsa¦"y11Ra−0#Qa and y11 "23#
0 0 0 d y11¦ ¦ Ea¦ Ea Ca Rm Rf dt −
"28#
[
In this solution only
This equation simpli_es to q r¦
9 0
"22#
1
0 0 E − T 9[ CaRf a Rf f
¦
$ %
This state eqn "26# is the complete description of the zone[ It|s solution\ for given forcing functions\ gives the response of the zone to the forcing functions in terms of the energy stored in the air!mass and interior mass[ The format of eqns "25# and "26# is misleading since it indi! cates a second order system[ In fact\ the envelope oper! ators y10 and y11 are high order di}erential operators so that the actual order is much higher[ Equation "26# re~ects the simplicity of Fig[ 4[ Note that energy stored in the structure is not visible in the equation because of the black box hiding the complexity of the conduction process\ which is treated in a rigorously exact manner[ To calculate the passive temperature response of the zone eqn "26# must be integrated to _nd Ea and Em[ The interior air!temperature is then Ti Ea:Ca[ In the frequency domain this is straightforward[ Writing s jv for the frequency dependence and solving eqn "26#
Substituting this result in eqn "21# gives 0 y11 d 0 qr¦ Ea¦ Ea¦ Ea− E Ca dt CaRm CmRm m
0 9
0 L CmRm H H −0 H l CmRm
"27#
qa sCaTir qm
Tir−Tm Rm
Tm
Tir 0¦sRmCm
"33#
with qr as de_ned above[ The frequency domain solution involves ] calculating and storing the two envelope matrix operators y10"v# and
17
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
y11"v#\ _nding the frequency components of the forcing functions and calculating the frequency components of the solution through 39 or 31[ For a time domain implementation the y10 and y11 operators must be trans! lated to di}erence operators[ Various standard tech! niques exist as discussed in ð20\ 08\ 17\ 18\ 06Ł[ This leads to a method similar to the CTF!method\ except that the CTF!coe.cients are calculated for the whole envelope from the accurate matrix description of the envelope\ instead of for each surface[ 2[3[ Implementation For implementation of the scheme a modern object oriented approach is ideal[ In an OOP language such as C¦¦ a software object is created for each physical object[ Software objects are also de_ned to store and manipulate the data structures[ We have implemented an object hier! archy consisting of the following data objects ] 0[ A two!port object which allows matrix algebra and T \ Y transforms through operator overloading[ 1[ A time series object which can retrieve\ store and print a time sequence\ and transform it to the frequency domain using a mixed radix FFT algorithm suitable for sequences with length a multiple of 13 ð5Ł[ 2[ A spectral object which can store and print a sequence of frequency components and transform them to a time series representation with the same FFT algo! rithm[ The software representation of the zone consists of the following {physical| objects ] 0[ A layer object which contains material and dimen! sional data for a wall layer[ 1[ A surface object which contains and combines a sand! wich structure of layers with an external Sol!Air tem! perature series and an interior radiative load series[ 2[ An envelope object which combines a number of sur! faces in a zone envelope and _nds the equivalent envel! ope description[ 3[ A zone object which combines the envelope with the interior air\ interior mass\ ventilation and convective forcing functions in the manner indicated in Fig[ 4[ To verify the code exercise routines were written for each object and the answers were thoroughly checked by hand calculations[ All the calculations were veri_ed in this way by comparing the mean response and the swing for a single sinusoid forcing function[ The surface object was veri_ed by calculating the heat ~ow for the wall object of example 2\ Chap[ 15\ in ð2Ł[ In this example the cooling load is calculated to maintain a temperature of 13>C behind a wall constructed of 099 mm heavy concrete\ with 49 mm insulation\ and 19 mm indoor plaster ^ for a prescribed external Sol!Air temperature[ The same struc! ture was also calculated with the CTF method\ with
Fig[ 09[ Results for the wall of ASHRAE fundamentals\ Chap[ 15\ example 2\ as obtained by the method described in this paper "marked matrices#\ HVACSIM¦\ and BLAST[ Also shown are the results for the equivalent wall with tabulated coe.cients from ð2Ł and the results obtained with the single node RCR simpli_cation[
coe.cients determined by the CFTGEN procedure sup! plied with HVACSIM¦ ð15Ł\ as well as BLAST[1 All these methods give for all practical purposes the same answer[ The results are given in Fig[ 09[ In the _gure the ASH! RAE result\ which is based on an equivalent wall\ is also shown as well as the result obtained with the single node RCR lumped method[ The ASHRAE result in the _gure di}er from the results listed in the table in the book\ due to the fact that the U!factor correction for the b and c coe.cients of the equivalent wall\ as explained on page 15[07\ was not applied in the example[ To test the zone object a test routine is also available which calculates the passive response and convective load for a _xed interior temperature[ Results for the bench! mark commercial building in ASHRAE fundamentals\ Chap[ 15\ example 5 of ð2Ł are shown in Fig[ 00[ The building is a one storey small commercial building in the eastern United States\ constructed from brick with south facing windows[ The construction\ material properties and forcing functions are listed in detail in ð2\ Table 17\ p[ 15[25Ł[ To perform the simulation the interior radiant sources was distributed on the interior surfaces in the following manner ] radiation from lights and solar pen! etration*on the ~oor\ radiation from people*evenly on interior surfaces[ The ~oor {external| temperature was assumed equal to the mean outside air temperature[ In Fig[ 00 the load to maintain an interior temperature of 13>C\ as obtained with the method of this paper is marked ] matrix method[ The other traces are the Total Sensible Cooling Load "SCL#\ and Instantaneous Sensible Heat Gain "SHG#\ as listed in Table 17 of ASHRAE[ 1 Kindly computed by Prof[ J[ D[ Spitler\ School of Mechanical and Aerospace Engineering\ Oklahoma State University[
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
18
Acknowledgement The authors would like to thank Dr M[ G[ Davies\ of the School for Architecture and Building Engineering\ Liverpool\ who was asked to anonymously referee the paper by the editor|s assistant\ but decided to engage in open correspondence with the authors[ His encouraging comments and well founded critique helped clarify many important aspects and improved the paper[ We thank him for his e}ort and patience[
References Fig[ 00[ Results for the commercial building example of ASHRAE fundamentals\ Chap[ 15\ example 5 ð2Ł[ The load obtained by the method of this paper*marked matrix method*is indicated as well as the sensible heat gain "SHG# and the sensible cooling load "SCL# as calculated in the example[
Although the routines demonstrate the power and speed of the matrix method they are at present cum! bersome to use since the input data which describes the building and the loads must be supplied in a _le with a rigid format[ No post!processing to present the results are available[ A full blown implementation requires a carefully designed user interface[
3[ Conclusion A zone envelope thermal model is presented which makes no compromises except for assuming one!dimen! sional wall heat ~ow[ The conduction equation is solved exactly with matrices[ The matrices for the various sur! faces of the envelope can be combined to yield a single matrix description for the whole envelope if the interior radiant ~ux network between the walls can be simpli_ed[ This description can be combined with the heat balance or any other zone interior heat ~ow model[ A zone thermal model which includes the interior air mass and a lumped representation of the interior mass is also presented[ The lumped node for interior mass is not an essential feature*a distributed mass node can easily be incorporated[ This model is re_ned enough to be useful for simulating the response of actively controlled environ! ments\ while still simple and straightforward to explain[ The only important assumption still remaining is sim! pli_cation of the interior heat transfer\ by making use of the single interior temperature concept to characterise the interior space[ This implies "i# an assumption of mixed interior air and "ii# a less obvious simpli_cation of the radiant heat exchange between the surfaces in the room[ This last aspect will be the topic of a future investigation[
ð0Ł Achterbosch GGJ\ de Jong PPG\ Krist!Spit CE\ van der Meulen SF\ Verberne J[ The development of a convenient thermal dynamic building model[ Energy and Buildings 0874^7]072Ð85[ ð1Ł Anh TTC\ Martel A\ Camarero R[ A computer code for solving general network simulation in buildings[ In The Fourth Inter! national Symposium on the Use of Computers for Environmental Engineering Related to Buildings\ pp[ 114Ð29[ ð2Ł ASHRAE[ 0882 Fundamentals Handbook "SI#[ American Society of Heating\ Refrigerating and Air!Conditioning Engineers\ Inc[\ 0882[ ð3Ł Athienitis AK[ Application of network methods to thermal analy! sis of passive solar buildings in the frequency domain[ Ph[D[ thesis\ University of Waterloo\ Ontario\ Canada\ 0874[ ð4Ł Ayresss JM\ Stamper E[ Historical development of building energy calculations[ ASHRAE Journal 0884^36Ð44[ ð5Ł Burrus CS\ Eschenbacher PW[ An in!place\ in!order prime factor FFT algorithm[ IEEE Transactions on Acoustics\ Speech\ and Signal Processing 0870^ASSP!18"3#]795Ð06[ ð6Ł Clarke JA[ Building energy simulation ] the state!of!the!art[ Solar and Wind Technology 0878^5"3#]234Ð44[ ð7Ł Crabb JA\ Murdoch N\ Penman JM[ A simpli_ed thermal response model[ Building Serv Eng Res Technol 0876^7]02Ð8[ ð8Ł Davies MG[ An idealised model for room radiant exchange[ Build! ing and Environment 0889^14"3#]264Ð7[ ð09Ł Davies MG[ The basis for a room global temperature[ Phil Trans R Soc Lond 0881^228 A]042Ð80[ ð00Ł Davies MG[ Flaws in the environmental temperature model[ Build! ing Serv Eng Res Technol 0881^02"3#]198Ð04[ ð01Ł Davies MG[ De_nitions of room temperature[ Building and Environment 0882^17"3#]272Ð87[ ð02Ł Davies MG[ Wall parameters by time domain methods ] Part 0* response factors[ Building Serv Eng Res Technol 0884^05"2#]042Ð 6[ ð03Ł Davies MG[ Wall parameters by time domain methods ] Part 1* the conduction transfer coe.cients a\ b\ c and d[ Building Serv Eng Res Technol 0884^05"2#]048Ð53[ ð04Ł De Witt MH\ Driessen HH[ ELAN*a computer model for build! ing energy design[ Building and Environment 0877^12"3#]174Ð8[ ð05Ł Haghighat F\ Athientis A[ Comparison between time domain and frequency domain computer program for building energy analysis[ Computer!Aided Design 0877^19"8#]414Ð32[ ð06Ł Haghighat F\ Sander DM[ Experimental procedure for deter! mination of dynamic response using system identi_cation tech! niques[ Journal of Thermal Insulation 0876^00]019Ð21[ ð07Ł Kimura KI[ Scienti_c Basis of Air Conditioning[ Applied Science\ 0866[ ð08Ł Kusuda T[ Thermal response factors for multi!layered structures of various heat conduction systems[ In ASHRAE Transactions\ vol[ 64[ ASHRAE\ 0858[ 135Ð56[
29
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
ð19Ł Laret L[ Use of general models with a small number of parameters[ CLIMA 19\999 ] 6th Int Cong of Heating and Air Conditioning\ Budapest\ 0879[ ð10Ł Lord EA\ Wilson CB[ The prediction of the behaviour of building! environment systems[ Building and Environment 0872^07"0:1#]54Ð 73[ ð11Ł Mathews EH\ Richards PG[ A tool for predicting hourly air tem! peratures and sensible energy loads in buildings at sketch design stage[ Energy and Buildings 0878^03]50Ð79[ ð12Ł Mathews EH\ Richards PG\ Lombard C[ A _rst!order thermal model for building design[ Energy and Buildings 0883^10]022Ð34[ ð13Ł Minkowycz WJ\ Sparrow EM\ Schneider GE\ Pletcher RH[ Hand! book of Numerical Heat Transfer[ John Wiley and Sons\ 0877[ ð14Ł Newton D\ James R\ Bartholomew D[ Building energy simu! lation*a user|s perspective[ Energy and Buildings 0877^09]130Ð6[ ð15Ł Park C\ Clark RC\ May WB[ HVACSIM Buildings Systems and Equipment Simulation Program ] Building Loads Calculation[ Gaithersburg\ MD 19788\ 0875[ ð16Ł Rao KR\ Chandra P[ Digital computer determination of thermal frequency response of building sections[ Building Science 0855^0]188Ð296[ ð17Ł Seem JE et al[ Comprehensive room transfer functions for e.cient calculation of the transient heat transfer process in buildings[ In ] Kuehn TH\ Hickox CE "Eds[#[ Heat Transfer in Buildings and Structures\ Presented at the 13th National Heat Transfer
ð18Ł
ð29Ł
ð20Ł
ð21Ł ð22Ł
ð23Ł ð24Ł ð25Ł
Conference and Exhibition\ Pittsburgh\ Pennsylvania ] 0876^67] 24Ð34[ Seem JE et al[ Transfer functions for e.cient calculation of multi! dimensional transient heat transfer[ In ] Kuehn TH\ Hickox CE "Eds[#[ Heat Transfer in Buildings and Structures\ Presented at the 13th National Heat Transfer Conference and Exhibition\ Pitts! burgh\ Pennsylvania ] 0876^67]14Ð22[ Applied Science Division Simulation Research Group\ Center for Building Science[ DOE!1 Basics "LBL!18039#[ Lawrence Berkeley Laboratory\ University of California\ Berkeley\ CA 83619\ August 0880[ Stephenson DG and Mitalas GP[ Calculation of heat conduction transfer functions for multi!layer slabs[ ASHRAE Transactions 0860^66"1#]006Ð15[ Tindale A[ Third!order lumped parameter simulation method[ Building Serv Eng Res Technol 0882^03"2#]76Ð86[ Tuddenham D[ Computers in air!conditioning load estimation[ In ] Sherrat AFC "Ed[#[ Air!Conditioning System Design for Buildings[ McGraw!Hill\ 0871[ 85Ð000[ Van Valkenburg\ ME[ Network Analysis[ Prentice!Hall\ 0863[ Winkelman F[ Advances in building energy simulation in North America[ Energy and Buildings 0877^09]050Ð062[ Zmeureanu\ R\ Fazio P\ Haghighat F[ Analytical and inter! program validation of a building thermal model[ Energy and Build! ings 0876^09]010Ð22[
Building and Environment 23 "0888# 08Ð29
A two!port envelope model for building heat transfer C[ Lombard \ E[ H[ Mathews Department of Mechanical Engineering\ University of Pretoria\ South Africa\ 9991 Received 08 October 0884 ^ revised 08 August 0885 ^ accepted 8 October 0886
Abstract A zone envelope thermal model is derived from _rst principles\ based on the exact two!port solution to the heat conduction equation[ A single set of matrices\ one for each frequency of interest\ is used to describe the thermal behaviour of the total envelope[ The set of matrices must be computed once for a particular zone and cooling loads for various conditions\ or the passive response\ can then be obtained[ The method is based on the same principles underlying the CTF method[ However\ the transfer function is calculated for the whole envelope by combining the two!port matrices of the individual walls\ after simpli_cation of the room interior heat transfer network[ The conduction transfer coe.cients of the individual walls forming the zone are not required[ Except for the simpli_cation of interior radiation\ the method is completely exact and easy to implement[ It is also numerically e.cient[ For load prediction\ an implementation in the frequency domain\ based on the prime factor FFT algorithm is convenient[ For time domain simulation with active control\ the envelope conduction coe.cients are easily obtained by transforming to the Z!domain[ This represents a modern solution of the envelope conduction problem\ which is convenient to implement in modern object oriented computer programming languages[ Þ 0887 Elsevier Science Ltd[ All rights reserved[
Nomenclature A area of wall ðm1Ł A\ B\ C\ D elements of the transmission matrix C zone e}ective heat storage capacitance ðJ:KŁ Ca capacitance of interior air mass ðJ:KŁ Cm lumped capacitance of interior furnishings ðJ:KŁ Ea stored thermal energy in interior air mass ðJŁ Em stored thermal energy in furnishings ðJŁ F state transition matrix G state input matrix h convection coe.cient ðW:m1 = KŁ l thickness of layer ðmŁ k thermal conductance ðW:m = KŁ q heat ~ow ðWŁ Qc interior convective load ðWŁ Qcr system load for speci_ed interior temperature Tir ðWŁ Qir radiative loads on interior surfaces ðWŁ Qr radiative load ðWŁ Ra interior surface _lm resistance ðK:WŁ Rf e}ective resistance of interior forcing temperature ðK:WŁ Ri envelope e}ective interior resistance ðK:WŁ Rm resistance between interior air and furnishings ðK:WŁ Corresponding author[ Tel[ ] 99!16!01!319!320 ^ fax ] 99!16!01! 251420 ^ e!mail ] sto}el[LombardÝeng[up[ac[za
Ro envelope e}ective exterior resistance ðK:WŁ Rv ventilation resistance ðK:WŁ s Laplace domain variable ðrad:sŁ Tm transmission matrix Tmw wall transmission matrix Ti zone e}ective interior temperature ðKŁ Tir speci_ed interior temperature for system load Qcr ðKŁ To bulk temperature of outside air ðKŁ Tsa sol!air temperature of exterior surface ðKŁ Tf e}ective interior forcing temperature ðKŁ U state input vector X state vector Y thermal admittance matrix y00\ y01\ y10\ y11 elements of the admittance matrix Zo characteristic thermal impedance of a layer ðm1K:WŁ[ Greek symbols a thermal di}usivity ðm1:sŁ v radian frequency ðrad:sŁ u z"v:1a#"0¦i#l argument of two!port matrix[ Subscripts e external wall surface conditions i interior wall surface conditions s surface conditions t refers to the combined substructures[
S9259Ð0212:88:,*see front matter Þ 0887 Elsevier Science Ltd[ All rights reserved PII ] S 9 2 5 9 Ð 0 2 1 2 " 8 6 # 9 9 9 5 4 Ð 5
19
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Superscripts 0\ 1 indicating di}erent walls
0[ Introduction A very important part of a thermal simulation package for buildings is the zone thermal model[ It is also a very problematic part because the heat ~ow in a building is a complicated process involving conduction through the walls\ interior convection and radiation\ external insu! lation and convection\ penetrating radiation through fen! estration\ heat from equipment as well as interior occupant behaviour[ The graphical depiction of the vari! ous heat transfer paths by Clarke as well as his excellent discussion of the various issues are most instructive ð6Ł[ Clarke treats the development of four generations of thermal models[ These methods di}er mostly in the way they model the conduction of heat through the envelope[ In this study an extension of the highly simpli_ed model of Mathews\ Richards and Lombard ð12Ł is presented which is based on a rigorous treatment of the conduction process[ The method is based on representing the total zone envelope with one distributed heat ~ow element[ A thermal model for building designers can attempt to treat the heat ~ows in a building from _rst principles\ in a fairly rigorous\ exact manner[ This requires detailed knowledge of the _nal construction and environment ^ and also a good deal of heat transfer knowledge from the user of the model[ Another approach is to assume that few details are available and to model only the most important aspects[ This latter approach is viable if the answers so obtained are su.ciently valid to _x major passive design features\ such as window size and orien! tation[ However\ for prediction of the dynamic behaviour of HVAC systems\ the dynamics of the zone thermal response are important[ The simplest methods solve the conduction problem by lumping all the thermal mass in a single node ð0\ 04\ 7\ 19\ 21\ 11Ł\ so that a single! or second!order heat ~ow network is the result[ Solution of these models only require the integration of one or two di}erent equations\ instead of solving the partial di}er! ential equation for conduction[ The more accurate methods avoid lumping and tries to solve the conduction equation[ This article shows how a highly simpli_ed method can be extended to become a rigorous method\ without sacri_cing the bene_ts of simplicity\ by hiding the details of the conduction process[ The model of Mathews\ Richards and Lombard ð12Ł was derived with one overriding objective ^ to use a single node for all of the thermal mass[ This objective was motivated by the need to provide designers of passive buildings with a simple tool\ easy to understand\ and useful for evaluating design choices at the very beginning of the design stage when no details are available[ The validity of this highly simpli_ed method was established
Fig[ 0[ The single node model avoids the conduction equation and attempts to represent the building zone with a single heat storage capacitance[
by comparing predictions with actual measurements in 21 buildings[ The model represents the envelope with the Ro\ C and Ri elements in the network of Fig[ 0[ The other elements in the _gure are Tsa the outside e}ective sol! air temperature ðKŁ\ Qir radiative loads distributed over interior surfaces ðWŁ\ Ra interior surfaces e}ective _lm resistance ðK:WŁ\ Qc interior convective loads ðWŁ\ Rv ventilation resistance ðK:WŁ and To temperature of ven! tilating air ðKŁ[ This simple model is adequate for passive design pur! poses but it is not su.ciently detailed for simulation of the energy requirements of buildings with active climate control[ In an actively controlled building the system input may vary rapidly\ and only the innermost layer of the mass of the envelope often participates in the heat exchange[ Also\ HVAC system designers require detailed knowledge of the contribution of each surface to the total load and about the distribution of radiant energy across the interior surfaces[ For these purposes the single node model\ with no distinction between surfaces and all the mass lumped together\ is clearly inadequate[ 1[ Thermal models for buildings The prediction of the behaviour of building!environ! ment systems is discussed from _rst principles by Lord and Wilson ð10Ł[ The underlying principles are those of mass and energy conservation\ together with rate equa! tions from transport theory[ Building analysis models are necessarily simpli_cations of these fundamental prin! ciples[ Widely di}ering implementation and numerical strategies are possible[ A plethora of models and software tools exist[ Clarke discusses the state of the art in 0878 ð6Ł[ For a comparison of some North American work refer to Winkelmann ð24Ł and Tuddenham ð22Ł[ See also ð4\ 25\ 1Ł and the user|s perspective of Newton\ James and Bartholomew ð14Ł[ All these methods are attempts to solve the following two fundamental problems of thermal analysis of buildings ] 0[ The heat transmission through the building envelope must be calculated[ This requires the solution of the conduction equation\ which is a second order partial
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
10
di}erential equation in space and time coordinates ^ with appropriate boundary conditions[ 1[ The interior heat balance must be found[ The heat transfer between the interior air and interior surfaces\ the heat carried by ventilating air and air!conditioned air\ also by convective interior loads\ radiative interior loads as well as radiation penetrating through windows\ must be in balance with the energy stored in the bulk of the interior air[ The various methods di}er mostly with respect to the _rst of the two problems namely the solution of the heat conduction equation[ For the second problem only two approaches are clearly distinguishable ] the _rst is the heat balance method which solves the net of interior heat ~ows more or less exactly[ The other is to simplify the network by introducing the concept of an e}ective radi! ant or e}ective room temperature ð8\ 01\ 00\ 09Ł[ The _rst problem\ solving the conduction equation\ can be attacked by all of the very large body of analytical\ graphical and numerical methods available for solving partial di}erential equations[ Consequently\ the various thermal models for buildings are largely distinguishable by the methods and simpli_cations they employ for solv! ing the heat conduction equation[ It is probably fair to say that all the general thermal analysis programs for buildings simplify the conduction problem by assuming one!dimensional heat ~ow through walls[ This considerably reduces the problem as it allows surface temperature distributions to be approximated by a single representative temperature[ The assumption obviously fails for heat ~ow underneath buildings\ envel! ope corners and edges\ complicated envelope geometries and envelope structures with thermal bridges and special programs are available for dealing with these particular problems[ Nevertheless\ the simpli_cation to one!dimen! sional heat ~ow has been well justi_ed by experience over the past few decades[ Transverse heat transfer\ e[g[ from one part of a wall to another part of the same wall\ apparently\ plays a very small part[ If this possibility may be ignored\ it is always possible to divide a surface into smaller independent areas[ For instance\ a part of a wall which receives solar radiation through a window can be regarded as distinct from the rest of the wall which receives no radiation[ A feature of the one!dimensional conduction problem\ when applied to building walls\ is the existence of a num! ber of layers of material with di}erent thermal charac! teristics[ For strictly numerical methods this presents no complication at all ^ for analytical methods it implies the solution of the partial di}erential equation for every layer and then matching the boundary conditions[ Many approximate methods for solving the one!dimensional conduction equation change the partial di}erential equa! tion into a set of ordinary di}erential equations\ by lump! ing the distributed capacitance and resistance of the layer
Fig[ 1[ The capacitance of a massive element can be lumped to obtain an RC!ladder representation[
together at nodes\ as in Fig[ 1[ With these methods the degree of approximation is related to the number of nodes*which may be _xed or variable[ In a lumped method the temperatures between the nodes are unde! _ned[ The wall is represented with a discrete set of elements and a _nite number of temperatures[ It can be shown that most of the popular thermal models\ analytical and numerical\ make assumptions which reduce to lumping[ This is true of simpli_ed methods which explicitly lump\ and also of methods which use conduction transfer factors e[g[ the well!known ASHRAE!CTF method ð2Ł[ To construct the latter method the conduction equation was solved analytically\ and from this solution was determined a digital approxi! mation\ which is tabulated for use by designers ð08\ 20\ 02\ 03Ł[ The distributed element ideally requires an in_nite number of transfer factors[ In practice\ a _nite set is used which approximates the response of the in_nite set\ in similar manner to the approximation of the distributed wall element by a _nite number of RC!sections[ A popu! lar numerical method is to represent the partial deriva! tives of the conduction equation with _nite di}erences\ a process synonymous with lumping[ A numerical method which avoids lumping is the _nite element method[ With this method the temperature distributions of small elements are estimated by curve _tting while also ensuring continuity at the boundaries of the elements[ The _nite element method is very often used for three!dimensional heat ~ow in structures ^ but\ for computational reasons and because it is di.cult to use\ it is not favoured in building heat transfer[ For a general comparison of these methods refer to ð13Ł[ The one!dimensional heat conduction equation is an equation with space x and time t as independent vari! ables[ Lumping is a method for simplifying the spatial dimensions by concentrating the heat stored in the struc! ture at a few points in space[ A di}erent approach is to simplify the time coordinate instead[ The conduction equation is a linear equation\ hence complicated solutions can be constructed by superimposing simpler solutions[ Fourier methods are based on Fourier|s theorem that all
11
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
periodic solutions can be constructed from sinusoids[ Fourier|s theorem allows the construction of arbitrary periodic responses by adding the responses of the Fourier components[ Since no simpli_cation of the spatial coor! dinate is assumed the procedure is exact for periodic functions[ The method calculates the response under the assumption that the time coordinate is periodic[ This assumption is well suited to the quasi!periodic time dependencies experienced in practice[ Non!periodic responses can be calculated to any degree of accuracy by assuming a su.ciently long period[ However\ for very long periods\ and for short sudden chances\ it is less attractive from a computational point of view ^ requiring a large number of frequencies to be taken into account[ The suitability of this frequency domain approach has been most vigorously argued by Haghighat and Athi! enithis ð05Ł and long before them by Rao and Chandra ð16Ł[
Fig[ 3[ The two!port representation of a sandwich structure is simply the product of the two!port matrices for each layer[
2[ The matrix method The matrix method uses the solution for the con! duction equation with sinusoidal forcing functions[ The construction of the solution is well described in the litera! ture ð07\ p[ 49Ł\ ð08\ 20Ł[ It represents conduction through a layer with a 1×1 matrix as in Fig[ 2[ Mathematically the solution to the conduction equation is written in the form of a matrix operator which transform conditions at one surface into the corresponding conditions at the other surface[ With T0 the temperature and q0 the heat ~ow at surface 0 and T1\ q1 the corresponding quantities on the other the two!port model are T1
$ % $ q1
A B C
D
T0
%$ % ×
q0
"0#
where the two!port cascade "transmission# matrix is A cosh u B Zo sinh u C
0 sinh u D cosh u[ Zo
Fig[ 2[ Representing a massive layer with a two!port matrix[
"1#
In this equation u z"v:1a#"0¦i#l and Zo l:Aku are parameters of the layer which depend on the frequency v and physical properties ] A is the wall area\ l is the layer thickness and k is the thermal conductivity of the layer[ The matrix solution is very convenient for a layered structure[ The matrix for a combination of layers is sim! ply the product of the matrices for the individual layers as in Fig[ 3[ The matrix solution to the conduction equation characterises the heat transfer with four complex elements[ Of this four only three are independent and the fourth element is _xed by the condition that the determinant of the matrix be unity[ The CTF method\ as described in ASHRAE Fundamentals ð2Ł\ uses two of the elements of the matrix to calculate conduction from the external surface to the interior of the wall and to the wall interior surface[ To model the ~ux of heat from the interior surface to the interior mass of the wall the other two matrix elements are required[ Here the CTF method is inconsistent ] it combines the walls and uses a room transfer function "RTF# to take into account the storage of penetrating radiation in the walls\ while determining the sensible heat gain by conduction through the walls individually for each wall[ Besides being convenient for solving heat ~ow through a layered structure\ the matrices can also be used for adding the contributions from di}erent walls in a con! sistent manner[ If the model is consistent it will allow the decomposition of the total response into the con! tributions of the individual walls[ With this method the matrix for the whole zone envelope is found once and the matrix is then used as a complete representation of the envelope for all further calculations[ In the next section it is shown that a single 1×1 frequency dependent matrix characterises the whole building envelope and that it pro!
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
vides an exact solution for the heat ~ow through the envelope[ Most of the more accurate building thermal simulation programs are based on the matrix method[ The ASH! RAE!CTF method and programs like DOE1e\ BLAST\ HVACSIM¦ use the matrix method to calculate the heat transfer through a wall and from this exact solution an approximate di}erence equation is extracted[ This di}erence equation can be determined by approximating either the frequency response or the impulse response of the matrices ð20Ł[ The di}erence equation steps from the frequency domain back into the time domain[ For a speci_ed wall the CTF|s are determined beforehand and can then be used to determine the response for any\ known beforehand or not\ forcing functions[ ASHRAE Fundamentals tabulates the sets of CTF values of bn and dn for various walls and the user may look them up and apply them directly[ The CTF method treats each wall on its own and combines the heat!~ows through the walls with other loads to _nd a heat load to the interior air\ the conduction heat gain "CHG#[ The cooling includes the conduction heat paths to the interior air\ direct con! vective loads and also radiation penetrating through the windows and striking the interior surfaces[ This heat is initially absorbed by the structure and slowly released to the interior air[ The CTF procedure takes this into account by the use of another di}erence equation\ the room transfer function "RTF# ð17Ł[ The room transfer function represents the combined e}ect of all the walls of the zone[ The sets of RTF values of v and w are also tabulated in ASHRAE for a small number of rep! resentative zones[ An argument against the Fourier series method is that it cannot be applied to systems simulation[ In simulation a controller intervenes to regulate the system[ The inter! vention is based on the instantaneous interior tem! perature so that the response of the system is not known beforehand[ For Fourier analysis\ the forcing functions must be known beforehand to be able to calculate the spectra[0 But\ for passive analysis\ and to _nd the required system response that will maintain a speci_ed interior temperature\ the frequency domain method is accurate and easy to implement[ To simulate ~oating interior tem! peratures it is better to transform to the time domain and to characterise the envelope with a di}erence equation[ The ASHRAE!CTF method\ as presented in the Fun! damentals ð2\ Chap[ 15Ł\ is also limited in this regard[ The tabulated coe.cients do not include the full set required 0 This argument is not entirely valid[ For some control laws it is possible to solve the set of equations in the frequency domain and to determine both the response of the building and that of the system simultaneously[ This is sometimes done by control engineers to _nd a control law which optimises the system response ^ e[g[ minimum energy controller[ However\ solving the equations for actual systems is not often feasible because of the strong non!linearities in practical control laws[
12
Fig[ 4[ Equivalent building zone thermal model with walls represented by a two!port matrix[
for a time dependent interior temperature\ only the summed c series "Sci# is given[ This means the interior temperature must be constant[ The simulation program HVACSIM¦ calculates the full set of CTF values via the Laplace transform method[ This means the solution for each wall is obtained in the frequency domain and then transformed to the time domain\ for solving the interior energy balance[ Davies ð03Ł also presents a purely time domain method for obtaining the full set[
2[0[ A two!port zone model The building thermal models used by HVACSIM¦ and in the ASHRAEÐCTF method are complicated networks of energy ~ow through the walls and between the interior surfaces of the walls[ It is much easier to explain building heat transfer in terms of a single envelope through which heat enters or leaves the zone[ This approach emphasises the importance of designing a building shell with the correct overall thermal properties[ To obtain this simple representation of building heat ~ow it is necessary to _nd a single element description of the envelope and to combine the forcing functions into external and internal e}ective functions[ This can be accomplished without sacri_cing any accuracy in the manner explained below[ Figure 4 is an extension of Fig[ 0[ In the _gure a single two!port matrix characterises the envelope[ The forcing functions on the exterior surfaces of the walls are com! bined in a single e}ective forcing function Tsa[ The matrix includes both exterior and interior _lm coe.cients[ Ti represents the e}ective interior temperature of the zone with a single node[ The other elements are ] Ca which represents the lumped mass of the interior air\ Cm the lumped mass of furniture and other interior structures\ Rm the thermal resistance between the interior e}ective temperature and the interior mass\ Tf and Rf represent direct heat ~ow to the interior air caused by ventilation\ air conditioning and convective heating as explained in Fig[ 5 Tf
RvRsQc¦RsTo¦RvTs Rs¦Rv
"2#
13
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Fig[ 5[ De_nition of direct heat gain to interior air source and resistance[
Rf
RvRs [ Rv¦Rs
"3#
Qr represents an e}ective radiative heat ~ow to the interior surface of the envelope[ Ra is the e}ective _lm resistance between the interior envelope surface and the e}ective interior node[ The de_nitions of Ra and Ti are based on simpli_cation of the interior heat transfer[ Dav! ies ð09Ł has published an extensive analysis of the issues involved[ He presents a model which characterises the room interior with two temperatures\ the rad!air model\ which is logically founded and allows the calculation of interior comfort temperature[ Furthermore\ it enables additional simpli_cation by allowing all the interior radi! ation sources to be connected only to the two interior nodes[ For the purposes of this paper it is su.cient to assume that it is possible to use this representation with some degree of accuracy[ In our implementation we have used the approximation Ra 0:"7A#\ where A is the area of the wall\ and we assume that Ti will closely correspond with the interior air temperature[ The inclusion of a more re_ned model awaits further investigation[ In Fig[ 4 the resistance Ra is shown absorbed in the two!port for the envelope[ The model in Fig[ 4 is a direct extension of the single node model of Mathews\ Richards and Lombard ð12Ł as in Fig[ 0[ It represents an attempt to retain the simplicity and ease of understanding of the single node model but to extend its accuracy and applicability by re_ning the model of the envelope and the inclusion of an interior air mass node and a separate node for interior massive structures[ Although it is convenient to have a single matrix rep! resentation for a zone envelope it is not a requirement[ It is possible to retain the matrix description for individual surfaces and with these matrices construct a heat ~ow graph which can be solved with standard network methods[ This is the approach described by Athienitis ð3Ł[ His method is a full frequency domain method which represents the zone with a network of admittances[ The solution procedure involves _nding an admittance matrix which combines the heat ~ows from the various walls with other loads to determine the _nal response[ This admittance matrix characterises the whole zone and not just the envelope[ In Figure 4 the simpli_cations are ]
0[ Massive elements in the interior of the zone are rep! resented with a lumped capacitance and an air contact resistance[ This approximation is not essential[ Anal! ogous to the two!port representation of the envelope a one!port representation of interior mass is possible[ 1[ The radiant energy exchange between interior surfaces is simpli_ed[ A single node is used which is involved in all energy exchanges\ both convective and radiative[ This simpli_cation is often employed since it allows the simpli_cation of the network of interior heat ~ows as in Fig[ 6 to a single node[ For a full discussion of its bene_ts and problems see Davies ð01\ 09Ł[ The same approximation is used to DOE1e ð29Ł to calculate RTF coe.cients and in the ASHRAE!CTF method[ In this model no lumping of envelope elements takes place[ The envelope is represented by the two!port matrix[ The outside forcing function is also completely exact[ The envelope can always be decomposed again into individual walls and the contributions of the walls will correctly sum to the total response[ However\ in the numerical procedures for solving the model of Fig[ 4 inaccuracies may arise[ These inaccuracies have no direct relation with the heat ~ow model[ For calculating the passive response to given outside con! ditions and speci_ed internal loads\ and for _nding the required cooling load for a speci_ed interior temperature\ a full frequency domain method is convenient[ The full frequency domain method _nds a Fourier series rep! resentation of the forcing functions\ calculates a matrix for each frequency and combine the responses for each frequency to obtain the _nal response[ For simulation and intelligent control a time domain approach is more convenient[ This requires transforming the matrix to a set of di}erence equations which can be used in a similar way as in the CTF method[ However\ in contrast with the CTF method\ a consistent approach is followed[ The various elements of the same matrix are used both for calculating the conduction through the walls and the response to interior radiative loads[ In mathematical terms the di}erence equations represent a transform of the matrix operator\ which represents the envelope\ from the Fourier domain to the Z!domain\ and all the tech! niques for accomplishing this are applicable[
Fig[ 6[ Transforming the interior heat ~ow network to the star rep! resentation yields a single node involved in all heat transfers[
14
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
2[1[ The envelope matrix
Ti qiRa¦Ts
In Fig[ 4 the envelope is represented with a single matrix[ It is convenient if this matrix includes both sur! face _lm resistances[ The radiative sources\ which work directly on the surfaces and not through the _lm resist! ances\ are then connected to points which are interior to the matrix[ Also\ the forcing functions on the surfaces of the walls are di}erent for each surface[ The single matrix representation requires that the forcing functions be reduced to two e}ective forcing functions ^ one for the outside and another for the inside surface of the envelope[ This reduction of the number of forcing functions is a crucial aspect of the simple representation[ On the outside surface the forcing functions can be combined in a single temperature forcing function by using the well established Sol!Air concept[ The external solar!\ long wave!radiation and convection with outside air is combined by _nding the Thevenin equivalent circuit of these forcing functions[ The Thevenin source is then the Sol!Air temperature and the Thevenin resistance the e}ective _lm resistance which can be included in the matrix[ For the interior surfaces\ which are not termin! ated\ but will be connected to other surfaces and the interior air node it is not so simple[ In Fig[ 7 the model for a single wall is given[ In this _gure the exterior e}ective "Thevenin# temperature is rep! resented by Te and the exterior surface _lm resistance has been absorbed into the two!port for the wall\
and
Tmw
$
Aw Bw Cw
%
Dw
Ts
$% $ qs
Aw Bw Cw
Ti
×
"qs−Qr#Ra¦Ts
$% $ qi
Te qe
0 Ra
9
×
0
Ts
Ra
−Qr
qs
0
[
"7#
But with eqn "4# conditions on the interior surface can be written in terms of exterior conditions ] Ti
0 Ra
$% $ %$ qi
9
×
0
Aw B w
Te
Ra
%$%
%$Cw
×
Dw
qe
−Qr
"8#
0
showing that despite the presence of the source Qr\ the matrix for the wall and interior _lm coe.cient is still the product of the two individual matrices and both _lm resistances can be included in the de_nition of the wall! matrix Tm
$
A B C
0 Ra
% $ %$
D
9
0
×
Aw Bw Cw
Dw
%
"09#
[
The _nal solution for the wall is
A B C
D
Te
Ra
%$% $% ×
qe
−Qr
0
"00#
[
To determine the heat ~ows through the surfaces of the wall eqn "00# must be solved for qe and qi in terms of the surface temperatures[ This gives qe
−ATe¦Ti¦RaQr 0 B "BC−AD#Te¦DTi¦"DRa−B#Qr
$% $ qi
%
"01#
which can conveniently be written qe
"4#
% $ %$%
%$qs−Qr
¦
qi
%
%$Dw
Hence
$% $
The heat ~ow across the exterior surface is qe[ Ts is the temperature on the inside surface of the wall and qs the heat ~ow across the inside surface[ Ra is the _lm resistance between the interior surface and the interior air\ which is at temperature Ti[ The heat ~ow to the interior surface from the interior air is qi[ Qr is radiation heat ~ow to the interior surface[ The matrix solution for the interior surface of the wall gives ]
"6#
qi qs−Qr[
Ti
[
"5#
−A 0 0 B BC−AD D
$% $ qi
%$
Te Ti¦QrRa
9
% $ % ¦Qr
−0
[ "02#
and from the network
In eqn "02# we see that the e}ective interior wall surface forcing temperature is Ti¦QrRa[ However\ it is more convenient to keep the forcing functions separate and to write qe
−A 0 0 B BC−AD D
$% $ qi
Fig[ 7[ The thermal model for a wall consists of a two!port matrix and exterior and interior forcing functions as de_ned in the text[
Te
Ra Qr [ B DRa−B
%$ % $ Ti
¦
%
"03#
In eqn "03# the two!port matrix for the wall is in the admittance or Y format[ If the B element of the Tm matrix is non!zero it can always be transformed to the
15
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
Y
$
y00 y01
qet
%
qit
format with y00 −
A 0 y01 B B "04#
y00 0 B y01 y01 y00y11 y11 D [ y01 y01
"05#
Qre
Ra y10Ra Qr Qr B DRa−B y11Ra−0
$ % $ Qri
% $
%
"06#
eqn "03# can be written compactly "07#
q YT¦Qr with the further de_nitions q
qe
$% qi
qi1
"19#
and T
Te
$% Ti
[
"10#
where subscript t is used to indicate variables which are representative of the combined structure[ This important equation gives the contribution of the various forcing functions on the individual walls to the total heat ~ow across the envelope[ To _nd a single two!port representation for the envel! ope eqn "10# must agree with qt YtTt¦Qrt[
"11#
This implies that
The T! and Y!representation of the two!port are math! ematically equivalent[ The T!representation is convenient when systems are in cascade while the Y!representation is used for systems in parallel where the ~ows add[ If a radiative source vector is de_ned by Qr
qe1
qt Y0T0¦Qr0¦Y1T1¦Qr1
If y01 is non!zero the inverse transform is
C y10−
qi0
¦
or
AD D y11 [ y10 C− B B
A −
qe0
$% $ %$ %
y10 y11
"08#
The building envelope consists of a number of surfaces which enclose the interior space[ The envelope heat ~ow is the sum of the heat ~ows through all the surfaces[ In Fig[ 8 the envelope\ for simplicity\ has only two walls[ The total heat ~ow through the envelope is
Yt Y0¦Y1
"12#
Yt Y−0 t "Y0T0¦Y1T1#
"13#
Qrt Qr0¦Qr1
"14#
which de_nes the e}ective envelope Y!matrix and forcing functions[ The equation for the combined heat ~ow of two walls is in exactly the same form as that for a single wall[ This makes it possible to extent the model from two walls to any number of walls by accumulating the walls one by one into a total representation for the whole envelope[ The Y!matrices of the walls are added together and the e}ective forcing functions are calculated from a weighted average where the weights are also obtained from the elements of the Y!matrices[ 2[2[ Solution With the envelope matrix known the circuit of Fig[ 4 can be solved with the envelope regarded as a black box described by eqn "11#[ For passive analysis the interior temperature is the unknown variable[ For calculating the HVAC system load the heat ~ow to the interior air for a speci_ed interior air temperature is the unknown[ To solve the network of Fig[ 4 standard nodal analysis ð23Ł is used[ For state variables we choose the thermal energy in the interior air mass Ea TiCa and the thermal energy in the massive interior structures Em TmCm[ At node i the heat ~ows must add to zero\ hence qi¦qa¦qm¦qf 9[
"15#
But qa
Fig[ 8[ Thermal model for an envelope consisting of only two surfaces[
d Ea dt
"16#
d Ti−Tm Em dt Rm
"17#
qm
16
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
qf
Ti−Tf Rf
"18#
and qi is supplied by eqn "11# "29#
qi y10Tsa¦y11Ti¦Qri
where the variables now refer to the totals\ summed for all walls forming the envelope[ It is convenient to de_ne qr y10Tsa¦Qri qi−y10Ti\
"20#
a quantity which is known as it involves only the known forcing functions[ If these are substituted in eqn "15# we have
0
qr¦y11
1
0
1
0 0 d d 0 E a ¦ Ea ¦ Em ¦ Ea− Tf 9[ Ca dt dt CaRf Rf "21#
At node m the heat balance is ] qm
d Ti−Tm 0 0 Ea− Em[ Em dt Rm CaRm CmRm
G
0
0
Ea
0
1
1
0 0 Em− Tf 9[ CmRf Rf
X
$ % Em
and U
$
Tf:Rf−qr 9
%
"39#
"25#
"30#
Tf−To −Qc Rf
Qcr
"31#
where "32#
and qf −"qi¦qa¦qm# qi qr¦y11Tir
d X FX¦GU dt
"26#
where
1
0 0 K −0 H Ca y11¦ Rm ¦ Rf F H H 0 k CaRm and
1
[
are required from the envelope so that only two matrix elements are calculated and stored[ In the above equa! tions Tf\ Ea\ qr\ Qa\ Tsa\ y11 and y10 are complex numbers giving the amplitude and phase at frequency v[ It is also possible to _nd the required extra convective load Qcr which will ensure a given interior temperature Tir[
Tf Tir−Rfqf
the state equation in standard vector form is
0
0
0 Rf 0¦Rmy11¦sRmCa− ¦Rm sRmCm¦0
"24#
Equations "24# and "22# are the state equations for the circuit[ If the state vector X\ and the source vector U\ are de_ned by Ea
CaRm"Tf−Rfqr#
qr y10 ( Tsa¦"y11Ra−0#Qa and y11 "23#
0 0 0 d y11¦ ¦ Ea¦ Ea Ca Rm Rf dt −
"28#
[
In this solution only
This equation simpli_es to q r¦
9 0
"22#
1
0 0 E − T 9[ CaRf a Rf f
¦
$ %
This state eqn "26# is the complete description of the zone[ It|s solution\ for given forcing functions\ gives the response of the zone to the forcing functions in terms of the energy stored in the air!mass and interior mass[ The format of eqns "25# and "26# is misleading since it indi! cates a second order system[ In fact\ the envelope oper! ators y10 and y11 are high order di}erential operators so that the actual order is much higher[ Equation "26# re~ects the simplicity of Fig[ 4[ Note that energy stored in the structure is not visible in the equation because of the black box hiding the complexity of the conduction process\ which is treated in a rigorously exact manner[ To calculate the passive temperature response of the zone eqn "26# must be integrated to _nd Ea and Em[ The interior air!temperature is then Ti Ea:Ca[ In the frequency domain this is straightforward[ Writing s jv for the frequency dependence and solving eqn "26#
Substituting this result in eqn "21# gives 0 y11 d 0 qr¦ Ea¦ Ea¦ Ea− E Ca dt CaRm CmRm m
0 9
0 L CmRm H H −0 H l CmRm
"27#
qa sCaTir qm
Tir−Tm Rm
Tm
Tir 0¦sRmCm
"33#
with qr as de_ned above[ The frequency domain solution involves ] calculating and storing the two envelope matrix operators y10"v# and
17
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
y11"v#\ _nding the frequency components of the forcing functions and calculating the frequency components of the solution through 39 or 31[ For a time domain implementation the y10 and y11 operators must be trans! lated to di}erence operators[ Various standard tech! niques exist as discussed in ð20\ 08\ 17\ 18\ 06Ł[ This leads to a method similar to the CTF!method\ except that the CTF!coe.cients are calculated for the whole envelope from the accurate matrix description of the envelope\ instead of for each surface[ 2[3[ Implementation For implementation of the scheme a modern object oriented approach is ideal[ In an OOP language such as C¦¦ a software object is created for each physical object[ Software objects are also de_ned to store and manipulate the data structures[ We have implemented an object hier! archy consisting of the following data objects ] 0[ A two!port object which allows matrix algebra and T \ Y transforms through operator overloading[ 1[ A time series object which can retrieve\ store and print a time sequence\ and transform it to the frequency domain using a mixed radix FFT algorithm suitable for sequences with length a multiple of 13 ð5Ł[ 2[ A spectral object which can store and print a sequence of frequency components and transform them to a time series representation with the same FFT algo! rithm[ The software representation of the zone consists of the following {physical| objects ] 0[ A layer object which contains material and dimen! sional data for a wall layer[ 1[ A surface object which contains and combines a sand! wich structure of layers with an external Sol!Air tem! perature series and an interior radiative load series[ 2[ An envelope object which combines a number of sur! faces in a zone envelope and _nds the equivalent envel! ope description[ 3[ A zone object which combines the envelope with the interior air\ interior mass\ ventilation and convective forcing functions in the manner indicated in Fig[ 4[ To verify the code exercise routines were written for each object and the answers were thoroughly checked by hand calculations[ All the calculations were veri_ed in this way by comparing the mean response and the swing for a single sinusoid forcing function[ The surface object was veri_ed by calculating the heat ~ow for the wall object of example 2\ Chap[ 15\ in ð2Ł[ In this example the cooling load is calculated to maintain a temperature of 13>C behind a wall constructed of 099 mm heavy concrete\ with 49 mm insulation\ and 19 mm indoor plaster ^ for a prescribed external Sol!Air temperature[ The same struc! ture was also calculated with the CTF method\ with
Fig[ 09[ Results for the wall of ASHRAE fundamentals\ Chap[ 15\ example 2\ as obtained by the method described in this paper "marked matrices#\ HVACSIM¦\ and BLAST[ Also shown are the results for the equivalent wall with tabulated coe.cients from ð2Ł and the results obtained with the single node RCR simpli_cation[
coe.cients determined by the CFTGEN procedure sup! plied with HVACSIM¦ ð15Ł\ as well as BLAST[1 All these methods give for all practical purposes the same answer[ The results are given in Fig[ 09[ In the _gure the ASH! RAE result\ which is based on an equivalent wall\ is also shown as well as the result obtained with the single node RCR lumped method[ The ASHRAE result in the _gure di}er from the results listed in the table in the book\ due to the fact that the U!factor correction for the b and c coe.cients of the equivalent wall\ as explained on page 15[07\ was not applied in the example[ To test the zone object a test routine is also available which calculates the passive response and convective load for a _xed interior temperature[ Results for the bench! mark commercial building in ASHRAE fundamentals\ Chap[ 15\ example 5 of ð2Ł are shown in Fig[ 00[ The building is a one storey small commercial building in the eastern United States\ constructed from brick with south facing windows[ The construction\ material properties and forcing functions are listed in detail in ð2\ Table 17\ p[ 15[25Ł[ To perform the simulation the interior radiant sources was distributed on the interior surfaces in the following manner ] radiation from lights and solar pen! etration*on the ~oor\ radiation from people*evenly on interior surfaces[ The ~oor {external| temperature was assumed equal to the mean outside air temperature[ In Fig[ 00 the load to maintain an interior temperature of 13>C\ as obtained with the method of this paper is marked ] matrix method[ The other traces are the Total Sensible Cooling Load "SCL#\ and Instantaneous Sensible Heat Gain "SHG#\ as listed in Table 17 of ASHRAE[ 1 Kindly computed by Prof[ J[ D[ Spitler\ School of Mechanical and Aerospace Engineering\ Oklahoma State University[
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
18
Acknowledgement The authors would like to thank Dr M[ G[ Davies\ of the School for Architecture and Building Engineering\ Liverpool\ who was asked to anonymously referee the paper by the editor|s assistant\ but decided to engage in open correspondence with the authors[ His encouraging comments and well founded critique helped clarify many important aspects and improved the paper[ We thank him for his e}ort and patience[
References Fig[ 00[ Results for the commercial building example of ASHRAE fundamentals\ Chap[ 15\ example 5 ð2Ł[ The load obtained by the method of this paper*marked matrix method*is indicated as well as the sensible heat gain "SHG# and the sensible cooling load "SCL# as calculated in the example[
Although the routines demonstrate the power and speed of the matrix method they are at present cum! bersome to use since the input data which describes the building and the loads must be supplied in a _le with a rigid format[ No post!processing to present the results are available[ A full blown implementation requires a carefully designed user interface[
3[ Conclusion A zone envelope thermal model is presented which makes no compromises except for assuming one!dimen! sional wall heat ~ow[ The conduction equation is solved exactly with matrices[ The matrices for the various sur! faces of the envelope can be combined to yield a single matrix description for the whole envelope if the interior radiant ~ux network between the walls can be simpli_ed[ This description can be combined with the heat balance or any other zone interior heat ~ow model[ A zone thermal model which includes the interior air mass and a lumped representation of the interior mass is also presented[ The lumped node for interior mass is not an essential feature*a distributed mass node can easily be incorporated[ This model is re_ned enough to be useful for simulating the response of actively controlled environ! ments\ while still simple and straightforward to explain[ The only important assumption still remaining is sim! pli_cation of the interior heat transfer\ by making use of the single interior temperature concept to characterise the interior space[ This implies "i# an assumption of mixed interior air and "ii# a less obvious simpli_cation of the radiant heat exchange between the surfaces in the room[ This last aspect will be the topic of a future investigation[
ð0Ł Achterbosch GGJ\ de Jong PPG\ Krist!Spit CE\ van der Meulen SF\ Verberne J[ The development of a convenient thermal dynamic building model[ Energy and Buildings 0874^7]072Ð85[ ð1Ł Anh TTC\ Martel A\ Camarero R[ A computer code for solving general network simulation in buildings[ In The Fourth Inter! national Symposium on the Use of Computers for Environmental Engineering Related to Buildings\ pp[ 114Ð29[ ð2Ł ASHRAE[ 0882 Fundamentals Handbook "SI#[ American Society of Heating\ Refrigerating and Air!Conditioning Engineers\ Inc[\ 0882[ ð3Ł Athienitis AK[ Application of network methods to thermal analy! sis of passive solar buildings in the frequency domain[ Ph[D[ thesis\ University of Waterloo\ Ontario\ Canada\ 0874[ ð4Ł Ayresss JM\ Stamper E[ Historical development of building energy calculations[ ASHRAE Journal 0884^36Ð44[ ð5Ł Burrus CS\ Eschenbacher PW[ An in!place\ in!order prime factor FFT algorithm[ IEEE Transactions on Acoustics\ Speech\ and Signal Processing 0870^ASSP!18"3#]795Ð06[ ð6Ł Clarke JA[ Building energy simulation ] the state!of!the!art[ Solar and Wind Technology 0878^5"3#]234Ð44[ ð7Ł Crabb JA\ Murdoch N\ Penman JM[ A simpli_ed thermal response model[ Building Serv Eng Res Technol 0876^7]02Ð8[ ð8Ł Davies MG[ An idealised model for room radiant exchange[ Build! ing and Environment 0889^14"3#]264Ð7[ ð09Ł Davies MG[ The basis for a room global temperature[ Phil Trans R Soc Lond 0881^228 A]042Ð80[ ð00Ł Davies MG[ Flaws in the environmental temperature model[ Build! ing Serv Eng Res Technol 0881^02"3#]198Ð04[ ð01Ł Davies MG[ De_nitions of room temperature[ Building and Environment 0882^17"3#]272Ð87[ ð02Ł Davies MG[ Wall parameters by time domain methods ] Part 0* response factors[ Building Serv Eng Res Technol 0884^05"2#]042Ð 6[ ð03Ł Davies MG[ Wall parameters by time domain methods ] Part 1* the conduction transfer coe.cients a\ b\ c and d[ Building Serv Eng Res Technol 0884^05"2#]048Ð53[ ð04Ł De Witt MH\ Driessen HH[ ELAN*a computer model for build! ing energy design[ Building and Environment 0877^12"3#]174Ð8[ ð05Ł Haghighat F\ Athientis A[ Comparison between time domain and frequency domain computer program for building energy analysis[ Computer!Aided Design 0877^19"8#]414Ð32[ ð06Ł Haghighat F\ Sander DM[ Experimental procedure for deter! mination of dynamic response using system identi_cation tech! niques[ Journal of Thermal Insulation 0876^00]019Ð21[ ð07Ł Kimura KI[ Scienti_c Basis of Air Conditioning[ Applied Science\ 0866[ ð08Ł Kusuda T[ Thermal response factors for multi!layered structures of various heat conduction systems[ In ASHRAE Transactions\ vol[ 64[ ASHRAE\ 0858[ 135Ð56[
29
C[ Lombard\ E[H[ Mathews:Buildin` and Environment 23 "0888# 08Ð29
ð19Ł Laret L[ Use of general models with a small number of parameters[ CLIMA 19\999 ] 6th Int Cong of Heating and Air Conditioning\ Budapest\ 0879[ ð10Ł Lord EA\ Wilson CB[ The prediction of the behaviour of building! environment systems[ Building and Environment 0872^07"0:1#]54Ð 73[ ð11Ł Mathews EH\ Richards PG[ A tool for predicting hourly air tem! peratures and sensible energy loads in buildings at sketch design stage[ Energy and Buildings 0878^03]50Ð79[ ð12Ł Mathews EH\ Richards PG\ Lombard C[ A _rst!order thermal model for building design[ Energy and Buildings 0883^10]022Ð34[ ð13Ł Minkowycz WJ\ Sparrow EM\ Schneider GE\ Pletcher RH[ Hand! book of Numerical Heat Transfer[ John Wiley and Sons\ 0877[ ð14Ł Newton D\ James R\ Bartholomew D[ Building energy simu! lation*a user|s perspective[ Energy and Buildings 0877^09]130Ð6[ ð15Ł Park C\ Clark RC\ May WB[ HVACSIM Buildings Systems and Equipment Simulation Program ] Building Loads Calculation[ Gaithersburg\ MD 19788\ 0875[ ð16Ł Rao KR\ Chandra P[ Digital computer determination of thermal frequency response of building sections[ Building Science 0855^0]188Ð296[ ð17Ł Seem JE et al[ Comprehensive room transfer functions for e.cient calculation of the transient heat transfer process in buildings[ In ] Kuehn TH\ Hickox CE "Eds[#[ Heat Transfer in Buildings and Structures\ Presented at the 13th National Heat Transfer
ð18Ł
ð29Ł
ð20Ł
ð21Ł ð22Ł
ð23Ł ð24Ł ð25Ł
Conference and Exhibition\ Pittsburgh\ Pennsylvania ] 0876^67] 24Ð34[ Seem JE et al[ Transfer functions for e.cient calculation of multi! dimensional transient heat transfer[ In ] Kuehn TH\ Hickox CE "Eds[#[ Heat Transfer in Buildings and Structures\ Presented at the 13th National Heat Transfer Conference and Exhibition\ Pitts! burgh\ Pennsylvania ] 0876^67]14Ð22[ Applied Science Division Simulation Research Group\ Center for Building Science[ DOE!1 Basics "LBL!18039#[ Lawrence Berkeley Laboratory\ University of California\ Berkeley\ CA 83619\ August 0880[ Stephenson DG and Mitalas GP[ Calculation of heat conduction transfer functions for multi!layer slabs[ ASHRAE Transactions 0860^66"1#]006Ð15[ Tindale A[ Third!order lumped parameter simulation method[ Building Serv Eng Res Technol 0882^03"2#]76Ð86[ Tuddenham D[ Computers in air!conditioning load estimation[ In ] Sherrat AFC "Ed[#[ Air!Conditioning System Design for Buildings[ McGraw!Hill\ 0871[ 85Ð000[ Van Valkenburg\ ME[ Network Analysis[ Prentice!Hall\ 0863[ Winkelman F[ Advances in building energy simulation in North America[ Energy and Buildings 0877^09]050Ð062[ Zmeureanu\ R\ Fazio P\ Haghighat F[ Analytical and inter! program validation of a building thermal model[ Energy and Build! ings 0876^09]010Ð22[