Radiant System Models[LINK]

Low Temperature Radiant System Model[LINK]

The input objects ZoneHVAC:LowTemperatureRadiant:ConstantFlow,
ZoneHVAC:LowTemperatureRadiant:VariableFlow, and
ZoneHVAC:LowTemperatureRadiant:Electric provide models for low temperature radiant heating and cooling systems that appear, on the surface, to be relatively simple systems. The system circulates hot or cold fluid through tubes embedded in a wall, ceiling, or floor or runs current through electric resistance wires embedded in a surface or a panel. Energy is thus either added to or removed from the space, and zone occupants are conditioned by both radiation exchange with the system and convection from the surrounding air that is also affected by the system. Unless specifically required for indoor air quality considerations, fans, ductwork, dampers, etc. are not needed.

Despite the relative simplicity of the low temperature radiant systems, the integration of such a system within an energy analysis program requires one to overcome several challenges. First, for systems with significant thermal mass, the conduction transfer function method for modeling transient conduction must be extended to include embedded heat sources or sinks. Second, one must integrate this formulation within an energy analysis program like EnergyPlus. Finally, one must overcome the fact that the radiant system is both a zone heat balance element and a conditioning system. Each of these issues will be addressed in the next several subsections.

One Dimensional Heat Transfer Through Multilayered Slabs[LINK]

One of the most important forms of heat transfer in energy analysis is heat conduction through building elements such as walls, floors, and roofs. While some thermally lightweight structures can be approximated by steady state heat conduction, a method that applies to all structures must account for the presence of thermal mass within the building elements. Transient one dimensional heat conduction through a homogeneous layer with constant thermal properties such as the one shown in Figure 1 is governed by the following equation:

$$\frac{{\partial }^{2}T}{\partial x^2}=\frac{1}{\alpha }\frac{\partial T}{\partial t}$$

where:

T is the temperature as a function of position and time

x is the position

t is the time

$\alpha$ is the thermal diffusivity of the layer material which is equal to $\frac{k}{\rho c_p}$

k is the thermal conductivity

$\rho$ is the density

$c_p$ is the specific heat.

This equation is typically coupled with Fourier’s law of conduction that relates the heat flux at any position and time to temperature as follows:

$$q^{\prime \prime }\left(x,t\right)=-k\frac{\partial T\left(x,t\right)}{\partial x}$$

While analytical solutions exist for the single homogeneous layer shown in Figure 1, the solution becomes extremely tedious for the multiple layered slab shown in Figure 2.

Time Series Solutions: Conduction Transfer Functions[LINK]

Equations [eq:1DPartialDiffEqinTt] and [eq:FouriersLawOfConduction] can be solved numerically in a variety of ways. As mentioned in the previous section, other models have used control theory and numerical methods such as finite difference and finite element. However, each of these methods have drawbacks which render them inappropriate for use within an energy analysis program which requires both accuracy and efficiency from the simulation.

Another possible modeling method is a time series solution. Several of the detailed energy analysis programs such as EnergyPlus use a time series solution to transient heat conduction. The most basic time series solution is the response factor equation which relates the flux at one surface of an element to an infinite series of temperature histories at both sides as shown by:

$$q_{i,t}^{\prime \prime }=\sum _{m=1}^{\infty }{X}_m{T}_{i,t-m+1}-\sum _{m=1}^{\infty }{Y}_m{T}_{o,t-m+1}$$

where q" is heat flux, T is temperature, i signifies the inside of the building element, o signifies the outside of the building element, and t represents the current time step.

While in most cases the terms in the series decay fairly rapidly, the infinite number of terms needed for an exact response factor solution makes it less than desirable. Fortunately, the similarity of higher order terms can be used to replace them with flux history terms. The new solution contains elements that are called conduction transfer functions (CTFs). The basic form of a conduction transfer function solution is shown by the following equation:

$$q_{i,t}^{\prime \prime }=\sum _{m=1}^M{X}_m{T}_{i,t-m+1}-\sum _{m=1}^M{Y}_m{T}_{o,t-m+1}+\sum _{m=1}^k{F}_m{q^{\prime \prime }}_{i,t-m}$$

where k is the order of the conduction transfer functions, M is a finite number defined by the order of the conduction transfer functions, and X, Y, and F are the conduction transfer functions. This equation states that the heat flux at the interior surface of any generic building element for which the assumption of one dimensional conduction heat transfer is valid is linearly related to the current and some of the previous temperatures at both the interior and exterior surface as well as some of the previous flux values at the interior surface. A similar equation holds for the heat flux at the exterior surface.

The final CTF solution form reveals why it is so elegant and powerful. With a single, relatively simple equation, the conduction heat transfer through an element can be calculated. The coefficients (CTFs) in the equation are constants that only need to be determined once. The only storage of data required is the CTFs themselves and a limited number of temperature and flux terms. The formulation is valid for any surface type and does not require the calculation or storage of element interior temperatures.

As the next several sections will detail, there are two main methods for calculating conduction transfer functions: the Laplace Transform method and the State Space method. Both methods are well suited for the main focus of this research, the extension of conduction transfer functions to include heat sources or sinks.

Laplace Transform Formulation[LINK]

The traditional method for calculating conduction transfer functions is described in detail by Hittle (1981). Beginning with the transient one dimensional heat conduction equation {Equation [eq:1DPartialDiffEqinTt]} and Fourier’s law of conduction {Equation [eq:FouriersLawOfConduction]}, the Laplace transform method is used to convert the governing equations into the s-domain for a single layer such as the one shown in Figure 1.

$$\frac{d^{2}T\left(x,s\right)}{d{x}^2}=\frac{s}{\alpha }T\left(x,s\right)$$

$$q^{\prime \prime }\left(x,s\right)=-k\frac{dT\left(x,s\right)}{dx}$$

The transformed equations are solved and then put in matrix form as shown below:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_2\left(s\right)\\ q_2\left(s\right)\end{array}\right]$$

where:

$T_1\left(s\right)$, $T_2\left(s\right)$, $q_1\left(s\right)$, and $q_2\left(s\right)$ are the temperature and flux terms in the Laplace domain

$$A_1\left(s\right)=\cosh \left({\ell }_1\sqrt{s/s{\alpha }_1{\alpha }_1}\right)$$

$$B_1\left(s\right)=\left(1/1k_1\sqrt{s/s{\alpha }_1{\alpha }_1}{k}_1\sqrt{s/s{\alpha }_1{\alpha }_1}\right)\sinh \left({\ell }_1\sqrt{s/s{\alpha }_1{\alpha }_1}\right)$$

$$C_1\left(s\right)=k_1\sqrt{s/s{\alpha }_1{\alpha }_1}\sinh \left({\ell }_1\sqrt{s/s{\alpha }_1{\alpha }_1}\right)$$

$$D_1\left(s\right)=\cosh \left({\ell }_1\sqrt{s/s{\alpha }_1{\alpha }_1}\right)$$

$k_1$ is the thermal conductivity of the layer

${\alpha }_1$ is the thermal diffusivity of the layer

${\ell }_1$ is the thickness of the layer.

The 2 x 2 matrix consisting of $A_1\left(s\right)$, $B_1\left(s\right)$, $C_1\left(s\right)$, and $D_1\left(s\right)$ is called the transmission matrix and contains all of the thermophysical properties of the layer necessary to calculate transient conduction heat transfer through it. It can easily be shown that a second layer could be characterized in a similar way as:

$$\left[\begin{array}{c}{T}_2\left(s\right)\\ q_2\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_3\left(s\right)\\ q_3\left(s\right)\end{array}\right]$$

where $A_2\left(s\right)$, $B_2\left(s\right)$, $C_2\left(s\right)$, and $D_2\left(s\right)$ are calculated using the properties of the second layer. This can be substituted into Equation to provide insight how the extension to multilayered slabs is achieved.

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_3\left(s\right)\\ q_3\left(s\right)\end{array}\right]$$

Thus, for a multilayered element as shown in Figure 2, each separate layer has a transmission matrix of $A_i\left(s\right)$, $B_i\left(s\right)$, $C_i\left(s\right)$, and $D_i\left(s\right)$ associated with it. The form of the matrix equation for the multilayered element is the same as the equation for a single layer:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}A\left(s\right)& B\left(s\right)\\ C\left(s\right)& D\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_{n+1}\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]$$

but the transmission matrix is replaced by:

$$\left[\begin{array}{cc}A\left(s\right)& B\left(s\right)\\ C\left(s\right)& D\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\cdots \left[\begin{array}{cc}{A}_n\left(s\right)& B_n\left(s\right)\\ C_n\left(s\right)& D_n\left(s\right)\end{array}\right]$$

Equation [eq:LaplaceMultiLayerMatrix] is typically rearranged as follows:

$$\left[\begin{array}{c}{q}_1\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]=\left[\begin{array}{cc}\frac{D\left(s\right)}{B\left(s\right)}& \frac{-1}{B\left(s\right)}\\ \frac{1}{B\left(s\right)}& \frac{-A\left(s\right)}{B\left(s\right)}\end{array}\right]\left[\begin{array}{c}{T}_1\left(s\right)\\ T_{n+1}\left(s\right)\end{array}\right]$$

which relates the flux at either surface of the element to the temperature histories at both surfaces. When the temperature histories are formulated as triangular pulses made up of simple ramp functions, the roots of this equation can be found and result in response factors. The response factors can be simplified as described above through the introduction of flux history terms to form conduction transfer functions. A simplified method of finding the roots of the Laplace domain equations is described by Hittle and Bishop (1983) and is used by the current version of BLAST.

State Space Formulation[LINK]

Recently, another method of finding conduction transfer functions starting from a state space representation has begun receiving increased attention (Ceylan and Myers 1980; Seem 1987; Ouyang and Haghighat 1991). The basic state space system is defined by the following linear matrix equations:

$$\frac{d\left[x\right]}{dt}=\left[A\right]\left[x\right]+\left[B\right]\left[u\right]$$

$$\left[y\right]=\left[C\right]\left[x\right]+\left[D\right]\left[u\right]$$

where x is a vector of state variables, u is a vector of inputs, y is the output vector, t is time, and A, B, C, and D are coefficient matrices. Through the use of matrix algebra, the vector of state variables (x) can be eliminated from the system of equations, and the output vector (y) can be related directly to the input vector (u) and time histories of the input and output vectors.

This formulation can be used to solve the transient heat conduction equation by enforcing a finite difference grid over the various layers in the building element being analyzed. In this case, the state variables are the nodal temperatures, the environmental temperatures (interior and exterior) are the inputs, and the resulting heat fluxes at both surfaces are the outputs. Thus, the state space representation with finite difference variables would take the following form:

$$\frac{d\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]}{dt}=\left[A\right]\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]+\left[B\right]\left[\begin{array}{c}{T}_i\\ T_o\end{array}\right]$$

