Inside Heat Balance[LINK]

The heart of the heat balance method is the internal heat balance involving the inside faces of the zone surfaces. This heat balance is generally modeled with four coupled heat transfer components: 1) conduction through the building element, 2) convection to the air, 3) short wave radiation absorption and reflectance and 4) longwave radiant interchange. The incident short wave radiation is from the solar radiation entering the zone through windows and emittance from internal sources such as lights. The longwave radiation interchange includes the absorption and emittance of low temperature radiation sources, such as all other zone surfaces, equipment, and people.

The heat balance on the inside face can be written as follows:

$$q_{LWX}^{\prime \prime }+q_{SW}^{\prime \prime }+q_{LWS}^{\prime \prime }+q_{ki}^{\prime \prime }+q_{sol}^{\prime \prime }+q_{conv}^{\prime \prime }=0$$

where:

$q_{LWX}^{\prime \prime }$ = Net longwave radiant exchange flux between surfaces in a zone or group of zones (enclosure).

$q_{SW}^{\prime \prime }$ = Net short wave radiation flux to surface from lights.

$q_{LWS}^{\prime \prime }$ = Longwave radiation flux from equipment in a zone or group of zones (enclosure).

$q_{ki}^{\prime \prime }$ = Conduction flux through the wall.

$q_{sol}^{\prime \prime }$ = Transmitted solar radiation flux absorbed at surface.

$q_{conv}^{\prime \prime }$ = Convective heat flux to zone air.

Each of these heat balance components is introduced briefly below.

Internal Long-Wave Radiation Exchange[LINK]

Long-wave radiation exchange is balanced within an enclosure which may be a single zone or a group of zones connected by air boundaries (see Construction:AirBoundary). Throughout this section, “Zone” refers to either type of enclosure.

LW Radiation Exchange Among Zone Surfaces[LINK]

There are two limiting cases for internal LW radiation exchange that are easily modeled:

• The zone air is completely transparent to LW radiation.

• The zone air completely absorbs LW radiation from the surfaces within the zone.

The limiting case of completely absorbing air has been used for load calculations and also in some energy analysis calculations. This model is attractive because it can be formulated simply using a combined radiation and convection heat transfer coefficient from each surface to the zone air. However, it oversimplifies the zone surface exchange problem, and as a result, the heat balance formulation in EnergyPlus treats air as completely transparent. This means that it does not participate in the LW radiation exchange among the surfaces in the zone. The model, which considers room air to be completely transparent, is reasonable physically because of the low water vapor concentrations and the short mean path lengths. It also permits separating the radiant and convective parts of the heat transfer at the surface, which is an important attribute of the heat balance method.

EnergyPlus uses a grey interchange model for the longwave radiation among zone surfaces. EnergyPlus offers two algorithms for modeling long wave radiation: The “ScriptF” method, and the “CarrollMRT” methods. Users can select between these two algorithms using the “PerformancePrecisionTradeoffs” object.

ScriptF[LINK]

The “ScriptF” algorithm was developed by Hottel (Hottel and Sarofim, Radiative Transfer, Chapter 3, McGraw Hill, 1967). This procedure relies on a matrix of exchange coefficients between pairs of surfaces that include all exchange paths between the surfaces. In other words all reflections, absorptions and re-emissions from other surfaces in the enclosure are included in the exchange coefficient, which is called ScriptF. The major assumptions are that all surface radiation properties are grey and all radiation is diffuse. Both assumptions are reasonable for building zone interchange.

The ScriptF coefficients are developed by starting with the traditional direct radiation view factors. In the case of building rooms and zones, there are several complicating factors in finding the direct view factors—the main one being that the location of surfaces such as thermal mass representing furniture and partitions are not known. The other limitation is that the exact calculation of direct view factors is computationally very intensive even if the positions of all surfaces are known. Accordingly, EnergyPlus uses a procedure to approximate the direct view factors. The procedure has two steps:

1) Determine the total area of other surfaces “seen” by a surface.

2) Approximate the direct view factor from surface 1 to surface 2 as the ratio of the area of surface 2 to the total area “seen” by surface 1.

The determination of the “seen” area has several constraints:

• No surface sees itself.

• All surfaces see thermal mass surfaces.

• No surface facing within 10 degrees of another surface is seen by the other surface.

• All surfaces see roofs, floors and ceilings (subject to the preceding facing direction constraint).

Because the approximate view factors may not satisfy the basic requirements of reciprocity (two surfaces should exchange equal amounts of heat in each direction), and completeness (every surface should have a direct view factor sum of 1.0), EnergyPlus does a view factor fix operation before they are used in the ScriptF determination. Normally both of the requirements are satisfied, but in some special situations they are not, and special rules are applied.

• If a user includes less than four surfaces in a zone, the enforcement of “completeness” is not strictly possible because it is not possible to have a completed enclosure with less than four surfaces. In this situation, reciprocity is enforced initially and then view factors from each surface are summed. When dealing with this few surfaces, the result of using approximate view factors is that once reciprocity is enforced that the largest surface in the zone will end up having view factors that sum up to greater than unity. This is not physically possible. So, to correct this and maintain completeness for one surface, the entire view factor matrix is divided by the largest summation of view factors for any surface in the zone. By dividing the entire view factor matrix by the same value, the reciprocity established by the standard algorithm is maintained and completeness for one surface is established. All of the other surfaces will not have view factors that sum up to unity and thus will not satisfy completeness. However, these other surfaces will not have view factors that sum up to greater than unity thus avoiding a physically impossible situation.

• When there are four or more surfaces in a zone but the area of one surface in a zone is greater than roughly (0.9) the sum of the areas of all other surfaces, reciprocity only is enforced. Sometimes, for these very large surfaces, that enforcement of completeness without resorting to the very large surface seeing itself becomes impossible. In these rare cases, the large surface is allowed to see itself. Once this adjustment has been made, the calculation of the Script F values can proceed as normal via iteration until both completeness and reciprocity are satisfied. Note that this adjustment might happen for very “flat” zones where one or two surfaces are have a much larger area than the other surfaces like a Trombe Walls, a large single-story interior areas that are modeled as a single zone, or a space where a very large internal mass element has been added. Thermal mass surfaces also participate in internal long-wave radiation exchange and thus must have view factors calculated.

