From GridLAB-D Wiki
Jump to: navigation, search

 *** WORKING DRAFT *** 
Please review, edit and comment as necessary.


Tech:Multizone_ETP_Linearization - Linearized solution of the Equivalent Thermal Parameters method for modeling multizone commercial buildings

Introduction

 TODO:  Describe purpose of multizone ETP solver and reason for linearization.

Nomenclature

[math]C_n[/math] 
The heat capacity of [math]n[/math]th node (Btu/degF).
[math]\Delta t[/math] 
The time-step taken for the next iteration (h).
[math]\Delta T_n[/math] 
The difference between the observed temperature and the setpoint temperature of the [math]n[/math]th node, if controlled.
[math]k[/math] 
The proportional control gain of the HVAC control system (unitless).
[math]N[/math] 
The number of nodes in a thermal network (integer).
[math]n[/math] 
A reference to the [math]n[/math]th of [math]N[/math] nodes (integer).
[math]m[/math] 
A reference to the [math]m[/math]th of [math]N[/math] nodes (integer).
[math]Mode_n[/math] 
The HVAC mode of the system interacting with the [math]n[/math]th zone (heat, cool, etc.).
[math]Q_n[/math] 
The HVAC heat flow into (positive for heating) or out of (negative for cooling) the [math]n[/math]th node, if controlled (Btu/h).
[math]Q_{nm}[/math] 
The heat flow from the [math]n[/math]th to [math]m[/math]th nodes (Btu/h).
[math]Q_{mode-cap,n}[/math] 
The heating (positive) or cooling (negative) capacity of the HVAC system interacting with the [math]n[/math]th node, if controlled (Btu/h).
[math]T_n[/math] 
The temperature of the [math]n[/math]th node (degF).
[math]T_{cool,n}[/math] 
The cooling setpoint temperature of the [math]n[/math]th node, if controlled (degF).
[math]T_{heat,n}[/math] 
The heating setpoint temperature of the [math]n[/math]th node, if controlled (degF).
[math]U_{nm}[/math] 
The thermal conductance of the path from the [math]n[/math]th and [math]m[/math]th nodes (Btu/degF/h).

Methodology

Figure 1. 6-Node Thermal Network

Known parameter values and boundary conditions are shown in red. Illustrated in Figure 1 is an arbitrary six-node thermal network. Note that the top horizontal branch is the original ETP circuit model, where [math]T_0[/math] is the outdoor air temperature (known boundary condition), [math]T_1[/math] is the indoor air temperature, and [math]T_2[/math] is the mass temperature. Node 3 might represent an unconditioned atrium space. Nodes 4 and 5 are present to illustrate series solutions of massless nodes, and might also represent an atrium with an additional buffer space to the outdoors.

Two additional massless paths are added to illustrate additional circumstances. A node n is characterized by its temperature Tn, its mass [math]C_n[/math] (which can be zero), and an exogenous heat gain [math]Q_n[/math] (from solar, internal, and/or HVAC, which is a boundary condition that can also be zero).

Nodes [math]m[/math] and [math]n[/math] are connected by conductance’s [math]U_{mn}[/math]. A conductance might represent the UA of a building, an air flow rate between two nodes, or the product of a convective surface heat transfer coefficient and the surface area, for example.

To solve this circuit explicitly, the heat flow from any node [math]m[/math] through the conductance to any other node [math]n[/math] at the time [math]t[/math]is denoted as [math]Q_{mn}(t)[/math]. Thus

[math]Q_{mn}(t)=U_{mn}\left[T_m(t)-T_n(t)\right][/math]
(1)

Note that this implies a sign convention for heat flow into a node as positive. Thus, by definition, [math]Q_{mn}=-Q_{nm}[/math].

Assume a "small" time step [math]\Delta t[/math] during which all temperatures change by a "small" amount relative to the temperature differences between the nodes, i.e. [math]Q_{mn}[/math] can be considered a constant. This is the linearizing assumption. The maximum time-step [math]\Delta t[/math] permissible that satisfies the linearizing assumption is limited by two considerations:

(1) How long until an HVAC [math]Mode_n[/math] changes or internal gains [math]Q_n[/math] change?
This is determined by the time to the next expected change in mode, which can be estimated for each node based on the type of control (bang-bang or proportional).
(2) How long until the change in the temperature difference between a two nodes [math]n[/math] and [math]m[/math] exceeds a preset limit?
This is determined by computing the rate of change of the temperature difference and estimating the time until that rate of change exceeds a preset limit.

Step 1: Compute temperatures of all massive nodes

Then the heat balance on a massive node [math]n[/math], i.e., [math]C_n \gt 0[/math], from a set of connected nodes [math]m=1[/math] to [math]M[/math], from time [math]t[/math] to time [math]t+\Delta t[/math] is:

[math]Q_n(t)+\sum_{m=1}^M Q_{mn}(t)=\frac{C_n\left[T_n(t+\Delta t)-T_n(t) \right]}{\Delta t}[/math]
(2)

where [math]Q_n(t)[/math] is assumed to be constant boundary condition over the interval.

Solving this for the temperature of the node at the next time step:

[math]\frac{C_n T_n(t+\Delta t)}{\Delta t} = Q_n(t)+\sum_{m=1}^M Q_{nm}(t)+\frac{C_n T_n(t)}{\Delta t}[/math]
(3)
[math]T_n(t+\Delta t) = \frac{\Delta t}{C_n} \left[Q_n(t)+\sum_{m=1}^M Q_{nm}(t)\right] + T_n(t)[/math]
(4)

Substituting the definition for [math]Q_{nm}[/math] from Equation (1):

[math]T_n(t+\Delta t) = \frac{\Delta t}{C_n} \left[Q_n(t)+\sum_{m=1}^M U_{mn}\left[T_m(t)-T_n(t)\right]\right] + T_n (t)[/math]
(5)

or

[math]T_n(t+\Delta t) = \frac{\Delta t}{C_n} Q_n(t) + \frac{\Delta t}{C_n} \sum_{m=1}^M U_{mn} T_m(t) + T_n(t) \left[1-\frac{\Delta t}{C_n}\sum_{m=1}^M U_{mn}\right][/math]
(6)

Step 2: Compute temperatures of all massless nodes

If node [math]n[/math] is massless, i.e., [math]C_n = 0[/math], then it must be in thermal equilibrium with all adjacent nodes at any time. For massless nodes, Equation (2) reduces to

[math]Q_n(t)+\sum_{m=1}^M Q_{nm}(t+\Delta t)=0[/math]
(7)

where over the time interval from [math]t[/math] to [math]t+\Delta t[/math] [math]Qn[/math] is a constant boundary condition.

Substituting the definition for [math]Q_{mn}[/math] from Equation (1):

[math]Q_n(t)+\sum_{m=1}^M U_{mn} \left[ T_m(t+\Delta t)-T_n(t+\Delta t) \right]=0[/math]
(8)

Solving for the temperature of the node at time t+Δt:

[math]T_n(t+\Delta t)=\frac{Q_n(t)+\sum_{m=1}^M U_{mn}T_m(t+\Delta t)}{\sum_{m=1}^M U_{mn}}[/math]
(9)
Figure 2. Reduced equivalent of 6-node thermal network

Computing the temperature of massless nodes in the thermal network assumes all temperatures of adjacent nodes are known at the end of a time-step, either as boundary conditions, or because they are massive and have had their temperatures computed in Step 1.

The simplest way to resolve this is to reduce the network to an equivalent network when massless nodes are in series, parallel, or "wye" configuration.

Series massless nodes configurations

Figure 3a - Series massless node reduction

A series configuration of two nodes, such as the thermal network in Figure 1, can be reduced as shown in Figure 2, where Node 5 is eliminated and the equivalent series conductance from Node 0 to Node 4 is

[math]U_{xy} = \frac{U_{xw} + U_{wy}}{U_{xw} U_{wy}}[/math]
(10a)

Parallel massless nodes configurations

Figure 3b - Parallel thermal path reduction

A parallel configuration of two nodes can be reduced to a single node thus:

[math]U_{xy} = U_{A} + U_{B} \![/math]
(10b)

Delta massless nodes configurations

A "delta" configuration of 3 nodes [math](x, y, z)[/math] can be transformed to a "wye" configuration of four nodes [math](x,y,z,w)[/math] as follows:

[math]\begin{align} U_{xw} &= U_{xy}U_{xz}\left[ \frac{1}{U_{xy}} + \frac{1}{U_{yz}} + \frac{1}{U_{xz}} \right] \\ U_{yw} &= U_{xy}U_{yz}\left[ \frac{1}{U_{xy}} + \frac{1}{U_{yz}} + \frac{1}{U_{xz}} \right] \\ U_{zw} &= U_{xz}U_{yz}\left[ \frac{1}{U_{xy}} + \frac{1}{U_{yz}} + \frac{1}{U_{xz}} \right] \end{align}[/math]
(10c)

Note that for the general "star-mesh" case, the transformation is:

[math] U_{ij} = U_iw U_jw \sum_{n=1}^N \frac{1}{U_n} \![/math]
(10d)

where [math]N=1[/math] is the dangling node case (which eliminates the node), [math]N=2[/math] is the series node case (see above) and [math]N=3[/math] is the "delta-wye" transformation case. The "star-mesh" transformation has no general inverse without additional constraints, so all mesh configurations must be simplified using a series of appropriate wye-delta transformations.