$$\left[\begin{array}{c}{q^{\prime \prime }}_i\\ {q^{\prime \prime }}_o\end{array}\right]=\left[C\right]\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]+\left[D\right]\left[\begin{array}{c}{T}_i\\ T_o\end{array}\right]$$

where $T_1$, $T_2$, …, $T_{n-1}$, Tn are the finite difference nodal temperatures, n is the number of nodes, Ti and To are the interior and exterior environmental temperatures, and $q_i^{\prime \prime }$ and $q_o^{\prime \prime }$ are the heat fluxes (desired output).

Seem (1987) shows that for a simple one layer slab with two interior nodes as in Figure 3 and convection at both sides the resulting finite difference equations are given by:

$$C\frac{d{T}_1}{dt}=hA\left(T_o-T_1\right)+\frac{T_2-T_1}{R}$$

$$C\frac{d{T}_2}{dt}=hA\left(T_i-T_2\right)+\frac{T_1-T_2}{R}$$

$$q_i^{\prime \prime }=h\left(T_i-T_2\right)$$

$$q_o^{\prime \prime }=h\left(T_1-T_o\right)$$

where:

$R$ is the thermal resistance which is equal to $\frac{\ell }{kA}$

$C$ is the thermal capacitance which is equal to $\frac{\rho c_p\ell A}{2}$

$T_o$ is the outside temperature

$T_i$ is the inside temperature

$T_1$ is the temperature of node 1

$T_2$ is the temperature of node 2

A is the area of the surface exposed to the environmental temperatures.

In matrix format:

$$\left[\begin{array}{c}\frac{d{T}_1}{dt}\\ \frac{d{T}_2}{dt}\end{array}\right]=\left[\begin{array}{cc}\frac{-1}{RC}-\frac{hA}{C}& \frac{1}{RC}\\ \frac{1}{RC}& \frac{-1}{RC}-\frac{hA}{C}\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]+\left[\begin{array}{cc}\frac{hA}{C}& 0\\ 0& \frac{hA}{C}\end{array}\right]\left[\begin{array}{c}{T}_o\\ T_i\end{array}\right]$$

$$\left[\begin{array}{c}{q^{\prime \prime }}_i\\ {q^{\prime \prime }}_o\end{array}\right]=\left[\begin{array}{cc}0& -h\\ h& 0\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]+\left[\begin{array}{cc}0& h\\ -h& 0\end{array}\right]\left[\begin{array}{c}{T}_o\\ T_i\end{array}\right]$$

The important aspect of the state space technique is that through the use of relatively simple matrix algebra the state space variables (nodal temperatures) can be eliminated to arrive at a matrix equation that gives the outputs (heat fluxes) as a function of the inputs (environmental temperatures) only. This eliminates the need to solve for roots in the Laplace domain. In addition, the resulting matrix form has more physical meaning than complex functions required by the Laplace transform method. The current version of EnergyPlus uses the state space method for computing CTFs.

The accuracy of the state space method of calculating CTFs has been addressed in the literature. Ceylan and Myers (1980) compared the response predicted by the state space method to various other solution techniques including an analytical solution. Their results showed that for an adequate number of nodes the state space method computed a heat flux at the surface of a simple one layer slab within 1% of the analytical solution. Ouyang and Haghighat (1991) made a direct comparison between the Laplace and state space methods. For a wall composed of insulation between two layers of concrete, they found almost no difference in the response factors calculated by each method.

Extension of Time Series Solutions to Include Heat Sources and Obtain Internal Temperatures[LINK]

Laplace Transform Formulation[LINK]

Degiovanni (1988) proposed two methodologies for including sources or sinks in the Laplace Transform Formulation. The first method shows how a source that varies as a function of time and location can be incorporated. The resulting equations involve some fairly complicated terms including spatial derivatives.

The second method that will be analyzed in more detail involves the addition of a source or sink at the interface between two layers. The derivation of the necessary equations is begun by analyzing the simple two layer element shown in Figure 4.

For the first layer, it was determined that in the Laplace domain:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_2\left(s\right)\\ q_2\left(s\right)\end{array}\right]$$

For the second layer:

$$\left[\begin{array}{c}{T}_2\left(s\right)\\ q_2\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_3\left(s\right)\\ q_3\left(s\right)\end{array}\right]$$

To link the two layers and include the heat source between them, the following substitution is made:

$$\left[\begin{array}{c}{T}_2\left(s\right)\\ q_2\left(s\right)\end{array}\right]=\left[\begin{array}{c}{T}_{2+}\left(s\right)\\ q_{2+}\left(s\right)\end{array}\right]+\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]$$

which results in:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left\{\left[\begin{array}{c}{T}_{2+}\left(s\right)\\ q_{2+}\left(s\right)\end{array}\right]+\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]\right\}$$

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left\{\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_3\left(s\right)\\ q_3\left(s\right)\end{array}\right]+\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]\right\}$$

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{cc}{A}_2\left(s\right)& B_2\left(s\right)\\ C_2\left(s\right)& D_2\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_3\left(s\right)\\ q_3\left(s\right)\end{array}\right]+\left[\begin{array}{cc}{A}_1\left(s\right)& B_1\left(s\right)\\ C_1\left(s\right)& D_1\left(s\right)\end{array}\right]\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]$$

While Degiovanni concludes with this formula, some insight into what the generic equation for an element that has n layers might look like is gained by working with Equation [eq:TwoLayerEquivalentQTFLaplaceForm]. If a layer is added to the left of the first layer, the entire right hand side of Equation [eq:TwoLayerEquivalentQTFLaplaceForm] is multiplied by the transmission matrix of the new layer. Conversely, if a layer is added to the right of the second layer in Figure 4, the vector containing the Laplace transform of the temperature and heat flux at interface 3 is replaced by the product of the transmission matrix of the new layer and the vector for temperature and heat flux at the next interface, and the term dealing with the heat source is not affected. The general equation for a building element with n layers and m layers between the left hand surface and the heat source can be derived as:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left(\prod _{i=1}^n\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]\right)\left[\begin{array}{c}{T}_{n+1}\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]+\left(\prod _{i=1}^m\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]\right)\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]$$

or in more compact form:

$$\left[\begin{array}{c}{T}_1\left(s\right)\\ q_1\left(s\right)\end{array}\right]=\left[\begin{array}{cc}A\left(s\right)& B\left(s\right)\\ C\left(s\right)& D\left(s\right)\end{array}\right]\left[\begin{array}{c}{T}_{n+1}\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]+\left[\begin{array}{cc}a\left(s\right)& b\left(s\right)\\ c\left(s\right)& d\left(s\right)\end{array}\right]\left[\begin{array}{c}0\\ q_{source}\left(s\right)\end{array}\right]$$

where: $\left[\begin{array}{cc}A\left(s\right)& B\left(s\right)\\ C\left(s\right)& D\left(s\right)\end{array}\right]=\prod _{i=1}^n\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]$ and $\left[\begin{array}{cc}a\left(s\right)& b\left(s\right)\\ c\left(s\right)& d\left(s\right)\end{array}\right]=\prod _{i=1}^m\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]$ .

Next, Equation [eq:LaplaceTransformDerivation649] must be rearranged to match the form of Equation [eq:LaplaceTransformDeriv631], which relates the heat flux at both sides of the element to the temperature at each side. The matrix equation that is obtained shows that:

$$\left[\begin{array}{c}{q}_1\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]=\left[\begin{array}{cc}\frac{D\left(s\right)}{B\left(s\right)}& \frac{-1}{B\left(s\right)}\\ \frac{1}{B\left(s\right)}& \frac{-A\left(s\right)}{B\left(s\right)}\end{array}\right]\left[\begin{array}{c}{T}_1\left(s\right)\\ T_{n+1}\left(s\right)\end{array}\right]+\left[\begin{array}{c}d\left(s\right)-\frac{D\left(s\right)b\left(s\right)}{B\left(s\right)}\\ \frac{b\left(s\right)}{B\left(s\right)}\end{array}\right]\left[q_{source}\left(s\right)\right]$$

This equation bears a striking resemblance to Equation [eq:LaplaceTransformDeriv631]. If the source term in Equation [eq:LaplaceTransformDerivation650] is dropped, then the equation is identical to Equation [eq:LaplaceTransformDeriv631]. This result conforms with the superposition principle which was used to develop the conduction transfer functions from the summation of a series of triangular pulses or ramp sets. Now, the effect of the heat source is simply added to the response to the temperature inputs.

While Equation [eq:LaplaceTransformDerivation650] is correct for any single or multilayered element, the first term in the heat source transmission matrix does not appear to match the compactness of the other terms in the matrix equation. It can be shown (see Strand 1995: Equations 32 through 42 which detail this derivation) that the heat source transmission term for a two-layer problem reduces to:

$$\left[\begin{array}{c}{q}_1\left(s\right)\\ q_3\left(s\right)\end{array}\right]=\left[\begin{array}{cc}\frac{D\left(s\right)}{B\left(s\right)}& \frac{-1}{B\left(s\right)}\\ \frac{1}{B\left(s\right)}& \frac{-A\left(s\right)}{B\left(s\right)}\end{array}\right]\left[\begin{array}{c}{T}_1\left(s\right)\\ T_3\left(s\right)\end{array}\right]+\left[\begin{array}{c}\frac{B_2\left(s\right)}{B\left(s\right)}\\ \frac{B_1\left(s\right)}{B\left(s\right)}\end{array}\right]\left[q_{source}\left(s\right)\right]$$

If this is extended to a slab with n layers and a source between the m and m+1 layers, the general matrix equation for obtaining heat source transfer functions using the Laplace transform method is:

$$\left[\begin{array}{c}{q}_1\left(s\right)\\ q_{n+1}\left(s\right)\end{array}\right]=\left[\begin{array}{cc}\frac{D\left(s\right)}{B\left(s\right)}& \frac{-1}{B\left(s\right)}\\ \frac{1}{B\left(s\right)}& \frac{-A\left(s\right)}{B\left(s\right)}\end{array}\right]\left[\begin{array}{c}{T}_1\left(s\right)\\ T_{n+1}\left(s\right)\end{array}\right]+\left[\begin{array}{c}\frac{\overline{b}\left(s\right)}{B\left(s\right)}\\ \frac{b\left(s\right)}{B\left(s\right)}\end{array}\right]\left[q_{source}\left(s\right)\right]$$

where:

$$\left[\begin{array}{cc}A\left(s\right)& B\left(s\right)\\ C\left(s\right)& D\left(s\right)\end{array}\right]=\prod _{i=1}^n\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]$$

$$\left[\begin{array}{cc}a\left(s\right)& b\left(s\right)\\ c\left(s\right)& d\left(s\right)\end{array}\right]=\prod _{i=1}^m\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]$$

$$\left[\begin{array}{cc}\overline{a}\left(s\right)& \overline{b}\left(s\right)\\ \overline{c}\left(s\right)& \overline{d}\left(s\right)\end{array}\right]=\prod _{i=m+1}^n\left[\begin{array}{cc}{A}_i\left(s\right)& B_i\left(s\right)\\ C_i\left(s\right)& D_i\left(s\right)\end{array}\right]$$

At first glance, the terms in the heat source transmission matrix may appear to be reversed. It is expected that only the layers to the left of the source will affect $q_1\left(s\right)$, but the presence of $\overline{b}\left(s\right)$ in the element multiplied by $q_{source}\left(s\right)$ to obtain $q_1\left(s\right)$ seems to be contradictory. In fact, the entire term, $\overline{b}\left(s\right)/\overline{b}\left(s\right)B\left(s\right)B\left(s\right)$ , must be analyzed to determine the effect of $q_{source}\left(s\right)$ on $q_1\left(s\right)$. In essence, the appearance of $\overline{b}\left(s\right)$ removes the effects of the layers to the right of the source from B(s) leaving only the influence of the layers to the left of the source. The form displayed by Equation [eq:LaplaceTransformDeriv652] is, however, extremely convenient because the terms in the heat source transmission matrix have the same denominators, and thus roots, as the terms in the temperature transmission matrix. Thus, the same roots that are calculated for the CTFs can be used for the QTFs, saving a considerable amount of computer time during the calculation of the transfer functions.

Once Equation [eq:LaplaceTransformDeriv652] is inverted from the Laplace domain back into the time domain, the combined CTF–QTF solution takes the following form:

$$q_{i,t}^{\prime \prime }=\sum _{m=1}^M{X}_m{T}_{i,t-m+1}-\sum _{m=1}^M{Y}_m{T}_{o,t-m+1}+\sum _{m=1}^k{F}_m{q^{\prime \prime }}_{i,t-m}+\sum _{m=1}^M{W}_m{q}_{source,t-m+1}$$

This relation is identical to Equation [eq:LaplaceTransformDeriv623] except for the presence of the QTF series that takes the heat source or sink into account.

State Space Formulation[LINK]

The two-node example introduced by Seem (1987) can be utilized to examine the extension of the state space method to include heat sources or sinks. Figure 5 shows the simple two node network with a heat source added at node 1.

The nodal equations for the finite difference network shown in Figure 5 are:

$$C\frac{d{T}_1}{dt}=hA\left(T_o-T_1\right)+\frac{T_2-T_1}{R}+q_{source}A$$

$$C\frac{d{T}_2}{dt}=hA\left(T_i-T_2\right)+\frac{T_1-T_2}{R}$$

$$q_i^{\prime \prime }=h\left(T_i-T_2\right)$$

$$q_o^{\prime \prime }=h\left(T_1-T_o\right)$$

In obtaining the matrix equivalent for this set of equations, it should be noted that the source term is not a constant but rather an input that varies with time. Thus, it must be grouped with the environmental temperatures as inputs. The resulting matrix equations take the following form:

$$\left[\begin{array}{c}\frac{d{T}_1}{dt}\\ \frac{d{T}_2}{dt}\end{array}\right]=\left[\begin{array}{cc}\frac{-1}{RC}-\frac{hA}{C}& \frac{1}{RC}\\ \frac{1}{RC}& \frac{-1}{RC}-\frac{hA}{C}\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]+\left[\begin{array}{ccc}\frac{hA}{C}& 0& \frac{A}{C}\\ 0& \frac{hA}{C}& 0\end{array}\right]\left[\begin{array}{c}{T}_o\\ T_i\\ q_{source}\end{array}\right]$$

$$\left[\begin{array}{c}{q^{\prime \prime }}_1\\ {q^{\prime \prime }}_2\end{array}\right]=\left[\begin{array}{cc}0& -h\\ h& 0\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]+\left[\begin{array}{ccc}0& h& 0\\ -h& 0& 0\end{array}\right]\left[\begin{array}{c}{T}_o\\ T_i\\ q_{source}\end{array}\right]$$

Equation [eq:LaplaceTransformDeriv659] appears to suggest that the source term has no direct effect on the heat flux at either side of the element because its coefficients are zero. This is not the case. Equation [eq:LaplaceTransformDeriv659] only relates variables that have a direct influence on heat flux. So, while $T_i$ has no direct influence on $q_o^{\prime \prime }$, it does have an indirect influence through the nodal network. The same would hold for the influence of $q_{source}$.

If this analysis is extended to a finite difference network with n nodes, the corresponding matrix equations can be shown to be:

$$\frac{d\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]}{dt}=\left[A\right]\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]+\left[B\right]\left[\begin{array}{c}{T}_o\\ T_i\\ q_{source}\end{array}\right]$$

$$\left[\begin{array}{c}{q^{\prime \prime }}_i\\ {q^{\prime \prime }}_o\end{array}\right]=\left[C\right]\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]+\left[D\right]\left[\begin{array}{c}{T}_o\\ T_i\\ q_{source}\end{array}\right]$$

The influence of the heat source is also confirmed by the final solution form, which is identical to the Laplace transform result shown in Equation [eq:LaplaceTransformDeriv653] As with the Laplace solution method, the state space method results in a set of QTFs that relate the heat source at the current time step and several previous time steps to the current heat flux at the surface of the element.

Other similarities between the two solution methods are evident. It is interesting to note that as with the Laplace method there is no alteration of the CTFs calculated by the state space method. Thus, the principle of superposition is still valid. Furthermore, the introduction of the source term did not substantially increase the computing effort required to calculate the additional transfer functions. In the Laplace method, this was shown by the common roots, B(s), shared by both the CTFs and the QTFs. In the state space method, it can be noted that the A matrices in Equations [eq:QTFDeriv640] and [eq:QTFDeriv658] are identical. Since the state space method requires the inversion and the exponentiation of the A matrix only, the additional QTF terms will not require a substantial amount of additional computing time for their calculation.

Two-Dimensional Solutions for Radiant Systems[LINK]

One distinct advantage of the State Space method presented in the previous section over the Laplace Transform method is that it can be adapted to more than one dimension. In fact, as long as one can apply a network of nodes to a problem, the State Space method can be adapted to it. For EnergyPlus, the biggest implication is that conduction through a construction can be expanded from one-dimensional in nature to two-dimensional. This is particularly important in a hydronic radiant system where the presence of water tubes is clearly more than one-dimensional.

The modeling of two-dimensional conduction used when ConstructionPoperty:InternalHeatSource objects are attached to surfaces is thus possible by expanding the network of nodes shown in the above examples in the direction perpendicular to the main direction of heat transfer. Essentially, when a user requests the two-dimensional solution for a ConstructioniProperty:InternalHeatSource surface, the network of nodes is expanded perpendicular to the direction of heat transfer. While this is not an overly complex process, here are some rules and limitations to the application of two-dimensional solutions in EnergyPlus:

• Even though the solution internal to the surface is two-dimensional, in order to work within the confines of the standard EnergyPlus heat balance, the boundary condition at both the inside and outside face of the surface is that the temperature across the surface is the same at all points or that the surfaces are still isothermal. The point of this is so that once the conduction transfer functions are calculated that the the surface heat balance formulations remain the same as any other surface that is using a one-dimensional assumption in EnergyPlus.

• The domain being modeled for a two-dimensional approach goes from one face of the surface to the other face as a one-dimensional model. In the second dimension, perpendicular to the main direction of heat transfer, the domain goes from the line that goes through the tubing to the line of symmetry at the mid-point between adjacent tubing. The perpendicular distance for this domain is given by the tubing spacing user input in the ConstructionProperty:InternalHeatSource object. This tubing spacing is halved to determine the perpendicular distance for the solution domain.

• While the number of nodes used for each layer of a construction is determined by the thermo-physical properties of the material for the layer to maintain solution stability, the number of nodes in the perpendicular direction is fixed for all layers of the construction. Currently in EnergyPlus, the number of nodes in the perpendicular direction is fixed at 7. This number was chosen as a result of testing with an evaluation version of one of the precedessor programs of EnergyPlus. This number of nodes was a balance between accuracy requirements and solution speed. Because the speed of the process required to calculate conduction transfer functions increases greatly as the number of nodes used in the model increases, increasing the number of nodes in the perpendicular direction too much will result in an unacceptible increase in the time required to calculate the conduction transfer functions.

• The heat source or sink is applied evenly over the entire node where the user defines the location of the source through the ConstructionProperty:InternalHeatSource input. For calculating the conduction transfer functions, the model ignores the presence of the tubing and fluid and simply assumes that the entire layer consists of the appropriate material as defined by the Construction:InternalSource input. The radiant model does take into account heat transfer between the material and the fluid being circulated through the system.