Warning messages are produced for both of these cases, and the results should be examined very carefully to ascertain that they are reasonable. The suggested action for the second case (the extra-large surface) is to divide the large surface into several smaller surfaces; then the enclosure will be treated as normal.

Once the ScriptF coefficients are determined, the longwave radiant exchange is calculated for each surface using:

$$q_{i,j}=A_i{F}_{i,j}\left(T_i^4-T_j^4\right)$$

where F${}_{i,j}$ is the ScriptF between surfaces i and j.

CarrollMRT[LINK]

The Carroll method is an approximation of gray-body long-wave radiation exchange within an enclosure that simplifies the surface-to-surface radiation exchange by using a single, mean radiant temperature node, $Tr$, that acts as a clearinghouse for the radiation heat exchange between surfaces. Instead of solving a dense-matrix, linear algebra problem, the mean radiant temperature can be calculated using a single equation, and subsequently used to determine the net long-wave radiation to/from each surface. Unlike the $O\left(n^2\right)$ complexity of the current dense-matrix solution, this approach has linear complexity.

The mean radiant temperature is calculated using three steps:

1. Calculation of the mean radiant temperature “view factor”, $Fi$. These view factors represent each surface’s “view” to the mean radiant temperature node as though all surfaces were part of a spherical enclosure (i.e., they all have equal view of the node regardless of their orientation to each other). $Fi$ is calculated as:

   $$F_i=\frac{1}{1-\frac{A_i{F}_i}{{\sum }_1^n{A}_j{F}_j}}$$

   Because of the circular reference in this equation, the collection of all “view factors” must be solved iteratively, but only once per simulation as surface areas do not change throughout. This converges for realistic enclosures but won’t necessarily converge for “enclosures” having only two or three surfaces, particularly if there are large area disparities.

2. Calculating the gray-body radiation resistance, $F^{\prime }i$. This calculation must be computed every time surface emissivity changes. $F^{\prime }i$ is calculated as:

   $$F_i^{\prime }=\frac{\sigma {\epsilon }_i}{\frac{{\epsilon }_i}{F_i}+1-{\epsilon }_i}$$

3. Finally, the mean radiant temperature, $Tr$, is:

   $$T_r=\frac{{\sum }_1^n{A}_i{F}_i^{\prime }{T}_i}{{\sum }_1^n{A}_i{F}_i^{\prime }}$$

Once the mean radiant temperature is known, the net radiation heat transfer for each surface can be calculated as:

$$q=F_i^{\prime }{A}_i(T_r^4-T_i^4)$$

Carroll, J. A., 1980, “An ‘MRT Method’ of Computing Radiant Energy Exchange in Rooms,” Proceedings of the Second Systems Simulation and Economic Analysis Conference, San Diego, CA.

Carroll, J. A., 1980a, “An MRT method of computing radiant energy exchange in rooms,” Proceedings of the 2nd Systems Simulation and Economic Analysis Conference, San Diego, CA.

Carroll, J. A., 1981, “A Comparison of Radiant Interchange Algorithms,” Proceedings of the 3rd Annual Systems Simulation and Economics Analysis/Solar Heating and Cooling Operational Results Conference, Reno. Solar Engineering, Proceedings of the ASME Solar division.

Thermal Mass and Furniture[LINK]

Furniture in a zone has the effect of increasing the amount of surface area that can participate in the radiation and convection heat exchanges. It also adds participating thermal mass to the zone. These two changes both affect the response to temperature changes in the zone and also affect the heat extraction characteristics.

The proper modeling of furniture is an area that needs further research, but the heat balance formulation allows the effect to be modeled in a realistic manner by including the furniture surface area and thermal mass in the heat exchange process.

LW Radiation From Internal Sources[LINK]

The traditional model for this source is to define a radiative/convective split for the heat introduced into a zone from equipment. The radiative part is then distributed over the surfaces within the zone in some prescribed manner. This, of course, is not a completely realistic model, and it departs from the heat balance principles. However, it is virtually impossible to treat this source in any more detail since the alternative would require knowledge of the placement and surface temperatures of all equipment.

Internal Short-Wave Radiation[LINK]

SW Radiation from Lights[LINK]

The short wavelength radiation from lights is distributed over the surfaces in the zone in some prescribed manner.

Transmitted Solar[LINK]

Transmitted solar radiation is also distributed over the surfaces in the zone in a prescribed manner. It would be possible to calculate the actual position of beam solar radiation, but that would involve partial surface irradiation, which is inconsistent with the rest of the zone model that assumes uniform conditions over an entire surface. The current procedures incorporate a set of prescribed distributions. Since the heat balance approach can deal with any distribution function, it is possible to change the distribution function if it seems appropriate.

Convection to Zone Air[LINK]

The convection flux is calculated using the heat transfer coefficients as follows:

$$q_{conv}^{\prime \prime }=h_c(T_s-T_a)$$

The inside convection coefficients (h${}_c$) can be calculated using one of many different models. Currently the implementation uses coefficients based on correlations for natural, mixed, and forced convection.

Interior Conduction[LINK]

This contribution to the inside surface heat balance is the wall conduction term, $q_{ki}^{\prime \prime }$ shown in Equation [eq:InsideFaceHeatBalanceEquation]. This represents the heat transfer to the inside face of the building element. Again, a CTF formulation is used to determine this heat flux.

Interior Convection[LINK]

There are many different modeling options available in EnergyPlus for inside convection coefficients, h${}_c$. There are four different settings to direct how EnergyPlus managers select h${}_c$ models during a simulation. There are numerous individual model equations for h${}_c$ in EnergyPlus to cover different situations that arise from surface orientations, room airflow conditions, and heat flow direction. Additionally, in many cases multiple researchers have developed competing models for the same situations that differ and there is no way to declare one is better than another. An overall default for the simulation is selected in the “SurfaceConvectionAlgorithm:Inside” object and can be overridden by selecting a different option in a zone description. These models are explained in the following sections. In addition to the correlation choices described below, it is also possible to override the convection coefficients on the inside of any surface by using the “SurfaceProperty:ConvectionCoefficients” object in the input file to set the convection coefficient value on the inside of any surface. The values can be specified directly or with schedules. Specific details are given in the Input Output Reference document.

