Equivalent Thermal Parameter Model
The model approach that is used to estimate thermal loads is called an equivalent thermal parameter (ETP) modeling approach. Because the ETP approach has been proven to provide reasonably model the residential (and small commercial building) loads and energy consumption and also because it is based on first principles, this modeling approach has been chosen for the current work (Sonderegger 1978; Subbarao 1981; Wilson et al. 1985; Pratt and Taylor 1994).
House
Thermal loads are implemented using the ETP model, which determines the heat flow between an inner and outer air mass and one solid mass embedded within the inner air mass.
The basic ETP parameters are
- \(U_a\): the conductance between the inner and outer air mass
- \(C_a\): the heat capacity of the internal air mass
- \(C_m\): the heat capacity of the internal solid mass
- \(U_m\): the conductance between the inner air mass and the inner solid mass
- \(T_o\): the outer air temperature
- \(T_i\): the inner air temperature
- \(T_m\): the inner mass temperature
- \(Q_a\): the heat flux to the interior air mass
- \(Q_m\): the heat flux to the interior solid mass (this is not supported yet)
The solution to the ETP model is based on the two ordinary differential equations (ODEs)
The general first-order ODEs are found by inspection:
where the constants \(c_1\) through \(c_6\) are
- \(c_1 = -\frac{U_a + U_m}{C_a}\)
- \(c_2 = \frac{U_m}{C_a}\)
- \(c_3 = \frac{Q_a + U_a T_0}{C_a}\)
- \(c_4 = \frac{U_m}{C_m}\)
- \(c_5 = -\frac{U_m}{C_m}\)
- \(c_6 = \frac{Q_m}{C_m}\)
The general form of the second-order ODE is given by
where
- \(p_1 = \frac{1}{c_2}\)
- \(p_2 = -\frac{c_1+c_5}{c_2}\)
- \(p_3 = \frac{c_1 c_5}{c_2}-c_4\)
- \(p_4 = \frac{c_3 c_5}{c_2}-c_6\)
where
- \(r_1\) and \(r_2\) are the roots of \(p_1 r^2 + p_2 r + p_3 = 0\),
- \(k_1 = \frac{r_2 T_{i,0} + r_2 p_4 / p_3 - \dot T_{i,0}}{r_2-r_1}\),
- $ k_2 = \frac{T_{i,0} -r_1 k_1}{r_2}$
- \(t\) is the elapsed time,
- \(T_{i,0}\) is the initial value of \(T_i\), and
- \(\dot T_{i,0} = c_1 T_{i,0} + c_2 T_{m,0} + (c_1 + c_2) T_0 + c_7\)
with
- \(c_7 = \frac{Q_a}{C_a}\), and
- \(T_{m,t} = k_1 \frac{r_1-c_1}{c_2} + k_2 \frac{r_2-c_1}{c_2} - \frac{p_4}{p_3} + \frac{c_6}{c_2}\)
Modeling Assumptions
The house model assumes that only envelope characteristics, solar gain through windows and internal gain contribute to the HVAC load. Only air conditioners and heat pumps are modeled in the current implementation.
Modeling Approach
The electric circuit analog of an ETP model used to simulate heating and cooling loads in a typical residence is shown in Figure 3. The heat transfer properties are represented by equivalent electrical components with associated parameters for modeling the thermostatically controlled heating, ventilation, and air-conditioning (HVAC) system.
where,
- \(C_{air}\) – air heat capacity (Btu/°F or J/°C)
- \(C_{mass}\) – mass (of the building and its content) heat capacity (Btu/°F or J/°C)
- \(UA_{wall}\) – the gain/heat loss coefficient (Btu/°F.h or W/°C) to the ambient
- \(UA_{mass}\) – the gain/heat loss coefficient (Btu/°F.h or W/°C) between air and mass
- \(T_o\) – CLTDc + Tair (°F or °C)
- \(T_{ambient}\) – ambient temperature (°F or °C)
- \(T_{air}\) – air temperature inside the house (°F or °C)
- \(T_{mass}\) – mass temperature inside the house (°F or °C)
- \(Q_{HVAC}\) – heat rate for HVAC (Btu/hr or W)
- \(Q_{internal}\) – heat rate from other appliance, plug loads, lights and people in the residence (Btu/h or W)
- \(Q_{solar}\) – heat gain from solar (Btu/h or W)
A state space description of the ETP model is:
where (there's something important missing here: A, B, C, D are not described)
- \(R_1 = 1/UA_{insul}\)
- \(R_2 = 1/UA_{mass}\)
- \(Q = Q_{HVAC} + Q_{solar} + Q_{internal}\)
Because there is wide diversity in the thermal parameters (UA, mass, efficiency of equipment, and over sizing factor) used to compute the load across homes, ranges can be provided for these parameters. Because of regional difference in construction practices and wide variation even within a region, a range of values can be assigned for critical parameters based on metering studies (Pratt et al. 1990). While estimating energy consumption of individual homes, the thermal parameters are randomly (from the established range using known distribution, for example, uniform) assigned to each home.
The internal gains in the home are computed external to the house object and are assumed to be inputs to the ETP model. These gains represent heat from other appliances (such as range, microwave), plug loads, lights, and people.
Solving the ETP model (simultaneously for Tair and Tmass), we can obtain the cooling/heating load of an individual home as a function of time. The energy consumption to meet the comfort needs in the home can then be computed by converting the thermal loads by using typical manufacturer-provided air conditioner part load performance data. Using the same simulation model, the energy consumption of a population of homes can also be computed. While simulating a population of homes, the energy consumption for each home is computed simultaneously at each time step. Therefore, we get an accurate representation of distribution feeder load when electric loads are aggregated.
A single ETP model with different thermal parameters can be used to represent all homes in the population for simplicity, or if there is a need, multiple ETP models with different thermal parameters can be used. Therefore, in addition to changing input parameters while simulating a population of homes, the ETP model can also be changed to accurately represent a given building stock. This ensures that it accurately reproduces the effects of load diversity.
Limitation and Future Improvements
The approach described in the previous section only accounts for sensible loads. This does not lead to a significant error in heating load calculation (because most heating load is sensible); it does lead to some error in the estimation of the cooling load.
ETP closed form solution
The methodology described below is the same as used in the GridLAB-D™ code base (see Caveats), although, over time additions have been made to allow for greater flexibility. Additionally modifications have been made to speed up the solution process as implemented in code. The solution methodology, as described here and implemented in GridLAB-D™ was validated against the referenced spreadsheet.
For the thermal circuit in Figure 1, a heat balance (conservation of energy) can be written for the air temperature node (\(T_A\)) as:
\(Q_A - U_A (T_A-T_O) - H_M (T_A-T_M) - C_A \frac{dT_A}{dt} = 0 \qquad\qquad(1)\)
The heat balance for the mass temperature node (\(T_M\)) can be written as:
\(Q_M - H_M (T_M-T_A) - C_M \frac{dT_M}{dt} = 0\qquad\qquad(2)\)
As shown in the ETP closed form solution, Equation (1) can be solved for \(T_M\), differentiated with respect to time to provide \(dT_M/dt\), and both of these substituted into (2) to form a second order linear differential equation in \(T_A\) of the form
\(a \frac{d^2T_A}{dt^2} + b \frac{dT_A}{dt} + c\ T_A = d\qquad\qquad(3)\)
where:
- \(a = \frac{C_M C_A}{H_M}\)
- \(b = \frac{C_M (U_A + H_M)}{H_M} + C_A\)
- \(c = U_A\)
- \(d = Q_M + Q_A + U_A T_O\)
which has the solution with known, constant boundary conditions \(T_O\), \(Q_A\), and \(Q_M\) and initial conditions at time \(t=0\) of \(T_{A_o}\) and \(dT_{A_o}/dt\)
\(T_A = A_1 e^{r_1 t} + A_2 e^{r_2 t} + \frac{d}{c}\qquad\qquad(4)\)
where:
- \(r_1 = \frac{-b + \sqrt{b^2-4ac}}{2a}\)
- \(r_2 = \frac{-b - \sqrt{b^2-4ac}}{2a}\)
- \(A_1 = \frac{r_2 T_{A_o} - \frac{dT_{A_o}}{dt} - r_2 \frac{d}{c}}{( r_2 - r_1 )}\)
- \(A_2 = T_{A_o} - \frac{d}{c} - \frac{r_2 T_{A_o} - \frac{dT_{A_o}}{dt} - r_2 \frac{d}{c}}{( r_2 - r_1 )}\)
The initial condition \(T_{A_o}\) is known as the final condition of \(T_A\) from previous time step. However, at any time step at which the boundary conditions \(T_O\), \(Q_A\), or \(Q_M\) have changed (i.e. the weather, internal gains, heating/cooling output) from the previous time interval, then the new air temperature trajectory at the beginning of the time step can be derived from Equation (1) as
Then, differentiating Equation (4) and substituting it and Equation (4) into Equation (1) yields a solution for \(T_M\) of the form:
where:
- \(g = \frac{Q_M}{H_M}\)
- \(A_3 = \frac{r_1 C_A}{H_M} + \frac{U_A + H_M}{H_M}\)
- \(A_3 = \frac{r_2 C_A}{H_M} + \frac{U_A + H_M}{H_M}\)
Initial Room and Mass Air Temperature
When initializing a House_E simulation, the temperature and weather history prior to the first time step is unknown. First assume the house is at the equilibrium air temperature (at steady state with the heating and cooling system off, i.e. in balance). Equilibrium is defined by \(dT_A/dt_o = 0\). Then, by differentiating (4), at the beginning of the initial time step, it can be shown that
and the corresponding mass temperature at equilibrium is
If the heating system would be “on” (based on the thermostat heating set point; see the section Heating/Cooling Thermostat Operations) at the condition \(T_A = T_{A_{eq}}\), then the initial conditions are best approximated as
\(T_{A_o} = T_{M_o} = T_{set\ heat}\qquad\qquad(9)\) # Heating system “on” at \(T_A = T_{A_{eq}}\)
If the cooling system would be “on” (based on the thermostat cooling set point) at the condition \(T_A = T_{A_{eq}}\), then the initial conditions are best approximated as
\(T_{A_o} = T_{M_o} = T_{set\ cool}\qquad\qquad(10)\) # Cooling system “on” at \(T_A = T_{A_{eq}}\)
The time of day when either of these approximate initial conditions is most correct is when the conditions in the house have been stable for a long period of time. So, starting the simulation at midnight may be a good choice. An often better choice would be the earlier of sunrise or just prior to a morning thermostat change.
If neither the heating or cooling system would be “on” at the condition \(T_A = T_{A_{eq}}\), then the initial conditions are best approximated as
\(T_{A_o} = T_{A_{eq}}\qquad\qquad(11)\) # Heating/cooling system “off” at \(T_A = T_{A_{eq}}\) \(T_{M_o} = T_{M_{eq}}\qquad\qquad(12)\)
Derivation of ETP second-order ODE
Solving (2) for \(T_A\) gives:
Differentiating (2.2):
Solving (1.2) for \(T_M\) gives:
Differentiating (1.3):
Rearranging (2.1):
Substituting (1.3) and (1.4) into (2.4) gives:
Rearranging into the form of a differential equation gives:
or
or
with
Indoor air temperature
The indoor air temperature is found by solving equation (3.2). A different solution method is used depending on whether \(d\) is constant with time.
Fixed conditions
Let \(T_A = e^{r t}\). Then:
and
The homogeneous (d=0) solution to (3.4) is
so
From (3.3b):
From (3.3a) and (3.3c):
So
which guarantees that \(r_1\) and \(r_2\) are real and distinct numbers.
Then we have
and
Substituting into (3.2)
The general solution is:
With the known boundary conditions at $ t=0$, \(T_A(0)\) and \(\frac{d}{dt}T_A(0)\)
Then \(A_1\) and \(A_2\) can be found by subtracting (5) from (4) to eliminate \(A_2\):
From (6.1) and (4)
Forced conditions
Given at the time \(t\):
\(T_A(t)\) is the indoor air temperature, and \(T_M(t)\) is the building mass temperature \(T_O(t)\) is the outdoor air temperature, \(Q_A(t)\) is the heat added to the indoor air, \(Q_M(t)\) is the heat added to the building mass,
and the constant values
\(U_A\) is the building envelop conductivity to the indoor air, \(C_A\) is the heat capacity of the indoor air, \(H_M\) is the building mass conductivity to the indoor air, and \(C_M\) is the heat capacity of the building mass.
The heat balance on \(T_A\) is
Rearranging (1a) and (2a) in the form for solving a differential equation:
From above the second-order ODE is
where
Separate (3.3e) into a constant and a linear function of time
where
\(g_0 = U_A T_O(0) + \Delta g \frac{C_M}{H_M} \,\) \(\Delta g = U_A \frac{T_O(\Delta t)-T_O(0)}{\Delta t}\)
and \(T_O(0)\) and \(T_O(\Delta t)\) are the outdoor air temperatures at time \(t=0\) and time \(t=\Delta t\), respectively.
Separate (3.3e) into a constant, a step \(u(t)\), and an impulse \(\delta(t)\) at time 0
where
Define the initial constant
We define the initial indoor air temperature conditions as
\(x_0 = T_A(0) - \frac{1}{U_A}f_0\) and \(\dot x_0 = \frac{d}{dt}T_A(0)\)
Use Laplace transforms convert equation (8) to \(s\)-domain
Rearrange to solve for \(X(s)\)
Expand the partial fractions to get
Match the terms in equations (9) and (9.1)
Solve for \(A\), \(B\), \(C\), and \(D\)
So we have
\(A = a x_0 - \frac{a}{c}\Delta h + \frac{ab}{c^2}\Delta g\) \(B = b x_0 + a \dot x_0 + \frac{C_M}{H_M} - \frac{b}{c}\Delta h + \left ( \frac{b^2}{c^2} - \frac{a}{c} \right ) \Delta g\) \(C = \frac{1}{c}\Delta h - \frac{b}{c^2}\Delta g\) \(D = \frac{1}{c}\Delta g\)
From above we know the quantity \((b^2-4ac)>0\) so we factor
with
and use inverse Laplace transforms on (9.1) to get
which simplifies to
where
Mass temperature
The mass temperature is found by solving equation (1.3) after the air temperature is determine. A different solution method is used depending on whether \(d\) is constant with time.
Fixed conditions
Solving (1.2) for \(T_M\) at any time gives
Substituting (3.8) and (3.6.1) for \(T_A\) and \(dT_A/dt\) gives
Rearranging and simplifying
Let
where
(1.5.2) gives \(T_M\) for the next time step, i.e., at \(t = t_0 + \Delta t\). Note the solution is discontinuous at any time \(t_s\) when \(Q_M\), \(Q_A\) or \(T_O\) change. So at the time \(t_s^+\) immediately after a change (1.3) must be solved for
Forced conditions
Substituting (3.3g) and (3.3h) into (1.5) with $ f_0=g_0+h_0$
From (10) we have
Substitute this and (10) into (11)
Rearrange this to get
Implementation Notes
Initial conditions
The initial indoor air temperature condition \(T_A(0)\) is available from the initial values given in the GLM model. However, the initial derivative of the indoor air temperature is not usually given. The value of \(dT_A/dt(0)\) can be computed based on the initial state of the HVAC system by solving the equation (1.3) for \(dT_A/dt\) at time \(t=0\):
where the conditions \(T_M(0)\), \(T_O(0)\), and \(Q_A(0)\) are all known states of the ETP model at the end of the last time increment.
Solving for time
Solving equation (10) for time must be done using a numerical solver that implements Newton's method. GridLAB-D™ has an ETP solver function built-in for this purpose.
Etp
The data structure for the ETP solver is as follows:
struct etp {
double t; // the time solution to the ETP equation
double a; // the a parameter
double n; // the n parameter
double b; // the b parameter
double m; // the m parameter
double c; // the c parameter
double p; // the precision of the solution
double *e; // the pointer to the error estimate
};
To use the ETP solver in a module, you must define the macro USE_GLSOLVERS in each source file that will makes solver call, including the init.cpp file (which loads and initializes the solver):
#define USE_GLSOLVERS
#include "gridlabd.h"
When the module is initialized, you must load the ETP solver (which will automatically initialize it).
glsolver *etp_solver;
EXPORT CLASS *init(CALLBACKS *fntable, MODULE *module, int argc, char *argv[])
{
if ( set_callback(fntable)==NULL ) return NULL;
etp_solver = new glsolver("etp");
// ...
return first_class;
}
To set/get parameters for the ETP solver, use the set/get calls:
int iteration_limit;
if ( etp_solver->get("max_iterations",&iteration_limit,NULL)==0 )
gl_error("unable to get iteration limit");
if ( etp_solver->set("max_iterations",iteration_limit,NULL)==0 )
gl_error("unable to set iteration limit");
To call the ETP solver, you must construct the ETP data structure required by the solver. For example:
struct { double t,a,n,b,m,c,p,*e; } etp = {0,0.1,-0.1,0.2,-0.2,0.3,1e-6,NULL};
if ( etp_solver->solve((void*)&etp)==0 || isnan(t) )
gl_error("no solution found");
Where max_iterations is used to control the maximum number of iterations used by the solver.
Caveats
Only the indoor air and mass temperature solution for fixed outdoor air and heat gain condition is currently implemented in GridLAB-D™. This requires that outdoor condition and heat gain conditions be linearized for each solution and can lead to small errors when integrating over variable timesteps (which GridLAB-D™ does quite often). This can introduce discrepancies when comparing simulations where the timing of solution updates differ for otherwise identical temperature and heat gain conditions. This problem could be largely remedied by using the forced condition solution, but as noted above the derivation of this solution appears to be incorrect and cannot be implemented at this time.
Sources
-
RG Pratt, ETP Solution , PNNL Engineering Worksheet , 5/27/2008
-
RG Pratt, ETP Solution , Excel Workbook , 5/27/2008 (edited later)
Heat Loss Coefficient ( \(U_A\))
Compute exterior surface areas:
Define the following, based on rectangular geometry. Let x = width, y = depth. Then the aspect ratio is \(R = y / x\) The floor area is \(A = x y n\)
and the volume is \(V = A h\)
The gross exterior wall area (\(A_{wt}\)) can be derived by introducing the perimeter (p), as follows:
| Equation | Explanation |
|---|---|
| \(A_{wt} = 2 n h (1 + R) \sqrt{\frac{A}{nR}}\) | The gross exterior wall area (\(A_{wt}\)). |
| \(A_g = WWR\ A_{wt}\ EWR\) | The gross window area (\(A_g\)). The window area is used to calculate net_exterior_wall_area, envelope_UA, solar_heatgain_factor, and house_content_heat_transfer_coeff. This is an automatically calculated value, not a variable that the user may set. The default window_area is 244.95 square feet. It is calculated by window_area = gross_wall_area * window_wall_ratio * exterior_wall_fraction. Using default values: \(1633\text{ft}^2 \cdot 0.15 \cdot 1 = 244.95\text{ft}^2\). |
| \(A_d = n_d\ A_{1d}\) | The total door area (\(A_d\)), used to calculate net_exterior_wall_area, envelope_UA, and house_content_heat_transfer. This is an automatically calculated value, not a variable that the user may set. The default door area is 78 square feet. It is calculated by door_area = number_of_doors * 3.0 * (78.0/12.0). |
| \(A_w = (A_{wt}-(A_g + A_d))\ EWR\) | The net exterior wall area (\(A_w\)). This is used to calculate envelope_UA and house_content_heat_transfer_coeff. The default net exterior wall area is about 1310 square feet. It is calculated by \(1 \cdot 1633\text{ft}^2 - 244.95\text{ft}^2 - 78\text{ft}^2 = 1310.05\text{ft}^2\). |
| \(A_c = \frac{A}{n} ECR\) | The net exterior ceiling area (\(A_c\)). This is an automatically calculated value, not a variable that the user may set. The default exterior ceiling area is 2500 sqft. It is calculated by exterior_ceiling_area = floor_area * exterior_ceiling_fraction / number_of_stories. |
| \(A_f = \frac{A}{n} EFR\) | The net exterior floor area (\(A_f\)). This is an automatically calculated value, not a variable that the user may set. The default exterior floor area is 2500 sqft. It is calculated by exterior_floor_area = floor_area * exterior_floor_fraction / number_of_stories. |
The total heat loss coefficient (conductance), \(U_A\), for the house (the last term is for air infiltration); the defaults produce \(U_A\) = 522.1 Btu/°F.hr
Note
0.018 is the volumetric heat capacity of air at standard conditions (Btu/°F.ft³std-air-pressure)
Interior Mass Surface Conductance ( \(H_m\))
Surface area is estimates as total exterior walls (less doors and windows) + interior walls + ceilings
Total “Air” Mass ( \(C_a\))
Based on tuning to typical home heating system cycling times, the “air mass” seems to be well approximated as 3 times the volumetric capacitance of the interior air volume.
\(C_a = 3\ (0.018 A\ h)\) # Short-cycle thermal mass
Total Thermal Mass ( \(C_m\))
\(C_m = A\ m_f - 2\ (0.018 A\ h)\) # Thermal mass (daily cycle), less that added to the “air” mass