• The location of the heat source and the calculation of an additional point that can be used for controlling the slab (see next section) are defined by input provided in the ConstructionPropery:InternalHeatSource input. The location of the tubing where heat is added or subtracted from the slab is always defined at one side of the solution domain. The location of the additional user temperature request is controlled by both the fields for the location of the user temperature request (between two layers) and the perpendicular direction for this temperature. This temperature can then be used for controlling a radiant system. (See section on Low Temperature Radiant System Controls for more information about controlling a radiant system.

These assumptions are applied to a low temperature radiant system within EnergyPlus and their impact can be seen in the next three diagrams. First, one can see the symmetry inherent in a low temperature radiant system in Figure 6 below.

Given the lines of symmetry in Figure 6, the solution domain can be narrowed down to what is shown in Figure 7 below.

Once the solution domain has been set, a network of nodes can be applied to this solution domain for the State Space method using the assumptions given above. An example of this is shown in Figure 8 below. Note that in EnergyPlus there are actually 7 nodes in perpendicular direction, not 5 as shown in this figure.

As mentioned in the figure above, the source or sink is applied at a node along the left side of this solution domain and the depth established by user input. The number of nodes in the main direction of heat transfer for each layer of the construction is calculated in the same manner as for a standard layer with 1-D heat conduction. The number of nodes chosen for the perpendicular direction is set to 7 in EnergyPlus as a balance between calculation accuracy and solution speed.

Determination of Internal Temperatures[LINK]

One aspect of low temperature radiant systems that has not been addressed to this point is the appropriateness of specifying the effect of the system on slab response via a heat source term. For a heating system that employs electrical resistance heating, the use of a heat source as the input variable is logical. The heat produced by such a system can easily be related to the current passing through the heating wire. However, for a hydronic heating or cooling system, the known quantity is not heat but rather the temperature of the water being sent to the building element.

The use of a temperature to simulate the presence of a heating or cooling system presents one major obstacle. When fluid is not being circulated, there is no readily available temperature value available for use as an input variable.

In a hydronic system, a link between the fluid temperature being sent to the slab and the heat delivered to the slab exist. The most effective way of relating these two variables is to consider the slab to be a heat exchanger. Using heat exchanger relationships, an equation could then be formulated to obtain the heat delivered to the slab based on the inlet fluid temperature.

Most heat exchangers are used to thermally link two fluids. In the case of a hydronic radiant system, there is only one fluid and a stationary solid. Presumably, if the inlet fluid temperature, the system geometry, and the solid temperature are known, then the outlet temperature and thus the heat transfer to the building element can be computed. This leads to an interesting question: what is the solid temperature?

By definition, for one dimensional conduction heat transfer, the solid temperature is the temperature of the building element at the depth where the hydronic loop is located. Typically, this temperature is not known because it is not needed. The goal of both methods of calculating CTFs was the elimination of internal temperatures that were not needed for the simulation. For a hydronic system, it is necessary to extract this information to solve for the heat source term. Two methods of accomplishing this are described below.

Returning to the two layer example shown in Figure 4, it can be shown that the final solution form in the time domain for the slab with a source at the interface between the two layers is:

$$q_{1,t}^{\prime \prime }=\sum _{m=1}^M{X}_{k,m}{T}_{1,t-m+1}-\sum _{m=1}^M{Y}_{k,m}{T}_{3,t-m+1}+\sum _{m=1}^k{F}_m{q^{\prime \prime }}_{1,t-m}+\sum _{m=1}^M{W}_m{q}_{source,t-m+1}$$

A similar equation could be written for the response of the first layer in absence of any source term and is given by:

$$q_{1,t}^{\prime \prime }=\sum _{m=1}^M{x}_{k,m}{T}_{1,t-m+1}-\sum _{m=1}^M{y}_{k,m}{T}_{2,t-m+1}+\sum _{m=1}^k{f}_m{q^{\prime \prime }}_{1,t-m}$$

While the current temperature at the interface is not known, presumably the previous values of this parameter will be known. In addition, the temperatures and the flux histories at surface 1 are also know. The unknowns in Equation [eq:QTFDeriv663] are the current heat flux at surface 1 and the temperature at surface 2. However, Equation [eq:QTFDeriv662] does define the current value of the heat flux at surface 1 based on temperature, heat flux, and heat source histories. Thus, if this value is used in Equation [eq:QTFDeriv663], the only remaining unknown in this equation is the current temperature at surface 2, the surface where the heat source or sink is present. Rearranging Equation [eq:QTFDeriv663] provides an equation from which the temperature at the source location may be calculated:

$$T_{2,t}=\sum _{m=1}^M{\overline{X}}_{k,m}{T}_{1,t-m+1}-\sum _{m=1}^{M-1}{\overline{Y}}_{k,m}{T}_{2,t-m}+\sum _{m=1}^{k+1}{\overline{F}}_m{q^{\prime \prime }}_{1,t-m+1}$$

where the new coefficients are obtained from the standard conduction transfer functions for the first layer via the following equations:

$${\overline{X}}_{k,m}=\frac{x_{k,m}}{y_1}\left(m=1,\cdots ,M\right)$$

$${\overline{Y}}_{k,m}=\frac{y_{k,m+1}}{y_1}\left(m=1,\cdots ,M-1\right)$$

$${\overline{F}}_1=\frac{1}{y_1}$$

$${\overline{F}}_m=\frac{f_{m-1}}{y_1}\left(m=2,\cdots ,k+1\right)$$

This system for backing out an internal temperature through the use of a second, rearranged CTF equation is valid regardless of whether the Laplace transform or state space method is utilized to calculate the CTFs and QTFs. The state space method, however, offers a more direct method of obtaining an internal temperature through its definition as an additional output variable.

Consider again the state space example shown in Figure 5. Two output variables were defined for this example: $q_i^{\prime \prime }$ and $q_o^{\prime \prime }$ . The temperature of the node where the source is present can also be defined as an output variable through the identity equation:

$$T_1=T_1$$

When this equation for $T_1$ is added to Equation [eq:LaplaceTransformDeriv659], the resulting output matrix equation for the heat flux at both surfaces and the internal temperature is:

$$\left[\begin{array}{c}{q^{\prime \prime }}_i\\ {q^{\prime \prime }}_o\\ T_1\end{array}\right]=\left[\begin{array}{cc}0& -h\\ h& 0\\ 1& 0\end{array}\right]\left[\begin{array}{c}{T}_1\\ T_2\end{array}\right]+\left[\begin{array}{ccc}0& h& 0\\ -h& 0& 0\\ 0& 0& 0\end{array}\right]\left[\begin{array}{c}{T}_i\\ T_o\\ q_{source}\end{array}\right]$$

The only difference between this relation and Equation [eq:LaplaceTransformDeriv659] is the presence of $T_1$ on both the right and left hand side of the equation. The dual role of $T_1$ as a state variable and an output parameter may seem to contradict the goal of the state space method of eliminating the state variables. However, due to the flexibility of the formulation, nodal temperatures can be extracted in the same manner that any other output quantity would be obtained. For an element with n layers, Equation [eq:QTFDeriv670] becomes:

$$\left[\begin{array}{c}{q^{\prime \prime }}_i\\ {q^{\prime \prime }}_o\\ T_s\end{array}\right]=\left[C\right]\left[\begin{array}{c}{T}_1\\ ⋮\\ T_n\end{array}\right]+\left[D\right]\left[\begin{array}{c}{T}_i\\ T_o\\ q_{source}\end{array}\right]$$

where $T_s$ is the temperature of the node where the heat source or sink is present. The transfer function equation for the calculation of $T_s$ that results from Equation [eq:QTFDeriv671] is identical in form to Equation [eq:LaplaceTransformDeriv653]:

$$T_{s,t}=\sum _{m=1}^M{x}_{k,m}{T}_{i,t-m+1}-\sum _{m=1}^M{y}_{k,m}{T}_{o,t-m+1}+\sum _{m=1}^k{f}_m{T}_{s,t-m}+\sum _{m=1}^M{w}_m{q}_{source,t-m+1}$$

Instead of the flux at either side of the element characterized as a function of temperature, flux, and source history terms, the temperature at the source location is related to source and temperature histories including histories of $T_s$. The validity of these internal temperature calculation methods as well as heat source transfer functions in general will be discussed in the next chapter.

Low Temperature Radiant System Controls[LINK]

The use of this equation allows the low temperature radiant system to be handled like any other surface within the heat balance framework. Heat balances at the inside and outside surfaces take on the same form as other surfaces, and the participation of the radiant system in the radiation balance within the space and thermal comfort models is automatically included. Thus, the radiant system model is fully integrated into the heat balance, and any improvements that are made in areas such as convection coefficients, shading models, etc. are immediately available to the radiant system as part of the overall heat balance solution.

Once the transient nature of the system is accounted for, one must then turn to the next difficult issue: controls. Controls are problematic for almost any simulation program. The problem is not whether something can be simulated because typically a simulation program offers the ability to experiment with many different control strategies. Rather, the problem is typically the diversity of controls that are implemented and keeping the controls that can be simulated up to date. EnergyPlus offers two different flow control schemes: variable flow (ZoneHVAC:LowTemperatureRadiant:VariableFlow) and constant flow/variable temperature (ZoneHVAC:LowTemperatureRadiant:ConstantFlow). The control strategies are different enough that they were developed as separate system types. More details of the controls are described below.

The controls for variable flow low temperature radiant systems within EnergyPlus are fairly simple though there is some flexibility through the use of schedules. The program user is allowed to define a setpoint temperature as well as a throttling range through which the system varies the flow rate of water (or current) to the system from zero to the user defined maximum flow rate. The flow rate is varied linearly with the flow reaching 50% of the maximum when the controlling temperature reaches the setpoint temperature. Setpoint temperatures can be varied on an hourly basis throughout the year if desired.

The controlling temperature can be the mean air temperature, the mean radiant temperature, the operative temperature of the zone, the outdoor dry-bulb temperature, the outdoor wet-bulb temperature, the surface inside face temperature, or the surface interior temperature. The choice of controlling temperature is left to the user’s discretion and set by input as described in the Input Output Reference. For radiant system controls, the operative temperature is calculated as the average of MAT and MRT. The surface inside face temperature is the temperature of the surface in which the radiant system is embedded at the inside face (the side facing the zone being conditioned).

When the user opts to control the radiant system on the surface interior temperature, this temperature is inside the slab itself, and its location is defined using input that describes the construction of the slab (ConstructionProperty:InternalHeatSource–see the Input Output Reference for more details). Note that this user defined temperature still must be at the interface between two layers, but this is easy to overcome by splitting any material into two separate layers. When the user elects to perform a two-dimensional solution, an additional input parameter in the ConstructionProperty:InternalHeatSource object allows the user to chose the location of the user requested temperature in the direction perpendicular to the main direction of heat transfer. This location can be in-line with the hydronic tubing, at the mid-point between two adjacent pipes, or at any node/point in between. Due to the State Space method and a fixed number of nodes in the direction (currently seven), the user’s decimal input for location between one side of the solution domain and the other is converted to a specific node. However, this does allow the user quite a bit of flexibility in the solution and also defining a point on the interior of a radiant system that can be used for controlling a radiant system.

Since flow rate is varied in a variable flow radiant system, there is no explicit control on the inlet water temperature or mixing to achieve some inlet water temperature in a hydronic system. However, the user does have the ability to specify on an hourly basis through a schedule the temperature of the water that would be supplied to the radiant system. Graphical descriptions of the controls for the low temperature radiant system model in EnergyPlus are shown in Figure 9 for a hydronic system. In a system that uses electric resistance heating, the power or heat addition to the system varies in a manner similar to mass flow rate variation shown in Figure 9.

In the constant flow/variable temperature systems, the controls are also considered piecewise linear functions, but in this case the user selects both the control temperatures and the water temperatures via schedules. This offers greater flexibility for defining how the radiant system operates though it may not model every situation. Figure 10 shows how the “desired” inlet water temperature is controlled based on user schedules. The user has the ability to specify the high and low water and control temperature schedules for heating and cooling (separately; a total of eight temperature schedules). Note that this inlet temperature is a “desired” inlet temperature in that there is no guarantee that the system will provide water to the system at that temperature. The model includes a local loop that attempts to meet this demand temperature through mixing and recirculation. Finally, the control temperature options are the same for a constant flow/variable temperature system as they are for the variable flow and electric radiant systems: mean air temperautre, mean radiant temperature, operative temperature, outdoor dry-bulb temperature, outdoor wet-bulb temperature, surface inside face temperature, and surface interior temperature.

The constant flow (variable temperature) low temperature radiant system model is actually a combination of mixing valves, a pump (constant speed, but the maximum flow can be modified by a schedule), and the radiant system (surface, panel, or group of surfaces/panels). This is connected to the main loop through the standard inlet connections as shown in Figure 11. The system controls determine the desired inlet temperature and system flow rate while loop controls determine the flow rate and temperature of the loop. Note that pump heat also factors into the model through a simple constant speed pump model and user input.

There are four possible conditions (separate for heating and cooling). First, if the loop has adequate temperature and flow to meet system requests, then the model sets the radiant system inlet temperature and controls to the desired values based on the controls and simulates. This is the best condition and recirculation and bypass amounts are adjusted accordingly based on radiant system outlet temperatures. Second, if the loop temperature is adequate but the loop flow rate is less than the radiant system flow rate, we may or may not be able to meet the desired inlet temperature since recirculation might lower the temperature below the desired temperature. In this second case, the model first simulates the radiant system with the desired conditions and then resimulates it to solve for the actual inlet temperature (see later in this section) if it cannot achieve the desired inlet temperature. Third, if the loop flow is greater than the radiant flow but the temperature of the loop is not adequate, then there is no amount of mixing that will solve this problem. All of the radiant flow comes from the loop and the loop temperature (after pump heat addition) becomes the radiant system inlet regardless of the temperature controls. Finally, if both the temperature and the flow of the loop are inadequate, then the model simply solves for the actual radiant system inlet temperature and does not try to meet the controls (merely tries to get as close as physically possible given the loop conditions).

One remaining challenge is the merging of the low temperature radiant system model with an integrated building simulation program. In the past, most simulation programs have simulated the building envelope, the space conditioning systems, and the central plant equipment in three separate steps. While this had some advantages and was partly due to a lack of computing capacity, the large drawback for this arrangement is that there is no feedback from the space conditioning system or central plant response to the building conditions. Thus, if the system or plant was undersized, it was reported as an “unmet load” and does not affect the temperatures experienced within the building. IBLAST, a predecessor (Taylor 1991) to EnergyPlus, resolved this issue by integrating all three major components of a building simulation and thus allowing feedback between the equipment and the building envelope.

This integration was not a trivial task and required that the systems be simulated at shorter time steps in some cases to maintain solution stability. In essence, the system simulation will shorten its time step whenever it senses that conditions are changing too rapidly. While this is effective in maintaining solution stability, it can present problems for a radiant system. The radiant system has either a direct or an indirect impact on the surfaces within a building. So, it must be simulated with the building envelope. Yet, it is also a space conditioning system that must act on the space like any other system and thus must also be simulated at the system time step, which can be less than the building time step and can also vary within EnergyPlus.

This issue was handled using a multi-step approach. In EnergyPlus, the heat balance is always simulated first. When this happens, the radiant system is temporarily shut-off to find how the building would respond if there was no heat source/sink. Then, as the system and plant are simulated at multiple shorter time steps, the radiant system is allowed to operate per the controls specified by the user. Flow rate is allowed to vary at each system time step, and the radiant system model is simulated at each time step as if the current flow rate was being used throughout the entire zone time step. This means that each time the heat source/sink in the radiant system is varied during the system simulation the zone heat balance must be recomputed to see what the reaction of the rest of the zone is to this change in the conditions of one (or more) of the surfaces.

In reality, this is not physically correct because each change in the flow rate throughout the system simulation will have an impact on the system time steps remaining before the heat balance is simulated during the next zone time step. Yet, other approaches to solving the mismatch between the system and the zone response of radiant systems are not feasible. One could force the system to run at the same time step as the zone, but this could result in instabilities in other types of systems that might be present in the simulation. On the other hand, one could try to force the zone to run at the shorter time steps of the system, but this could lead to instability within the heat balance due to limits on the precision of the conduction transfer function coefficients.

Despite the fact that the simulation algorithm described above may either over- or under-predict system response dependent on how the system has been controlled in previous system time steps, it is reasonable to expect that the effect of these variations will balance out over time even though it might lead to slightly inaccurate results at any particular system time step. The long-term approach is also in view in the final simulation step at each zone time step. After the system has simulated through enough system time steps to equal a zone time step, the radiant system will rerun the heat balance using the average heat source/sink over all of the system time steps during the past zone time step. This maintains the conservation of energy within the heat balance simulation over the zone time steps and defines more appropriate temperature and flux histories at each surface that are critical to the success of a conduction transfer function based solution. A graphical picture of this somewhat complex multiple step simulation is shown in Figure 12.

Heat Exchanger Formulation for Hydronic Systems[LINK]

As has been mentioned previously, the actual heat transferred between the building element and the hydronic loop is related to the temperature of the building element at the source location as well as the water system inlet and outlet temperatures. In EnergyPlus, it is assumed that the inlet temperature to the slab (defined by a user schedule and the plant simulation) and the mass flow rate (determined by the control algorithm) are known and that the remaining parameters must be calculated. However, the heat balance equations require the heat transferred to the building element from the water loop in order to calculate the heat transferred from the element to the building environment.

Even though systems defined by this model can vary somewhat, the same characteristic link between the system variables exist. For modeling purposes, the overall water/slab system can be thought of as a heat exchanger. While in principle there are two alternative heat exchanger methodologies, it is more convenient to use the effectiveness-NTU method in this case.

Several assumptions will be incorporated into the heat exchanger analysis. It is assumed that the building element that contains the hydronic loop is stationary and that its temperature along the length of the tubing is constant. The latter part of this assumption stems from assumptions made in both the one and two dimensional heat source transfer function derivations. In either case, the source was added at a single node that was characterized by a single temperature. For consistency, this assumption must be made again in the heat exchanger analysis. Another assumption for the current EnergyPlus model is that the fluid in the tubing is water. Additionally, it is assumed that the thermal properties of the water do not vary significantly over the length of the tubing and that the water flows at a constant flow rate. Finally, the temperature at the inside surface of the water tubing is assumed to be equal to the temperature at the source location. In other words, it is assumed that the water tubing itself has no appreciable effect on the heat transfer process being modeled.

Using these assumptions and the effectiveness-NTU heat exchanger algorithm, several equations can be defined which establish the relationship between the heat source and the water temperatures. First, a heat balance on the water loop results in:

$$q={\left(\dot{m}{c}_p\right)}_{water}\left(T_{wi}-T_{wo}\right)$$

where q is the energy transferred between the water loop and the building element, $\dot{m}$ is the mass flow rate of the water, $c_p$ is the specific heat of the water, $T_{wi}$ is the inlet water temperature, and $T_{wo}$ is the outlet water temperature.

The maximum amount of heat transfer that is possible according to the Second Law of Thermodynamics is:

$$q_{max}={\left(\dot{m}{c}_p\right)}_{water}\left(T_{wi}-T_s\right)$$

where $q_{max}$ is the maximum amount of energy transfer that is possible and $T_s$ is the temperature at the source location.

The effectiveness of the heat exchanger, $\epsilon$, is defined as the ratio of the actual energy transfer to the maximum possible, or:

$$\epsilon \equiv \frac{q}{q_{max}}$$

For a heat exchanger where one fluid is stationary, the effectiveness can be related to NTU, the number of transfer units, by the following equation (Incropera and DeWitt 1985):

$$\epsilon =1-e^{-NTU}$$

where NTU is defined by:

$$NTU\equiv \frac{UA}{{\left(\dot{m}{c}_p\right)}_{water}}$$

EnergyPlus includes two possible methods for calculating the heat exchange between the circulating fluid and the rest of the system and these are controlled by user input: ConvectionOnly and ISOStandard. Both are described in more detail below.

ConvectionOnly Heat Transfer Model. In this model, the water tubes were assumed to have no effect on the heat transfer process. So, the only term present in the overall heat transfer coefficient, UA, is a convection term. Thus, the equation for UA is:

$$UA=h\left(\pi DL\right)$$

where h is the convection coefficient, D is the interior tube diameter, and L is the total length of the tube.

The convection coefficient can be obtained from internal flow correlations that relate the Nusselt dimensionless number to other flow properties. For laminar flow in a tube of constant surface temperature, the Nusselt number is defined by:

$$N{u}_D=\frac{hD}{k}=3.66$$

where k is the thermal conductivity of the water.

For turbulent internal flow, the Colburn equation can be used to define the Nusselt number:

$$N{u}_D=\frac{hD}{k}=0.023R{e}_D^{4/5}P{r}^{1/3}$$

where Pr is the Prandtl number of water and $R{e}_D$ is the Reynolds number which is defined by:

$$R{e}_D=\frac{4\dot{m}}{\pi \mu D}$$

The parameter m is the absolute viscosity of water. For internal pipe flow, the flow is assumed to be turbulent for ReD $\ge$ 2300.

ISOStandard Heat Transfer Model. In this model, the convection between the fluid and the tubes and the conduction resistance of the tubes themselves are assumed to have an impact on the overall heat transfer. The method used in EnergyPlus for this model is based on ISO Standard 11855-2(2012) which is entitled “Building environment design — Design, dimensioning, installation and control of embedded radiant heating and cooling systems — Part 2: Determination of the design heating and cooling capacity“. The key equations from the ISO Standard are summarized in Appendix B of Standard 11855-2 and are outlined below.

In ISO Standard 11855-2, the heat convection for turbulent flow between the circulating fluid and the inside wall of the tube is characterized by the following equation:

$$R_{conv}=\frac{W^{0.13}}{8\pi }{\left(\frac{D_{out}-Din}{\dot{m}L}\right)}^{0.87}$$

where W is the tube separation in the direction perpendicular to the main heat transfer direction, $D_{out}$ is the outside tube diameter, $D_{in}$ is the inside tube diameter, $\dot{m}$ is the fluid mass flow rate in the tube, and L is the tube length.

The resistance of the pipe wall is characterized by this equation:

$$R_{tube}=\frac{Wln\frac{D_{out}}{D_{in}}}{2\pi k_{tube}}$$

where $k_{tube}$ is the thermal conductivity of the tube material.

Finally, the U-value for this model is calculated using:

$$U=\frac{1}{R_{conv}+R_{tube}}$$

This U-value is then used in Equation [eq:NTURadSysEquationGeneric] to calculate NTU in the heat exchanger formulation shown above.

Knowledge of the flow conditions allows Equations [eq:RadSysHX675] through [eq:RadSysHXReD681] to be calculated. This essentially eliminates $\epsilon$ as an unknown. The controls and the plant define the water mass flow rate and the inlet water temperature, leaving two equations (Equations [eq:RadSysHX673] and [eq:RadSysHX674]) and three unknowns. The third equation that can be used in conjunction with Equations [eq:RadSysHX673] and [eq:RadSysHX674] is Equation [eq:QTFDeriv672], which is the CTF/QTF equation for the temperature at the source location.

Knowing the inlet water temperature and water mass flow rate, the calculation procedure is somewhat involved and requires, in addition to Equations [eq:QTFDeriv672], [eq:RadSysHX673], and [eq:RadSysHX674], the use of a modified form of Equation [eq:LaplaceTransformDeriv653]. Equation [eq:LaplaceTransformDeriv653] is the standard conduction transfer function formula for a building element with an embedded source/sink of heat. In EnergyPlus, the surface flux on the left hand side of the equation is replaced with a surface heat balance:

$$\left[\begin{array}{c}Surface\\ Heat\\ Balance\end{array}\right]=\sum _{m=1}^M{X}_{k,m}{T}_{1,t-m+1}-\sum _{m=1}^M{Y}_{k,m}{T}_{3,t-m+1}+\sum _{m=1}^k{F}_m{q^{\prime \prime }}_{1,t-m}+\sum _{m=1}^M{W}_m{q}_{source,t-m+1}$$

The surface heat balance includes terms for incident solar energy, radiation heat transfer from internal heat sources such as lights and electrical equipment, radiation between surfaces using Hottel’s Gray Interchange concept, and convection to the surrounding air. The presence of the surface temperature in the heat balance does not pose any problems since Equation [eq:RadSysHX682] will be rearranged to solve for this temperature. Since the radiation heat balance is dependent on conditions at the other surfaces, an iteration loop is required to provide a more accurate estimate of the radiative exchange within the building. This is not the case with the mean air temperature. An assumption of the heat balance is that the mean temperature of the surrounding air is equal to the final air temperature of the previous time step. Using this estimate in the heat balance avoids a second iterative loop around the radiative iteration loop.

Thus, the terms of the heat balance on the left hand side of the equation have been set with the only unknown quantity being $T_i$, the inside surface temperature at the current time step. On the right hand side of Equation [eq:RadSysHX682], most of the terms are already defined since they depend on known values from previous time steps (temperature, flux, and source histories). The only terms which are not defined are the inside surface temperature ($T_i$), outside surface temperature ($T_o$), and internal heat source/sink ($q_{source}$) of the current time step.

The outside surface temperature will depend on the type of environment to which it is exposed. For example, if the surface is a slab on grade floor, the outside surface temperature is defined as ground temperature and does not require an outside surface heat balance. If the element is an interior surface which has both surfaces exposed to the same air space, the outside surface temperature is redefined to be equal to the inside surface temperature. In cases where the outside surface temperature is not simply defined such as a surface exposed to the exterior environment, a heat balance similar to Equation [eq:RadSysHX682] is required to define the outside surface temperature. However, to again avoid iteration, the heat balance equation for the outside surface assumes that conditions at the inside surface were the same as the previous time step. In most cases, since the influence of the current inside surface temperature on the outside surface temperature is very small, this is a valid assumption. In cases were the inside surface temperature has a significant effect, an approximate inside surface heat balance which defines the inside surface temperature is used. This approximate inside balance uses mean air and radiant temperatures from the previous time step.

At this point in the simulation algorithm then, all of the terms in Equation [eq:RadSysHX682] have been defined except the value at the current time step of the inside surface temperature and the heat source/sink. Thus, Equation [eq:RadSysHX682] can be rewritten in a simpler form:

$$T_{i,t}=C_1+C_2{q}_{source,t}$$

where the variable $C_1$ includes surface heat balance and past history terms as well as the influence of the current outside temperature. The term $C_2$ will depend on the heat source transfer function term and the coefficients of terms linked directly to $T_{i,t}$.

Equation [eq:QTFDeriv672], which was the CTF/QTF equation for the temperature at the source location, can be simplified in a similar manner. Grouping the temperature and source history terms which are known quantities together with the effect of the outside surface temperature which is defined as described above, the original equation (Equation [eq:QTFDeriv672] repeated here):

$$T_{s,t}=\sum _{m=1}^M{x}_{k,m}{T}_{i,t-m+1}-\sum _{m=1}^M{y}_{k,m}{T}_{o,t-m+1}+\sum _{m=1}^k{f}_m{T}_{s,t-m}+\sum _{m=1}^M{w}_m{q}_{source,t-m+1}$$

can be reduced to:

$$T_s=C_3+C_4{q}_{source,t}+C_5{T}_{i,t}$$

where $C_3$ includes all of the history terms and the effect of the current outside temperature, $C_4$ is the heat source transfer function for the current time step, and $C_5$ is the conduction transfer function for the inside surface temperature at the current time step.

Substituting Equation [eq:RadSysHX683] into Equation [eq:RadSysHX684] and noting that $q_{source,t}$ is the same quantity as q in Equations [eq:RadSysHX673] and [eq:RadSysHX674] results in:

$$T_s=C_3+C_{4}q+C_5\left(C_1+C_{2}q\right)$$

When this equation is combined with Equation [eq:RadSysHX674], the heat source, which results from a known water inlet temperature, can be shown to be:

$$q=\frac{T_{wi}-C_3-C_1{C}_5}{\frac{1}{\epsilon {\left(\dot{m}{c}_p\right)}_{water}}+C_4+C_2{C}_5}$$

With both q and $T_{wi}$ known, it is a trivial matter to calculate $T_{wo}$ and$T_s$ from Equations [eq:RadSysHX683] and [eq:RadSysHX684] , respectively. Even though the coefficients in Equation [eq:RadSysHX686] are fairly complex, the final equation relating the heat source directly to inlet water temperature is compact and does not require any iteration. As with flux control, once the heat source/sink is defined, the inside surface heat balance can be performed to determine the surface temperatures.

It should be noted that Equations [eq:RadSysHX683] through [eq:RadSysHX686] are a slight simplification of the actual implementation in EnergyPlus. The development shown above follows the heat balance conventions that assume previous values of the inside temperature to calculate the outside temperature. This, in reality, is not necessary and since the radiant system can be significantly influenced by the delay that such an assumption might cause, the initial implementation of radiant systems in EnergyPlus used a development (shown below) that does not lag either the inside or the outside surface temperature. In effect, we can establish three basic equations for the temperature at the inside and outside surface as well as at the location of the heat source/sink:

$$T_{inside}=C_a+C_b{T}_{outside}+C_C{q}^{\prime \prime }$$

$$T_{outside}=C_d+C_e{T}_{inside}+C_f{q}^{\prime \prime }$$

$$T_{source}=C_g+C_h{q}^{\prime \prime }+C_i{T}_{inside}+C_j{T}_{outside}$$

where:

T${}_{inside}$ is the temperature at the inside surface

T${}_{outside}$ is the temperature at the outside surface

T${}_{source}$ is the temperature within the radiant system at the location of the source/sink

C${}_a$ is all of the other terms in the inside heat balance (solar, LW exchange, conduction history terms, etc.)

C${}_b$ is the current cross CTF term

C${}_c$ is the QTF inside term for the current heat source/sink

C${}_d$ is all of the other terms in the outside heat balance (solar, LW exchange, conduction history terms, etc.)

C${}_e$ is the current cross CTF term (should be equal to C${}_b$)

C${}_f$ is the QTF outside term for the current heat source/sink

C${}_g$ is the summation of all temperature and source history terms at the source/sink location

C${}_h$ is the QTF term at the source/sink location for the current heat source/sink

C${}_i$ is the CTF inside term for the current inside surface temperature

C${}_j$ is the CTF outside term for the current outside surface temperature.

Equations [eq:RadSysHX687] and [eq:RadSysHX688] above can be solved to remove the other surface temperature. Substituting the new equations for T${}_{inside}$ and T${}_{outside}$ as a function of C and q" into the equation for T${}_{source}$ and simplifying results in the following equation:

$$T_{source}=C_k+C_l{q}^{\prime \prime }$$

where:

$$C_k=C_g+\frac{C_i\left(C_a+C_b{C}_d\right)+C_j\left(C_d+C_e{C}_a\right)}{1-C_e{C}_b}$$

$$C_l=C_h+\frac{C_i\left(C_c+C_b{C}_f\right)+C_j\left(C_f+C_e{C}_c\right)}{1-C_e{C}_b}$$

Combining this with heat exchanger analysis as shown above, we eventually arrive at the following equation to relate the flux to the slab to the water inlet temperature and mass flow rate:

$$q^{\prime \prime }=\frac{T_{water,in}-C_k}{\frac{C_l}{A}+\frac{1}{\epsilon {\left(\dot{m}{c}_p\right)}_{water}}}$$

which includes all of the inside and outside heat balance terms (“hidden” in the C${}_k$ and C${}_l$ coefficients). Once the flux to the slab is known, the remaining terms of interest (outlet water temperature, inside and outside surface temperatures, etc.) can be calculated using the relatively simpler equations shown above.

Note that the above development is valid for both the hydronic (variable flow) low temperature radiant system and for constant flow (variable temperature) low temperature radiant systems where the inlet temperature is known (based on controls). However, when due to loop conditions and the presence of recirculation, it is not possible to know the inlet temperature to the radiant system without simulating it, we must either iterate or perform more mathematics to arrive at the inlet temperature. The implementation in EnergyPlus chose to avoid iteration and solved for the inlet temperature as shown in the next paragraphs.

The previous equation combines with the following equation which is valid for an surface in the current radiant system:

$$q_j={\dot{m}}_j{c}_p(T_{water,in}-T_{water,out,j})$$

where q${}_j$ is the heat transfer to the jth surface in the radiant system, ${\dot{m}}_j$ is the mass flow rate only to this surface, and T${}_{water,out,j}$ is the outlet temperature for the jth surface. Combining the previous two equations results in:

$$(T_{water,in}-T_{water,out,j})=\frac{{\epsilon }_j\left(T_{water,in}-C_{k,j}\right)}{1+\frac{{\epsilon }_j{\dot{m}}_j{c}_p{C}_{l,j}}{A}}$$

If for each surface in the radiant system, we let:

$$C_{m,j}=\frac{{\epsilon }_j}{1+\frac{{\epsilon }_j{\dot{m}}_j{c}_p{C}_{l,j}}{A}}$$

then the previous equations because the slightly less complex:

$$(T_{water,in}-T_{water,out,j})=C_{m,j}\left(T_{water,in}-C_{k,j}\right)$$

Rearranging to obtain the outlet temperature for the jth surface:

$$T_{water,out,j}=\left(1-C_{m,j}\right)T_{water,in}+C_{m,j}{C}_{k,j}$$

The overall outlet temperature from the radiant system is just a simple mixing of all of the surface outlet temperatures based on flow fraction and results in T${}_{water,out}$.

An energy balance on the mixing valve-pump group results in the following equation that relates the radiant system inlet temperature (T${}_{water,in}$) to the loop inlet water temperature (T${}_{loop,in}$), the radiant system outlet temperature (T${}_{water,out}$), and the pump heat addition:

$$T_{water,in}=\frac{{\dot{m}}_{loop}}{{\dot{m}}_{system}}{T}_{loop,in}+\frac{{\dot{m}}_{recirc}}{{\dot{m}}_{system}}{T}_{water,out}+\frac{PumpHeat}{{\dot{m}}_{system}{c}_p}$$

Plugging in the definition of T${}_{water,out}$ based on the summation of T${}_{water,out,j}$ equations results in:

$$T_{water,in}=\frac{LoopTerm+RecircTerm}{TwiCoeff}$$

where:

$$LoopTerm=\frac{{\dot{m}}_{loop}}{{\dot{m}}_{system}}{T}_{loop,in}+\frac{PumpHeat}{{\dot{m}}_{system}{c}_p}$$

$$RecircTerm=\frac{{\dot{m}}_{recirc}}{{\dot{m}}_{system}}\sum _{j}FlowFractio{n}_j{T}_{water,out,j}$$

$$TwiCoeff=\left(1-\frac{{\dot{m}}_{recirc}}{{\dot{m}}_{system}}\right)\sum _j\left[FlowFractio{n}_j\left(1-C_{m,j}\right)\right]$$

Once the actual water inlet temperature is calculated with this equation, it is then possible to calculate individual outlet temperatures for each surface, the overall outlet temperature, and finally all of the necessary flow and loop quantities. This procedure avoids iteration but is somewhat complex to follow. However, this second mathematical process is only needed for select cases of the constant flow radiant system when the inlet temperature is not known explicitly. With the proper establishment of input data, it can be avoided.

NOTE: In some “low-energy” applications, it is possible that during cooling mode the elevated water temperatures may result in a heat source (or net heating) to the radiant system. To avoid heating the slab when cooling is desired, EnergyPlus performs a temperature limit check. If heating would result during cooling mode or cooling during heating mode, the model will cut-off the flow rate until the inlet water temperature will produce the proper effect.

High Temperature Radiant Heater Model[LINK]

The input object ZoneHVAC:HighTemperatureRadiant provides a model for a high temperature radiant heater that is intended to encapsulate an entire class of heating devices that seek to heat the occupants within a zone by direct radiation. This encompasses a wide variety of heaters including both gas-fired and electric. In most cases, the heater appears much like a lamp or a tube that is suspended from the ceiling of a space, and the surface temperatures are high enough that the heaters must be a safe distance away from the occupied portion of the space for safety concerns.

In EnergyPlus, the high temperature radiant heater model allows the user a reasonable amount of flexibility. Rather than specifying an exact location for the radiant heater(s), the user is allowed to specify the percentage of heat leaving the heater as radiation and then on which surfaces this radiation is incident. In addition, the user is also allowed the ability to define what fraction of radiation leaving the heater is incident directly on a person within the zone for thermal comfort purposes. This amount of heat is then used in the thermal comfort models as shown in Equations [eq:FangerQr] for Fanger, [eq:PierceQr] for Pierce Two-Node and [eq:KSU2NodeQdry] for KSU Two-node. These equations are similar in form to the equation promoted by Fanger (1970). The input parameters for the high temperature radiant heater model are shown below.

$$T_{radiant}={\left[\left(T_{MRT}^4\right)+\left(\frac{Q_{heater\to person}}{\sigma A_{person}}\right)\right]}^{0.25}$$

HIGH TEMP RADIANT SYSTEM,! Program keyword for high temp. radiant heaters

Zone 1 Radiant Heater, ! Zone name

Radiant Operation, ! Availability schedule

SHOP ZONE, ! Zone name (name of zone system is serving)

10000, ! maximum power input (in Watts)

GAS, ! type of heater (either gas or electric)

0.85, ! combustion efficiency (ignored for elec. heaters)

0.75, ! fraction radiant

0.05, ! fraction latent

0.05, ! fraction lost

OPERATIVE, ! temperature control type (MAT, MRT also possible)

2.0 , ! heating throttling range (in C)

Heating Setpoints, ! schedule of heating setpoint temperatures

0.05, ! fraction of radiant energy to people

Zn001:Flr001, 0.75, ! surface/fraction of radiant energy incident on it

Zn001:Wall001, 0.05, ! surface/fraction of radiant energy incident on it

Zn001:Wall002, 0.05, ! surface/fraction of radiant energy incident on it

Zn001:Wall003, 0.05, ! surface/fraction of radiant energy incident on it

Zn001:Wall004, 0.05; ! surface/fraction of radiant energy incident on it

The input for the high temperature radiant heater has two additive relationships that are assumed. First, the fractions of radiant, convective, latent, and lost heat must sum to unity. The user is required to enter the fractions radiant, latent, and lost with the remainder assumed to be convective energy. The fraction latent is added to the latent energy balance and will affect moisture levels within the zone. The fraction lost is assumed to have no impact on the energy balance of the zone and is assumed to be lost or vented to the exterior environment.

The second additive relationship is within the distribution of the radiant fraction. This energy is distributed to people and to the surfaces within the zone. The sum of all of these distribution fractions (the last six lines of input shown above) must sum to unity. Note that each high temperature radiant heater is allowed to distribute energy to up to 20 surfaces and that radiant energy placed on a surface using these distribution fractions is assumed to be completely absorbed. Thus, the distribution fractions should also take into account any differences in long wavelength absorptivity among the surfaces.

Several things should be noted about the fraction of heat that is radiated directly to people. This parameter is somewhat sensitive and will have a direct impact on the thermal comfort models. This is exactly the intent of the high temperature radiant heaters; however, one must use caution when determining this fraction since overestimation of this number might lead to predictions of thermal comfort where in fact it does not exist. In addition, this fraction of radiant energy to people does not have a direct impact on any of the surface heat balances. The thermal comfort energy balance is completely separate from and has no bearing on the zone air or the surface heat balances. Thus, in order to not “lose” this amount of energy from the perspective of the zone air or the surface heat balances, the model assumes that any radiation from the high temperature radiant heater that is incident directly on people is accounted for in the thermal comfort model using Equation [eq:HighTempRadSysTradiant] but is also assumed to be added to the zone air heat balance via convection from people to the surrounding air. This guarantees that the people within the space feel the direct radiative effect of the heaters and that this quantity of energy is not “lost” within the heat balance routines.

Many of the control and integration aspects of the high temperature radiant system model in EnergyPlus are very similar to the low temperature radiant system model. The controls are the same as shown in “Figure 9. Variable Flow Low Temperature Radiant System Controls” where the amount of heat generated by the radiant heater varies as a function of the difference between the controlling and the setpoint temperatures. As with the low temperature radiant system, the controlling temperature is allowed to be the mean air, the mean radiant, or the operative temperature, and the setpoint temperature is allowed to vary hourly based on a user defined schedule. (Operative temperature for radiant system controls is the average of MAT and MRT.) Also, since the high temperature radiant heater has a direct impact on the surfaces within a zone, the surface heat balances are recalculated to determine an approximate response to the radiation from the heater. A final “average” heat balance calculation is done after all of the system time steps have been simulated to maintain continuity within the surface heat balances. The algorithm shown in “Figure 12. Resolution of Radiant System Response at Varying Time Steps is also used for high temperature radiant heaters.

Simple Cooling Panel[LINK]

Overview[LINK]

The simple chilled ceiling model (ZoneHVAC:CoolingPanel:RadiantConvective:Water) is similar to the high temperature radiant system and the hot water radiant/convective baseboard models and was developed to make the modeling of certain radiant systems easier. The existing chilled beam model in EnergyPlus, while useful for modeling active chilled beam systems, is a convective heat transfer only model that is not suitable for all radiant panels or passive chilled beam systems.

Modeling Considerations/Assumptions[LINK]

The model is intended to be simple so several assumptions have been made regarding the performance of this system type:

• Water is flowing through a pipe and that the water has a known inlet temperature defined by the plant loop. The outlet water temperature is unknown and must be calculated as part of the simulation process.

• No information regarding the shape or material of the chilled ceiling panel is required as this is intended to be a simple model only.

• The air conditions of the space are known quantities and though the presence of this system will potentially impact this temperature, it will be assumed that the air temperature is “constant” for the time step in question.

• The entity being modeled is essentially a heat exchanger with flowing water on one side and air on the other.

• Since the water is changing temperature but the air is “not”, we effectively have a similar sort of heat exchanger formulation as the low temperature radiant heating system, but here the water temperature varies while the air temperature is “constant”. Thus, an effectiveness-NTU heat exchanger formulation is best and allows the calculation of total heat gain to the panel and outlet water temperature.

• The total net heat loss to the zone is determined using a user defined coefficient to calculate a zone temperature. This “fraction” radiant parameter defines what portion of the zone temperature is based on the zone MRT with the remainder determined by the zone MAT. While this user defined coefficient remains constant throughout the simulation in keeping with the simple nature of this model, it is not the same as defining a fixed radiant/convective split since this will vary with changing MRT and MAT. More information on this is given below.

• The load met by the system is included in the zone air heat balance just as with any other system. The impact of the system includes the convective load, the radiant gain that is distributed to people (assumed that it is converted into a convective loss), and the increased convection loss from the zone surfaces.

• The radiant loss/gain is distributed to the surfaces as per the radiant surface list that lists which surfaces the radiant cooling panel will interact with and what fraction of radiation each surface will “see”. This is similar to the existing radiant-convective baseboard model and these fractions are user inputs.

• The model is an HVAC level model but interacts with the zone heat balances much like the radiant system and radiant/convective baseboard models.

• Controls are based on controlling water flow rate (like the hydronic radiant system) to meet the thermostatic setpoint temperature of the unit with control types options including MAT, MRT, OpTemp, ODB, or OWB much like the low temperature radiant systems. Controls based on standard air system thermostatic controls are also possible.

• In the event that condensation is predicted on the surface of a panel, the model uses the same condensation control options as the low temperature radiant system. These include: ignoring the problem and allowing the simulation to continue, shutting the system down when condensation is predicted, and a variable response that tries to keep the system running at a lower output while avoiding condensation conditions.

• In order to avoid any violations of the Second Law of Thermodynamics or heating by the system, the cooling panel will be shut off if the inlet water temperature is above the zone temperature.

Algorithm Details: Unit Performance and Calculation of the UA Value[LINK]

As stated above in the assumptions, the most prudent simple method for calculating the performance of a radiant cooling panel is to use a heat exchanger style formulation. Since the air within the space (and the other surfaces of the space) will presumably be at a “constant” temperature while the water flowing through the panel will change temperature, we have a situation where one fluid is changing temperature and one remains at the same temperature. This is very similar to the low temperature radiant system models which use the effectiveness-NTU method for calculating heat exchange between the water in the loop and the slab in which the pipes are embedded. A summary of how this is modeled for the radiant cooling panel follows in the next paragraphs.

Several heat balances can be constructed that relate to the situation at hand. First, a heat balance on the water loop results in:

$$q={\left(\dot{m}{c}_p\right)}_{water}(T_{wi}-T_{wo})$$

where q is the energy transferred between the water loop and the surrounding space to which the cooling panel is exposed, $\dot{m}$ is the mass flow rate of the water, $c_p$ is the specific heat of the water, $T_{wi}$ is the inlet water temperature and $T_{wo}$ is the outlet water temperature.

The maximum amount of heat transfer that is possible according to the Second Law of Thermodynamics is:

$$q_{max}={\left(\dot{m}{c}_p\right)}_{water}(T_{wi}-T_z)$$

where $q_{max}$ is the maximum possible energy transfer and $T_z$ is the temperature of the surroundings (assumes that the air and the unconditioned surfaces in the space are at the same temperature for the sake of this development and the calculation of a UA; more explanation of $T_z$ is given in the next section).

The ratio of the actual energy exchange to the maximum possible is known as the effectiveness of the heat exchanger, $\epsilon$:

$$\epsilon =\frac{q}{q_{max}}$$

For the case when one fluid is stationary in a heat exchanger, Incropera and DeWitt (1985) state that the effectiveness of a heat exchanger can be related to the number of transfer units (NTU) through the following relation:

$$\epsilon =1-e^{-NTU}$$

and NTU is defined by:

$$NTU=\frac{UA}{{\left(m{c}_p\right)}_{water}}$$

The above equations will be used to model the performance of the panel within EnergyPlus and also to obtain a “design” UA value for the panel.

Calculating the UA value of a cooling panel or coming up with a relationship for such a term is not a simple matter. However, most devices are either rated or tested by the manufacturer under certain conditions. From this rating information, a UA value can be calculated. The EnergyPlus model requests rating information (inlet water temperature, mass flow rate, “air” temperature, and panel output under rating conditions) and then uses it to calculate a UA value that is used to model the performance of the panel.

The development of the method for estimating the panel UA value uses the same equations shown above but simply with the purpose of calculating UA rather than the performance of the device, i.e. determining the output of the panel. The method for estimating the UA value of the panel begins with the panel rating data that includes the following values:

• Rating Inlet Water Temperature $\left(T_{rwi}\right)$

• Rating Air/Surroundings Temperature $\left(T_{rz}\right)$

• Rating Water Mass Flow Rate $\left({\dot{m}}_r\right)$

• Rated Output Capacity $\left(Q_r\right)$

To find UA, first rearrange the equations above to make NTU a function of effectiveness as follows:

$$\epsilon =1-e^{-NTU}$$

$$NTU=-\ln (1-\epsilon )$$

where $\epsilon =Q_r/q_{m}ax$ and $q_{max}$ is calculated using design conditions. This leads to:

$$NTU=-\ln (1-\frac{Q_r}{{\left({\dot{m}}_r{c}_p\right)}_{water}(T_{rwi}-T_{rz})})$$

Finally, substituting the definition of NTU as shown above, one can obtain an equation for UA that is calculated directly from the rated values that the user supplies:

$$UA=-{\left({\dot{m}}_r{c}_p\right)}_{water}\ln (1-\frac{Q_r}{{\left({\dot{m}}_r{c}_p\right)}_{water}(T_{rwi}-T_{rz})})$$

This UA value is used by the simple model to calculate the output of the radiant cooling panel based on the current water inlet temperature and mass flow rate as well as the space temperature (combination of air and radiant temperature) as described below.

Algorithm Details: Zone Temperature for Radiant Panel Output and “Fraction Radiant” Input[LINK]

The output from the panel will be based on the mean air temperature (MAT) and the mean radiant temperature (MRT) from the space. These two temperatures are combined to create a “zone” temperature that is used in the model to calculate the radiant cooling panel output. As shown in the previous section, a UA value is calculated using rated conditions. This UA value is assumed to be constant and used in conjunction with the effectiveness-NTU method described in the previous section to determine the output of the panel based on water inlet temperature, water mass flow rate, and space conditions. The space conditions take into account both the mean air temperature and the mean radiant temperature via the following equation:

$$T_z=x_{r}MRT+(1-x_r)MAT$$

where $x_r$ is the “Fraction Radiant” for the radiant cooling panel that the user provides.

This is then used in the equation for qmax above to calculate the maximum possible panel output. Knowing UA and knowing what the water mass flow rate for the system is based on controls, the NTU parameter can be calculated. Once NTU is calculated, the effectiveness can be computed and used with the maximum possible panel output to calculated the actual panel output (q) from:

$$q=\epsilon q_{max}$$

Thus, the modeling equation to find the panel output is:

$$q=(1-e^{\frac{-UA}{{\left(\dot{m}{c}_p\right)}_{water}}}){\left(\dot{m}{c}_p\right)}_{water}(T_{wi}-T_z)$$

Using the results of the above equation, the outlet water temperature can also be calculated using a rearranged form of the water loop heat balance equation given at the beginning of this section.

There has been some concern expressed that the “Fraction Radiant” value ($x_r$) used to calculate the zone temperature ($T_z$) is difficult to estimate and that the correct modeling of the performance of certain models, particularly the high temperature radiant heater, are highly dependent on having a good estimate for this value. The high temperature radiant heater is a high output device that is highly radiant and the control of these high output devices can be difficult to simulate. The radiant cooling panels discussed here are generally not high output devices. A reasonable estimate for $x_r$ generally falls within a range of 0.6-0.7 when calculating $T_z$ above. While this value could vary for different scenarios, choosing a value within this range should be a reasonable approximation of the performance of a radiant cooling panel and should not result in significant inaccuracies in the model.

References[LINK]

Ceylan, H.T. and G.E. Myers. 1980. Long-time solutions to heat conduction transients with time-dependent inputs. ASME Journal of Heat Transfer, Volume 102, Number 1, pp. 115-120.

Degiovanni, A. 1988. Conduction dans un “mur” multicouche avec sources: extension de la notion de quadripole. International Journal of Heat and Mass Transfer, Volume 31, Number 3, pp. 553-557.

Fanger, P.O. 1970. Thermal Comfort-Analysis and Applications in Environmental Engineering, Danish Technical Press, Copenhagen.

Hittle, D.C. 1981. Calculating building heating and cooling loads using the frequency response of multilayered slabs. Ph.D. Thesis, Department of Mechanical and Industrial Engineering, University of Illinois at Urbana-Champaign and Technical Manuscript E-169, United States Army Construction Engineering Research Laboratory, Champaign, IL.

Hittle, D.C. and R. Bishop. 1983. An improved root-finding procedure for use in calculating transient heat flow through multilayered slabs. International Journal of Heat and Mass Transfer, Volume 26, Number 11, pp. 1685-1693.

Hottel, H.C. and A.F. Sarofim. 1967. Radiative Transfer, McGraw-Hill, New York.

Incropera, F.P. and D.P. DeWitt. 1985. Introduction to Heat Transfer. New York: John Wiley & Sons.

Lee, J. and R.K. Strand. 2001. “An Analysis of the Effect of the Building Envelope on Thermal Comfort using the EnergyPlus Program”, submitted for publication in the proceedings of the 2001 ACSA (Association of Collegiate Schools of Architecture) Technology Conference, Austin, TX.

Liesen, R.J. and C.O. Pedersen. 1997. “An Evaluation of Inside Surface Heat Balance Models for Cooling Load Calculations”, ASHRAE Transactions, Volume 103, Part 2.

Maloney, D. 1987. “Development of a radiant heater model and the incorporation of thermal comfort considerations into the BLAST energy analysis program”, M.S. thesis, University of Illinois at Urbana-Champaign, Department of Mechanical and Industrial Engineering.

McClellan, T.M. and C.O. Pedersen. 1997. “Investigation of Outside Heat Balance Models for Use in a Heat Balance Cooling Load Calculation Procedure”, ASHRAE Transactions, Volume 103, Part 2.

Pedersen, C.O., D.E. Fisher, and R.J. Liesen. 1997. “Development of a Heat Balance Procedure for Cooling Loads”, ASHRAE Transactions, Volume 103, Part 2.

Pedersen, C.O., D.E. Fisher, J.D. Spitler, and R.J. Liesen. 1998. Cooling and Heating Load Calculation Principles, ASHRAE.

Seem, J.E. 1987. Modeling of heat transfer in buildings. Ph.D. Thesis, University of Wisconsin-Madison.

Strand, R.K., and C.O. Pedersen. 1994. “Analytical verification of heat source transfer functions”, First Joint Conference of International Simulation Societies, Zürich, Switzerland.

Strand, R.K. 1995. “Heat source transfer functions and their application to low temperature radiant heating systems”, Ph.D. dissertation, University of Illinois at Urbana-Champaign, Department of Mechanical and Industrial Engineering.

Strand, R.K. and C.O. Pedersen. 1997. “Implementation of a Radiant Heating and Cooling Model into an Integrated Building Energy Analysis Program”, ASHRAE Transactions, Volume 103, Part 1.

Strand, R.K. and C.O. Pedersen. 2001. “Modularization and Simulation Techniques for Heat Balance Based Energy and Load Calculation Programs: the Experience of the ASHRAE LOADS Toolkit and EnergyPlus”, International Building Performance Simulation Association, Conference Proceedings of Building Simulation 2001, Rio de Janeiro, Brazil.

Taylor, R.D., C.O. Pedersen, D. Fisher, R. Liesen, and L. Lawrie. 1991. “Impact of simultaneous simulation of building and mechanical systems in heat balance based energy analysis programs on system response and control”, International Building Performance Simulation Association, Conference Proceedings of Building Simulation 1991, Nice, France.