For exterior simple-glazing windows modeled with the WindowMaterial:SimpleGlazingSystem object, h${}_c$ is scaled with an adjustment ratio. This enables the modeling of simple windows with highly conductive frames (large input U values). The calculation of the adjustment ratio is detailed in Section [application-issues].

Adaptive Convection Algorithm[LINK]

Beausoleil-Morrison (2000, 2002) developed a methodology for dynamically managing the selection of h${}_c$ equations called adaptive convection algorithm. The algorithm is used to select among the available h${}_c$ equations for the one that is most appropriate for a given surface at a given time. As Beausoleil-Morrison notes, the adaptive convection algorithm is intended to be expanded and altered to reflect different classification schemes and/or new h${}_c$ equations. The implementation in EnergyPlus has been modified from the original in the following ways:

• An input mechanism is provided (see the
  SurfaceConvectionAlgorithm:Inside:AdapativeModelSelections object) so that the user can customize the specific selections of h${}_c$ equations that are applied for different flow regimes and surface orientations. The changes apply in a general way to the entire model (but can be overridden by setting surface properties).

• To avoid requiring additional user input on the position of ZoneHVAC-type equipment within a zone, there is no distinction between zones that have convective zone heater equipment located underneath the windows and those that have convective heaters located away from the windows. This applies to the air flow regime associated with convective zone heaters. Using Beausoleil-Morrison’s terminology, regimes B1 and B2 are combined into just one B regime.

• To avoid requiring additional user input on the position of ZoneHVAC-type equipment within a zone, there is no distinction between surfaces that are directly blown on the fan and those that are away from the fan for the air flow regime associated with mechanical circulation from a zone fan (ZoneHVAC type equipment).

• The correlation for horizontal free jet developed by Fisher (1995) is not used. Ceiling diffuser models are used for all mechanical circulation from central air system. This decision was made for two reasons: (1) to avoid requiring additional user input on the position of, and momentum generated by, air terminal units, and (2) because Fisher (1995) found that the Coanda effect is so significant that in practice a free horizontal jet is difficult to maintain and mechanical-driven room airflows generally attach to surfaces and tend to match the flow regime of a ceiling diffuser much more often than a free jet.

• EnergyPlus supports arbitrary geometry so surfaces can be tilted with respect to vertical or horizontal. Beausoleil-Morrison’s adaptive convection algorithm was originally structured to use h${}_c$ equations that have no functional dependence on surface tilt angle. However, tilted surfaces do perform differently than vertical or horizontal surface when buoyancy forces are significant. Therefore, the EnergyPlus implementation expands the structure of the algorithm to include additional categories for tilted surfaces. The h${}_c$ equations developed by Walton (1983) are selected as the defaults for tilted surfaces because they have a functional dependence on tilt angle.

• Fohanno and Polidari (2006) produced a new h${}_c$ equation for vertical walls inside buildings with a simple buoyancy flow regime. They used a theoretical approach based on integral formalism and uniform heat flux (rather than uniform temperature) that covers both laminar and turbulent flow situations. In EnergyPlus, this model is selected as the default in place of the model by Alamdari and Hammond (1983) for vertical walls.

• Karadag (2009) produced a new h${}_c$ equation for ceiling surfaces that are actively chilled. He used computation fluid dynamics and various sized rooms and temperature conditions. In EnergyPlus, this model is selected as the default for surfaces that have active, in-ceiling cooling (in place of the model by Alamdari and Hammond (1983) for unstable ceilings).

• International Standard Organization (ISO) completed Standard 15099-2003 which includes h${}_c$ equations for the inside face of windows. EnergyPlus strives to adhere to formal modeling Standards where possible. Therefore the implementation includes a larger structure for the adaptive algorithm that includes additional categories for windows in all flow regimes and ISO 15099-2003 models are used as the default for windows in natural convection flow regimes. The ISO 15099 model applies to various tilt angles.

• Goldstein and Novosalec (2010) produced new h${}_c$ equations for forced air situations with ceiling slot diffusers along perimeters with significant glazing fractions. They used experiments with full-sized test room. These new equations are selected as the default for windows, ceilings and floors when there is an active central air system.

• Interior mass surfaces are assigned the h${}_c$ equation that would apply (stable or unstable) to a horizontal, upward facing surface for each flow regime.

• The algorithm switches between forced, mixed, and natural flow regimes by calculating the Richardson number, Ri = Gr/Re^2, for the zone. Large values of Ri indicate buoyancy dominates, while small values indicate forced flows dominate. To distinguish between opposing Zone unit type equipment (with fans) are assumed to force air up walls, and central air type equipment (with diffusers) are assumed to force air down walls.

The adaptive convection algorithm implemented in EnergyPlus for the inside face has a total of 45 different categories for surfaces and 29 different options for h${}_c$ equation selections. The following table summarizes the categories and the default assignments for h${}_c$ equations. The individual h${}_c$ equations are documented below.

*Indicates the default selection for h${}_c$ model equation.

Inside Face Surface Classification[LINK]

The adaptive convection algorithm is based on classifying surfaces by flow regime and orientation so that the correct h${}_c$ equation can be chosen at a particular point in time during the simulation. The classification depends on user input with some aspects processed only once at the beginning and others during each timestep. There are also various parameters or inputs to the h${}_c$ equations that need static or dynamic processing.

For each surface, it and the zone it is attached to are processed for the following static characteristics.

• Characteristic height for convection is taken as the zone height.

• Surfaces listed as receiving heat from Zone HVAC equipment with radiative models are considered “near” the heater.