In such all cases, the temperature of the new node, [math]w[/math] is calculated as follows:

[math] T_w = \frac{\sum_{n=1}^N U_{nw} T_n}{\sum_{n=1}^N U_{nw}} \! [/math]
(10e)

Special Case for "Bang-Bang" HVAC Control

Nodes whose temperature is used to control [math]Q_n[/math] (typically heating or cooling from the HVAC system) must have a modified calculation procedure, because [math]Q_n[/math] is no longer a boundary condition but rather a function of the node temperature, [math]T_n[/math], and the thermostat setpoints [math]T_{cool,n}[/math] and [math]T_{heat,n}[/math].

Define [math]Q_n[/math] as the sum of the HVAC energy input (heating positive, cooling negative), and the sum of internal heat gains and solar heat gains, to Node [math]n[/math]:

[math]Q_n = Q_{hvac,n}+Q_{gains,n} \quad[/math]
(11)

Let [math]Q_{cool-cap,n}[/math] and [math]Q_{heat-cap,n}[/math] be the net cooling and heating capacity available to supply the zone, respectively (net of the fan power). Let [math]Q_{fan,n}[/math] be the fan power when the HVAC is off (i.e. if the fan runs continually then [math]Q_{fan,n} \gt 0[/math]).

Assume a heating and cooling thermostat with setpoints centered in a deadband with range [math]+\Delta T_n[/math] on either side of the setpoint. Note that setpoints must not overlap:

[math]T_{cool,n} - \Delta T_n \gt T_{heat,n} + \Delta T_n \quad[/math]
(12)

The thermostat sets [math] Q_{hvac,n} [/math], which persists in subsequent time steps until changed by the thermostat:

When [math] T_n \lt T_{heat,n} - \Delta T_n [/math]

[math] Mode_n = On \quad ; \quad Q_{hvac,n} = Q_{heat-cap,n} [/math].
(13)

When [math] T_n \gt T_{cool,n} + \Delta T_n [/math]

[math] Mode_n = On \quad ; \quad Q_{hvac,n} = Q_{cool-cap,n} [/math].
(14)

When [math] Mode_n=On [/math] and [math] T_{heat,n}+ \Delta T_n [/math][math] T_n [/math][math] T_{cool,n} - \Delta T_n [/math]

[math] Mode_n = Off \quad ; \quad Q_{hvac,n} = Q_{fan,n} [/math].
(15)

Then the node temperature at time [math]t+\Delta t[/math] can be calculated from Equation (6). Note that all HVAC control nodes with On/Off control must be massive. If one were massless, Equation (9) indicates there is no way to maintain a setpoint without proportional control of [math]Q_{hvac,n}[/math].

Special case for Proportional-Differential (PD) HVAC Control

In the case of proportional control for [math]Q_{hvac,n}[/math], then a proportional-differential control scheme:

When [math] T_{heat,n} - \Delta T_n \lt T_n(t) \lt T_{heat,n} + \Delta T_n \quad[/math]

[math] Q_{hvac,n}(t) = Q_{heat-cap,n}(t) \left[ \frac{T_{heat,n}+\Delta T_n-T_n(t)}{2\Delta T_n} - \frac{k}{2\Delta T_n} \left[ T_n(t)-T_n(t-\Delta t)\right]\right][/math].
(16)

When [math] T_{cool,n} - \Delta T_n \lt T_n(t) \lt T_{cool,n} + \Delta T_n \quad[/math]

[math] Q_{hvac,n}(t) = Q_{cool-cap,n}(t) \left[ \frac{T_n(t) - T_{cool,n} + \Delta T_n}{2 ∆T_n} - \frac{k}{2\Delta T_n} \left[ T_n(t-\Delta t)-T_n(t) \right] \right][/math].
(17)

When [math] T_{heat,n} + \Delta T_n \le T_n(t) \le T_{cool,n} - \Delta T_n[/math]

[math]Q_{hvac,n}(t) = Q_{fan,n} \quad[/math].
(18)

where [math]k[/math] is the proportional gain for the controller.

Step 3: Compute temperatures of reduced nodes

Massless nodes that have been reduced from the network are unnecessary to the simulation model of the network if

  1. the temperature of the node does not affect some non-linear aspect such as a thermostatic control, and
  2. the temperature is not needed as an output variable for some purpose.

In general, it may be best to assume that if the user specified “unneeded” nodes, that there was some purpose in mind, and their temperatures should be computed from Equation (9) as a final step.

Validation

 TODO:  Describe how to validate a numerical implementation of this method.

Authors

This method was developed by Robert G. Pratt and Lucy Huang at Pacific Northwest National Laboratory

See also