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_{}$$