• Zones are examined for low temperature radiant systems. The surfaces that contain the active elements are examined and the zone characterized to know if it has in-floor heating, in-ceiling cooling, or in-wall heating.

• A hydraulic diameter is calculated for horizontal surfaces for the entire zone.

and calculating various parameters needed by h${}_c$ equations. Selecting flow regime is done in the following manner. For each surface, we examine the zone on the inside face for the following:

• HVAC system type

• HVAC operating status

• HVAC system ACH

During initial zone and system sizing calculations, HVAC system characteristics are not known yet, so the flow regime is selected assuming an uncontrolled zone. During simulations, including HVAC sizing simulations, the HVAC system characteristics are used.

The surfaces are evaluated to determine:

• Surface classification: Floor, wall, roof, window, types

• Tilt angle

• Convective stability (sign of ΔT)

The individual h${}_c$ model equations and their respective references are listed in next by the keyword used to identify them.

ASHRAE Vertical Wall[LINK]

Walton adopted the following equation for natural convection from ASHRAE:

$$h=1.31{\left|\Delta T\right|}^{\frac{1}{3}}$$

This is also a component of the TARP overall algorithm described below.

Walton Unstable Horizontal Or Tilt[LINK]

Walton (1983) developed the following equation by fitting curves from various sources:

$$h=\frac{9.482{\left|\Delta T\right|}^{\frac{1}{3}}}{7.238-\left|\cos \Sigma \right|}$$

where $\Sigma$ is the surface tilt angle. Unstable refers to the direction of heat flow and the associated buoyancy relative to the surfaces. Unstable is when the natural tendency is to enhance flow in the sense that rising warmer air, or falling cooler air, is free to move away from the surface. This is also a component of the TARP overall algorithm described below.

Walton Stable Horizontal Or Tilt[LINK]

Walton (1983) developed the following equation by fitting curves from various sources:

$$h=\frac{1.810{\left|\Delta T\right|}^{\frac{1}{3}}}{1.382+\left|\cos \Sigma \right|}$$

where $\Sigma$ is the surface tilt angle. Stable refers to the direction of heat flow and the associated buoyancy relative to the surfaces. Stable is when the natural tendency is to retard flow in the sense that rising warmer air, or falling cooler air, is driven against the surface. This is also a component of the TARP overall algorithm described below.

Fisher Pedersen Ceiling Diffuser Walls[LINK]

Fisher and Pedersen (1997) developed the following equation from laboratory chamber measurements.

$$h=1.208+1.012\ast AC{H}^{0.604}$$

This is a component of the CeilingDiffuser overall algorithm described below.

Fisher Pedersen Ceiling Diffuser Ceiling[LINK]

Fisher and Pedersen (1997) developed the following equation from laboratory chamber measurements.

$$h=2.234+4.099\ast AC{H}^{0.503}$$

This is a component of the CeilingDiffuser overall algorithm described below.

Fisher Pedersen Ceiling Diffuser Floor[LINK]

Fisher and Pedersen (1997) developed the following equation from laboratory chamber measurements.

$$h=3.873+0.082\ast AC{H}^{0.98}$$

This is a component of the CeilingDiffuser overall algorithm described below.

Alamdari Hammond Stable Horizontal[LINK]

Alamdari and Hammond (1983) developed the following correlation for horizontal surfaces in stable thermal situation.

$$h=0.6{\left(\frac{\left|\Delta T\right|}{D_h^2}\right)}^{\text{1}/15\text{5}}$$

where,

$D_h=\frac{4A}{P}$ , hydraulic diameter of horizontal surface, A is area (m${}^2$) and P is the perimeter (m) of the entire zone.

Alamdari Hammond Unstable Horizontal[LINK]

Alamdari and Hammond (1983) developed the following correlation for horizontal surfaces in a buoyant thermal situation.

$$h={\left\{{\left[1.4{\left(\frac{\left|\Delta T\right|}{D_h}\right)}^{{}^{\text{1}/14\text{4}}}\right]}^6+{\left[1.63{\left|\Delta T\right|}^{\text{1}/13\text{3}}\right]}^6\right\}}^{\text{1}/16\text{6}}$$

Alamdari Hammond Vertical Wall[LINK]

Alamdari and Hammond (1983) developed the following correlation for vertical surfaces.

$$h={\left\{{\left[1.5{\left(\frac{\left|\Delta T\right|}{H}\right)}^{{}^{\text{1}/14\text{4}}}\right]}^6+{\left[1.23{\left|\Delta T\right|}^{\text{1}/13\text{3}}\right]}^6\right\}}^{\text{1}/16\text{6}}$$

where,

H is the characteristic height for the surface. In EnergyPlus this is the zone’s ceiling height (which could be larger than the height of an individual surface when wall are subdivided into more than one surface).

Khalifa Eq3 Wall Away From Heat[LINK]

Khalifa (1989) conducted experiments with test chambers and developed correlations for certain types of surfaces. One of them, identified as “Equation 3” in original reference, is for convectively heated zones and applies to the inside surfaces of walls away from the heat source:

$$h=2.07{\left|\Delta T\right|}^{0.23}$$

Khalifa Eq4 Ceiling Away From Heat[LINK]

Khalifa (1989) conducted experiments with test chambers and developed correlations for certain types of surfaces. One of them, identified as “Equation 4” in original reference, is for convectively heated zones and applies to the inside surfaces of ceilings away from the heat source:

$$h=2.72{\left|\Delta T\right|}^{0.13}$$

Khalifa Eq5 Wall Near Heat[LINK]

Khalifa (1989) conducted experiments with test chambers and developed correlations for certain types of surfaces. One of them, identified as “Equation 5” in original reference, is for convectively heated zones and applies to the inside surfaces of walls near the heat source:

$$h=1.98{\left|\Delta T\right|}^{0.32}$$

Khalifa Eq6 Non Heated Walls[LINK]

Khalifa (1989) conducted experiments with test chambers and developed correlations for certain types of surfaces. One of them, identified as “Equation 6” in original reference, is for heated zones and applies to the inside surfaces of walls that are not heated:

$$h=2.30{\left|\Delta T\right|}^{0.24}$$

Khalifa Eq7 Ceiling[LINK]

Khalifa (1989) conducted experiments with test chambers and developed correlations for certain types of surfaces. One of them, identified as “Equation 7” in original reference, is for heated zones and applies to the inside surfaces of ceilings:

$$h=3.10{\left|\Delta T\right|}^{0.17}$$

Awbi Hatton Heated Floor[LINK]

Awbi and Hatton (1999) conducted laboratory measurements using environmental chambers and developed the following correlation for floor surfaces that are being actively heated.

$$h=\frac{2.175{\left|\Delta T\right|}^{0.308}}{D_h^{0.076}}$$

where,

$D_h=\frac{4A}{P}$, hydraulic diameter of horizontal surface, A is area (m${}^2$) and P is the perimeter (m) of the entire zone (all of the adjacent floor surfaces if more than one in the zone).

Awbi Hatton Heated Wall[LINK]

Awbi and Hatton (1999) developed the following correlation for wall surfaces that are being actively heated.

$$h=\frac{1.823{\left|\Delta T\right|}^{0.293}}{D_h^{0.076}}$$

where,

$D_h=\frac{4A}{P}$, hydraulic diameter of wall surface, A is area (m${}^2$) and P is the perimeter (m) of the entire wall (all of the adjacent wall surfaces if more than one along the wall).

Beausoleil Morrison Mixed Assisted Wall[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for walls where the flow driving forces from mechanical forces are augmented by the driving forces from buoyancy.

$$h={({\{{\left[1.5{\left(\frac{\left|\Delta T\right|}{H}\right)}^{1/4}\right]}^6+{\left[1.23{\left|\Delta T\right|}^2\right]}^{1/6}\}}^{1/2}+{\left\{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right][-0.199+0.190\cdot AC{H}^{0.8}]\right\}}^3)}^{1/3}$$

where,

T${}_{SAT}$ is the supply air temperature at the diffuser.

Here the reference temperature is the zone air temperature rather than the diffuser supply air temperature.

Beausoleil Morrison Mixed Opposing Wall[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for walls where the flow driving forces from mechanical forces are opposed by the driving forces from buoyancy.

$$h=max\{\begin{array}{c}{({\{{\left[1.5{\left(\frac{\left|\Delta T\right|}{H}\right)}^{1/4}\right]}^6+{\left[1.23{\left|\Delta T\right|}^2\right]}^{1/6}\}}^{1/2}-{\{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot [-0.199+0.190\cdot AC{H}^{0.8}]\}}^3)}^{1/3}\\ 0.8\cdot {\{{\left[1.5{\left(\frac{\left|\Delta T\right|}{H}\right)}^{1/4}\right]}^6+\left[1.23{\left|\Delta T\right|}^2\right]\}}^{1/6}\\ 0.8\cdot \{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot [-0.199+0.190\cdot AC{H}^{0.8}]\}\end{array}$$

Beausoleil Morrison Mixed Stable Floor[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for floors where the flow driving forces include both mechanical forces and thermally stable buoyancy.

$$h={\left({\left\{0.6\cdot {\left(\frac{\left|\Delta T\right|}{D_H}\right)}^{\text{1}/15\text{5}}\right\}}^3+{\left\{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot \left[0.159+0.116AC{H}^{0.8}\right]\right\}}^3\right)}^{\text{1}/13\text{3}}$$

Beausoleil Morrison Mixed Unstable Floor[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for floors where the flow driving forces include both mechanical forces and thermally unstable buoyancy.

$$h={\left({\left\{{\left[1.4{\left(\frac{\left|\Delta T\right|}{D_h}\right)}^{{}^{\text{1}/14\text{4}}}\right]}^6+{\left[1.63{\left|\Delta T\right|}^{\text{1}/13\text{3}}\right]}^6\right\}}^{\text{3}/36\text{6}}+{\left\{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot \left[0.159+0.116AC{H}^{0.8}\right]\right\}}^3\right)}^{\text{1}/13\text{3}}$$

Beausoleil Morrison Mixed Stable Ceiling[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for ceilings where the flow driving forces include both mechanical forces and thermally Stable buoyancy.

$$h={\left({\left\{0.6\cdot {\left(\frac{\left|\Delta T\right|}{D_H}\right)}^{\text{1}/15\text{5}}\right\}}^3+{\left\{\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot \left[-0.166+0.484AC{H}^{0.8}\right]\right\}}^3\right)}^{\text{1}/13\text{3}}$$

Beausoleil Morrison Mixed Unstable Ceiling[LINK]

Beausoleil-Morrison (2000) used blending techniques to combine correlations originally developed by Alamdari and Hammond (1983) and Fisher and Pedersen (1997) to create the following correlation is for ceilings where the flow driving forces include both mechanical forces and thermally unstable buoyancy.

$$h={({({\left[1.4{\left(\frac{\left|\Delta T\right|}{D_h}\right)}^{1/4}\right]}^6+{\left[1.63{\left|\Delta T\right|}^{1/3}\right]}^6)}^{3/6}+{(\left[\frac{T_{surf}-T_{SAT}}{\left|\Delta T\right|}\right]\cdot [-0.166+0.484\cdot AC{H}^{0.8}])}^3)}^{1/3}$$

Fohanno Polidori Vertical Wall[LINK]

Fohanno and Polidori (2006) developed the following equation for h${}_c$ for vertical walls under simple buoyancy flow conditions.

$$h=\{\begin{array}{c}1.332{\left(\frac{\left|\Delta T\right|}{H}\right)}^{\text{1}/14\text{4}},R{a}_H^{\ast }\le 6.3\times {10}^9\\ 1.235e^{\left(0.0467H\right)}{\left|\Delta T\right|}^{0.316},R{a}_H^{\ast }>6.3\times {10}^9\end{array}$$

where,

$$R{a}_H^{\ast }=\frac{g{\beta }_f{q^{\prime \prime }}_c{H}^4}{k_f{\nu }_f^2}\underset{f}{Pr}$$

Karadag Chilled Ceiling[LINK]

Karadag (2009) used numerical methods to develop the following equation for h${}_c$ for chilled-ceiling surfaces.

$$h=3.1{\left|\Delta T\right|}^{0.22}$$

ISO 15099 Windows[LINK]

ISO Standard 15099-2003 includes the equations for h${}_c$ for room-side of windows and surfaces with any tilt angle and heat flow direction. The ISO 15099 correlation is for still room air angle and is determined in terms of the Nusselt number, $Nu$ , where

$$h_i=Nu\left(\frac{\lambda }{H}\right)$$

where,

$\lambda$ is the thermal conductivity of air, and

$H$ is the height of the window.

The Rayleigh number based on height, $R{a}_H$ , is calculated using:

$$R{a}_H=\frac{{\rho }^2{H}^{3}g{c}_p\left|T_{surf,i}-T_{air}\right|}{T_{m,f}\mu \lambda }$$

where,

$\rho$ is the density of air

$g$ is the acceleration due to gravity,

$c_p$ is the specific heat of air,

$\mu$ is the dynamic viscosity of air, and

$T_{m,f}$ is the mean film temperature in Kelvin given by,

$$T_{m,f}=T_{air}+\frac{1}{4}\left(T_{surf,i}-T_{air}\right)$$

There are four cases for the Nusselt correlation that vary by the tilt angle in degrees, $\gamma$ , and are based on heating conditions. For cooling conditions (where $T_{surf,i}>T_{air}$ ) the tilt angle is complemented so that $\gamma =180-\gamma$

Case A. $0^o\le \gamma <{15}^o$

$$Nu=0.13R{a}_H^{\text{1}/13\text{3}}$$

Case B. ${15}^o\le \gamma \le {90}^o$

$$R{a}_{cv}=2.5\times {10}^5{\left(\frac{e^{0.72\gamma }}{\sin \lambda }\right)}^{\text{1}/15\text{5}}$$

$$Nu=0.56{\left(R{a}_H\sin \gamma \right)}^{\text{1}/14\text{4}};forR{a}_H\le R{a}_{CV}$$

$$Nu=0.13\left(R{a}_H^{\text{1}/13\text{3}}-R{a}_{CV}^{\text{1}/13\text{3}}\right)+0.56{\left(R{a}_{CV}\sin \gamma \right)}^{\text{1}/14\text{4}};R{a}_H>R{a}_{CV}$$

Case C. ${90}^o<\gamma \le {179}^o$

$$Nu=0.56{\left(R{a}_H\sin \gamma \right)}^{\text{1}/14\text{4}};{10}^5\le R{a}_H\sin \gamma <{10}^{11}$$

Case D. ${179}^o<\gamma \le {180}^o$

$$Nu=0.58R{a}_H^{\text{1}/15\text{5}};R{a}_H\le {10}^{11}$$

The material properties are evaluated at the mean film temperature. Standard EnergyPlus psychrometric functions are used for $\rho$ and $c_p$ . Thermal conductivity is calculated using:

$\lambda =2.873\times {10}^{-3}+7.76\times {10}^{-8}{T}_{m,f}$ .

Kinematic viscosity is calculated using:

$\mu =3.723\times {10}^{-6}+4.94\times {10}^{-8}{T}_{m,f}$ .

This correlation depends on the surface temperature of the room-side glazing surface and is therefore included inside the window heat balance iteration loop.

Goldstein Novoselac Ceiling Diffuser Window[LINK]

Goldstein and Novoselac (2010) used laboratory chamber measurements to develop convection correlations for perimeter zones with highly glazed spaces served by overhead slot-diffuser-based air systems. The following are for bare windows in such spaces.

For WWR<50% with window in upper part of wall:

$$h=0.117{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

For WWR<50% with window in lower part of wall

$$h=0.093{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

For WWR > 50%

$$h=0.103{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

Where,

WWR is the window to wall ratio.

L is the length of exterior wall with glazing in the zone.

$\dot{V}$ is the air system flow rate in m${}^3$/s.

Goldstein Novoselac Ceiling Diffuser Walls[LINK]

Goldstein and Novoselac (2010) used laboratory chamber measurements to develop convection correlations for perimeter zones with highly glazed spaces served by overhead slot-diffuser-based air systems. The following are for exterior walls in such spaces.

For walls located below a window

$$h=0.063{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

For walls located above a window

$$h=0.093{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

Goldstein Novoselac Ceiling Diffuser Floor[LINK]

Goldstein and Novoselac (2010) used laboratory chamber measurements to develop convection correlations for perimeter zones with highly glazed spaces served by overhead slot-diffuser-based air systems. The following is for floors in such spaces.

$$h=0.048{\left(\frac{\dot{V}}{L}\right)}^{0.8}$$

Separate from the above model structure, there are also other comprehensive algorithm structures which are described below.

TARP Algorithm[LINK]

The comprehensive natural convection model, accessed using the keyword “TARP,” correlates the convective heat transfer coefficient to the surface orientation and the difference between the surface and zone air temperatures (where $\Delta$T = Surface Temperature - Air Temperature). The algorithm is taken directly from Walton (1983). Walton derived his algorithm from ASHRAE literature which can be found in the ASHRAE Handbook (for example , the ASHRAE Handbook of Fundamentals 2001, Table 5 on p. 3.12) which gives equations for natural convection heat transfer coefficients in the turbulent range for large, vertical plates and for large, horizontal plates facing upward when heated (or downward when cooled). A note in the text also gives an approximation for large, horizontal plates facing downward when heated (or upward when cooled) recommending that it should be half of the facing upward value. Walton adds a curve fit as a function of the cosine of the tilt angle to provide intermediate values between vertical and horizontal. The curve fit values at the extremes match the ASHRAE values very well.

For no temperature difference OR a vertical surface, the following correlation is used:

$$h=1.31{\left|\Delta T\right|}^{\frac{1}{3}}$$

It should be noted that h is zero when the temperature difference is zero.

For ($\Delta$T < 0.0 AND an upward facing surface) OR ($\Delta$T > 0.0 AND an downward facing surface), an enhanced convection correlation is used:

$$h=\frac{9.482{\left|\Delta T\right|}^{\frac{1}{3}}}{7.238-\left|\cos \Sigma \right|}$$

where $\Sigma$ is the surface tilt angle.

For ($\Delta$T > 0.0 AND an upward facing surface) OR ($\Delta$T < 0.0 AND an downward facing surface)], a reduced convection correlation is used:

$$h=\frac{1.810{\left|\Delta T\right|}^{\frac{1}{3}}}{1.382+\left|\cos \Sigma \right|}$$

where $\Sigma$ is the surface tilt angle.

Simple Natural Convection Algorithm[LINK]

The simple convection model uses constant coefficients for different heat transfer configurations, using the same criteria as the detailed model to determine reduced and enhanced convection. The coefficients are also taken directly from Walton (1983). Walton derived his coefficients from the surface conductances for e = 0.90 found in the ASHRAE Handbook (1985) in the table found on p. 23.2. The radiative heat transfer component was estimated at 1.02 * 0.9 = 0.918 BTU/h-ft2-F and then subtracted off. Finally the coefficients were converted to SI units to yield the values below.

For a vertical surface:

$$h=3.076$$

For a horizontal surface with reduced convection:

$$h=0.948$$

For a horizontal surface with enhanced convection:

$$h=4.040$$

For a tilted surface with reduced convection:

$$h=2.281$$

For a tilted surface with enhanced convection:

$$h=3.870$$

Ceiling Diffuser Algorithm[LINK]

The ceiling diffuser algorithm is based on empirical correlations developed by Fisher and Pedersen (1997). The correlation was reformulated to use the room outlet temperature as the reference temperature. The correlations are shown below.

For Floors:

$$h=3.873+0.082\ast AC{H}^{0.98}$$ (93)

The correlation for floors is illustrated in Figure 2.

For ceilings:

$$h=2.234+4.099\ast AC{H}^{0.503}$$

The correlation for ceilings is illustrated in Figure 3.

For Walls:

$$h=1.208+1.012\ast AC{H}^{0.604}$$

The correlation for walls is illustrated in Figure 4.

It should be noted that the conditions under which the correlations were developed included ACH values down to around 3.0. So, within EnergyPlus, for ACH values of 3.0 or greater, the ceiling diffuser correlations as shown above are used. Below ACH values of 0.5, the natural convection is determined using the TARP Algorithm (see above). Between ACH values of 0.5 and 3.0, convection is assumed to be a mix between natural and forced convection. Thus, between these two ACH values, EnergyPlus will use a linear interpolation between the natural convection algorithm value and the forced convection value at ACH=3.0 to avoid any discontinuities in the convection coefficient calculated and also maintain stability within the solution.

Trombe Wall Algorithm[LINK]

The Trombe wall algorithm is used to model convection in a “Trombe wall zone”, i.e. the air space between the storage wall surface and the exterior glazing. (See the later sections on Passive and Active Trombe Walls below for more information about Trombe walls.) The algorithm is identical to the convection model (based on ISO 15099) used in Window5 for convection between glazing layers in multi-pane window systems. The use of the algorithm for modeling an unvented Trombe wall has been validated against experimental data by Ellis (2003).

This algorithm gives the convection coefficients for air in a narrow vertical cavity that is sealed and not ventilated. This applies both to the air gap in between panes of a window or to the air gap between the Trombe wall glazing and the inner surface (often a selective surface). These convection coefficients are really the only difference between a normal zone and a Trombe zone. The rest of the zone heat balance is the same, e.g., transmitted solar, long-wave radiation between surfaces, etc.

For a vertical cavity, the correlation from ISO 15099 is:

$$\begin{array}{ccc}NU1=0.0673838R{a}^{\frac{1}{3}}& for& 5E4<Ra<1E6\\ NU1=0.028154R{a}^{0.4134}& for& 1E4<Ra<5E4\\ NU1=1+1.7596678E^{-10}R{a}^{2.2984755}& for& Ra\le 1E4\end{array}$$

$$NU2=0.242{\left(\text{{{rm{Ra}}}}/RaA\text{{rm{A}}}\right)}^{0.272}$$

$$NU=MAX(NU1,NU2)$$

where

Nu = Nusselt number

Ra = Rayleigh number

A = aspect ratio of cavity

This is then used in EnergyPlus as follows:

Net convection coefficient from glazing to wall is:

$$h_{net}=k\left(\text{{{rm{NU}}}}/NUL\text{{rm{L}}}\right)$$

where

k = conductivity of air

L = air gap thickness

Convection coefficient applied to each wall separately and actually used in the zone heat balance is:

$$h_c=2h_{net}$$

ASTM C1340 Algorithm[LINK]

The algorithm is based on ASTM C1340 Standard. This algorithm specifies to simulate heat flows in attic zone. The correlations account for the effects of surface-to-air temperature difference, heat flow direction, surface length along the heat flow direction, surface tilt angle, air speed past surface, and surface orientation. The correlations are in the form of a Nusselt number, Nu, as a function of a Rayleigh number, Ra, Grashof number, Gr, or a Reynolds number, Re.

$$Nu=hL/k$$ $$Ra=g\beta \rho C_p\Delta T{L}^3/\nu k$$ $$Gr=Ra/Pr$$ $$Pr=\nu /\alpha$$ $$Re=VL/\nu$$

h = convection heat transfer coefficient
L = characteristic length of surface
k = thermal conductivity of air
g = acceleration of gravity
$\beta$ = volume coefficient of expansion of air
$\rho$ = density of air
$C_p$ = specific heat of air
$\Delta$T = temperature difference between surface and air
$\nu$ = kinematic viscosity of air
Pr = Prandtl number for air
$\alpha$ = thermal diffusivity of air
V = velocity of air stream
The correlations used in the model are as below. (reference: Table 1 of ASTM C1340 standard)

1. Natural Convection:

   1. Horizontal surface, heat flow up
      $N{u}_n=0.54R{a}^{1/4}forRa<8\ast {10}^6$
      $N{u}_n=0.15R{a}^{1/3}forRa>8\ast {10}^6$

   2. Horizontal surface, heat flow down
      $N{u}_n=0.58R{a}^{0.2}$

   3. Vertical surface
      $N{u}_n=0.59R{a}^{1/4}forRa<1\ast {10}^9$
      $N{u}_n=0.10R{a}^{1/3}forRa>1\ast {10}^9$

   4. Nearly horizontal surface (tilt angle less than $2^o$), heat flow down
      $N{u}_n=0.58R{a}^{0.2}$

   5. Tilted surfaces (greater than $2^o$ tilt), heat flow down
      $N{u}_n=0.56\left(Rasin\right(\Phi ){)}^{1/4}\Phi =tiltangle$

   6. Tilted surfaces, heat flow up
      $N{u}_n=0.56\left(Rasin\right(\Phi ){)}^{1/4}forRa/Pr<G{r}_c$
      $N{u}_n=0.14(R{a}^{1/3}-(G{r}_{c}Pr{)}^{1/3})+0.56(G{r}_{c}Prsin\left(\Phi \right){)}^{1/4}forRa/Pr>G{r}_c$
      $G{r}_c=1\ast {10}^{6}for\Phi <{15}^o$
      $G{r}_c={10}^{\left(\Phi /\right(1.1870+0.0870\ast \Phi \left)\right)}for{15}^o<\Phi <{75}^o$
      $G{r}_c=5\ast {10}^{9}for\Phi >{75}^o$

2. Forced Convection:
   $N{u}_f=0.664P{r}^{1/3}R{e}^{1/2}forRe<5\ast {10}^5$
   $N{u}_f=P{r}^{1/3}(0.037R{e}^{0.8}-850)forRe>5\ast {10}^5$

Separate coefficients are calculated for natural and forced flow, and a mixed coefficient is calculated by taking the third root of the sum of the cubes of the two separate coefficients.

$$h_n=N{u}_n\ast k/L$$ $$h_f=N{u}_f\ast k/L$$ $$h_c={(h_n^3+h_f^3)}^{1/3}$$

More information regarding the implemented equations in the code is provided in detail in Appendix A, section 7.7 and 7.8 of the reference: Fontanini, A. D., Aguilar, J. L. C., Mitchell, M. S., Kosny, J., Merket, N., DeGraw, J. W., and Lee, E. (2018). Predicting the performance of radiant technologies in attics: Reducing the discrepancies between attic specific and whole-building energy models. Energy and Buildings, 169, 69-83.

Additional Heat Balance Source[LINK]

EnergyPlus also enables importing the pre-calculated results of other heat transfer processes to be considered in the inside heat balance calculation. An additional heat source term defined as a surface property would enable the consideration of these processes to be imported as schedules in the interior surface heat balance calculation in EnergyPlus.

The heat balance on the inside face can be written as follows:

$$q_{LWX}^{\prime \prime }+q_{SW}^{\prime \prime }+q_{LWS}^{\prime \prime }+q_{ki}^{\prime \prime }+q_{sol}^{\prime \prime }+q_{conv}^{\prime \prime }+q_{add}^{\prime \prime }=0$$

where: $q_{add}^{\prime \prime }$ = Pre-calculated results of the heat flux to the inside face from other heat transfer processes.

References[LINK]

Alamdari, F. and G.P. Hammond. 1983. Improved data correlations for buoyancy-driven convection in rooms. Building Services Engineering Research & Technology. Vol. 4, No. 3.

ASHRAE. 1985. 1985 ASHRAE Handbook – Fundamentals, Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.

ASHRAE. 2001. 2001 ASHRAE Handbook – Fundamentals, Atlanta: American Society of Heating, Refrigerating, and Air-Conditioning Engineers, Inc.

ASTM International. 2015. ASTM C1340 Standard Practice for Estimation of Heat Gain or Loss Through Ceilings Under Attics Containing Radiant Barriers by Use of a Computer Program

Awbi, H.B. and A. Hatton. 1999. Natural convection from heated room surfaces. Energy and Buildings 30 (1999) 233-244.

Beausoleil-Morrison, I. 2000. The adaptive coupling of heat and air flow modeling within dynamic whole-building simulations. PhD. Thesis. University of Strathclyde, Glasgow, UK.

Ellis, Peter G. 2003. Development and Validation of the Unvented Trombe Wall Model in EnergyPlus. Master’s Thesis, University of Illinois at Urbana-Champaign.

Fisher, D.E. and C.O. Pedersen. 1997. “Convective Heat Transfer in Building Energy and Thermal Load Calculations”, ASHRAE Transactions, Vol. 103, Pt. 2.

Fohanno, S., and G. Polidori. 2006. Modelling of natural convective heat transfer at an internal surface. Energy and Buildings 38 (2006) 548 - 553

Fontanini, A. D., Aguilar, J. L. C., Mitchell, M. S., Kosny, J., Merket, N., DeGraw, J. W., and Lee, E. (2018). Predicting the performance of radiant technologies in attics: Reducing the discrepancies between attic specific and whole-building energy models. Energy and Buildings, 169, 69-83.

Goldstein, K. and A. Novoselac. 2010. Convective Heat Transfer in Rooms With Ceiling Slot Diffusers (RP-1416). HVAC&R Research Journal TBD

Karadag, R. 2009. New approach relevant to total heat transfer coefficient including the effect of radiation and convection at the ceiling in a cooled ceiling room. Applied Thermal Engineering 29 (2009) 1561-1565

Khalifa AJN. 1989. Heat transfer processes in buildings. Ph.D. Thesis, University of Wales College of Cardiff, Cardiff, UK.

ISO. 2003. ISO 15099:2003. Thermal performance of windows, doors, and shading devices – Detailed calculations. International Organization for Standardization.

Walton, G. N. 1983. Thermal Analysis Research Program Reference Manual. NBSSIR 83-2655. National Bureau of Standards (now NIST). This is documentation for “TARP.”