Electric thermal storage (ETS) devices can be used for grid demand load-leveling and off-peak domestic space heating (DSH). A high-resolution three-dimensional finite element model of a forced air ETS heater core is developed and employed to create a general charge/discharge model. The effects of thermal gradients, air flow characteristics, material properties, and core geometry are simulated. A simplified general stove discharge model with a single time constant is presented based on the results of the numerical simulations. This simplified model may be used to stimulate economic/performance case studies for cold climate communities interested in distributed thermal energy storage.


Electric heat from resistive “baseboard” heaters is a common form of domestic space heat (DSH) in cold climates. These systems are relatively inexpensive to purchase, easy to install, simple to control, and clean (from the end-user's standpoint). However, they require an instantaneous supply of electricity from the grid. When a large number of units are installed, the aggregate has the effect of potentially creating large peaks in grid demand during cold spells and the utility may struggle to avoid a capacity deficit.

A previous modeling effort [1] investigated controlling the temperature set points of electric heaters without storage to reduce the peak behavior of thermostatic loads using programmable communicating thermostats. This method is less effective for smaller thermal masses and requires additional installed communication infrastructure.

Perhaps, the most widely accepted solution to electric heat demand load leveling is the use of electric thermal storage (ETS) heaters (see Fig. 1). ETS heaters allow the consumer to decouple their residence's heat demand from the availability of electrical capacity on the grid by storing heat in a dense core of ceramic storage media. Heat can be released thermostatically by forcing air (or a liquid heat transfer fluid) through the core and into the residence as needed. Charging may occur during off-peak hours, or in accordance with excess power availability. On small, wind-powered microgrids, charging will primarily occur during periods of high wind availability.

Typical ETS units will maintain a target state of charge (SOC), which can be manually or automatically set. In general, the SOC should be higher when more heat demand is anticipated (e.g., periods of colder weather). As the core's SOC drops with the average stove temperature, the ability for the ETS unit to provide heat at full output is diminished. Hence, a key design consideration is the discharge behavior of the stove as it cycles from a full SOC to a dead, empty state.

An important design consideration for islanded wind-powered microgrids is total ETS fleet size. Case studies can be found in the literature [35], which are conducted in order to adequately size a fleet of ETS units for a particular grid. A surplus of ETS units generally means a smaller fraction of off-peak power, or wind energy, per consumer. A deficit of ETS units reduces the aggregate storage capacity of the fleet. Though there are various methods for optimization, a good understanding of the transient behavior of the particular ETS device is an advantage in these analyses.

One particular study [6] integrates a one-dimensional finite difference core model into the performance analysis, avoiding additional details in exchange for computational efficiency. For an even simpler analysis, Sateriale [7] makes the assumption that the stove outlet power is constant and independent of the SOC. Some effort has been made in the literature outlining approaches to developing ETS discharge models [8], though they leave it to the designer to follow such methods. Clearly, a simple yet accurate charge/discharge model of a typical room ETS unit is needed to further research and implementation in this area.

This study develops a generalized lumped thermal model for a typical forced air room ETS core. Factors such as core geometry, material, and air flow characteristics are studied using a numerical finite-element model. The results are summarized and outline the transient behavior of a general ETS core, given any combination of these factors. Potential additional variables such as stove thermal insulation and air chamber mixing are omitted for reasons of generalization.


The primary objective of the present work is to develop a generalized performance model of a typical forced air room ETS heater core. With decoupled charge and discharge commands, there are an infinite number of potential input patterns. These patterns could result in numerous time–temperature histories within the storage media. The first and primary step of the analysis is to determine the degree to which certain factors affect the charge and discharge behavior of the core.

The charging behavior alone is not given extensive attention in this study, and the heat addition is modeled as a constant heat flux into an insulated simulation region. This is an accurate description as long as the voltage is supplied consistently and the impedance of the heating elements does not change with temperature. (Some slight increase in impedance will occur as the core temperature rises.) A linear accumulation of sensible thermal energy is the result. However, the discharge characteristics of the core are apt to be a function of the core material properties, the air flow characteristics, and the distribution of thermal energy within it.


Two important unitless quantities can be defined for use in the analysis: SOC and storage disparity (SD). In the case of sensible thermal energy, the SOC is defined as 

where T is the temperature at any point in the storage volume, V; T0 and Th represent the low and high limits, respectively, on the possible range of temperatures within the core. The density and specific heat of the solid are assumed not to be functions of temperature and, therefore, cancel in the ratio. The state of charge is, thus, a nondimensional quantity.

The storage disparity represents the level to which thermal gradients in the core exist. It is a function of the maximum and minimum temperatures within the solid as well as the lowest possible core temperature, T 

When the temperature at every point in the core is the same (perfect thermal equilibration), SD = 0. If the minimum temperature approaches T, SD approaches unity. This metric can be used as a means of quantifying the degree of thermal equilibrium as SD varies between 0 and 1.


A time-variant (transient) COMSOL© Multiphysics® model was developed based on the general layout and brick geometry of the Steffes 2102 room unit. Figure 2 shows a cross section of the 2102 with a box encircling the simulation region. The core consists of eight layers of ceramic storage media with paired air channels between each layer. Four sets of resistive electric heating elements reside in alternate air channels, providing heat flux to the core. A baffle divides the inlet side from the outlet, forcing inlet air through the depth of the core, around a 180 deg turn, and out through the opposite side.

A representative simulation region was chosen in order to reduce the size of the numerical model (Fig. 3). The region encompasses two insulated brick layers and two air channels (one with heating elements, one without). A periodic boundary condition was applied to the upper and lower surfaces of the region to represent a repeating stack of regions. The front, back, and side surfaces are insulated (Q˙  = 0) and heat is generated on the charging surfaces at a constant heat flux. Conduction in the ceramic media and convection in the air flow channels are the operative modes of heat transfer throughout the simulation region. Fluid velocity profiles inside the channel are specified.

Governing Equations.

Transient heat conduction in the ceramic core (subscript “c”) is governed by Fourier's law, expressed here in terms of the ceramic material properties (density (ρ), specific heat (Cp), thermal conductivity (k)), heat flux (Q˙), temperature (T), and time (t) 
Similarly, heat transfer in the air flow channels (subscript “a”) is described with an additional convective term 

where u is the fluid velocity vector. The fluid velocity field is specified, and various profiles are considered.

Boundary Conditions.

Apart from the insulated (Q˙ = 0) and periodic boundaries, the inlet, outlet, and charging surface boundary conditions are defined as follows: 

where T is the room air temperature and h is the convective heat transfer coefficient. If the thermostatically controlled blower is off, h = 0 at the inlet and the boundary is insulated. When the blower activates, drawing air into the channels, h is set to infinity in order to represent an isothermal boundary at Ta=T. Both channel outlet boundary conditions are advective outflow, which can also be modeled using h = 0 combined with a finite value of u.

At the charging surfaces, the heat flux is given as 

In this study, the specific heat flux is 2660 W/m2, which is the rated power of the Steffes 2102.

Material Properties.

Typical parameters for the materials in an actual Steffes 2102 were used for the initial simulation and are specified in Table 1. These values are varied between realistic low and high limits in the parameter sweeps of Sec. 3 and offer a representative starting point for analysis of the transient core SOC. Air properties including thermal conductivity, density, and specific heat are analytic functions of temperature provided within the software.


The initial part of the analysis used the baseline numerical model with the typical material properties from Table 1 in order to gain insight into the stove's general discharge behavior as a function of time.

Solution Linearization and Air Velocity Profile.

Given the nominal material parameters (Table 1), the output of the stove core can be modeled for varying degrees of bulk air flow, which will directly affect the discharge rate of the core. A range of 2–8 m/s was investigated based on laboratory measurements of the Steffes 2102. The initial temperatures are specified at every point in the storage material to the maximum allowable temperature in the core (900 K) for an initial SOC of 100%. Initially, a flat velocity profile is specified in the air channel.

The resulting surface temperature of the simulation region is shown in Fig. 4 (for u = 8 m/s), several minutes after the initiation of a discharge. Cold room air enters the left side and leaves on the right after doubling back in the rear of the stove. The coldest region is, expectedly, the area around the inlet channels. A localized hot spot can be seen in the corner opposite to the inlet side. A strong gradient is visible between the inlet and outlet channels of the brick. The effects of charging-induced gradients are not visible in the bottom channels versus the top, as they quickly disappear after the initiation of a discharge.

Given the numerical temperature solution at each tested air velocity, Eq. (1) can then be used to calculate SOC as a function of time. The resulting stove discharge curves follow an approximate exponential decay pattern (Fig. 5) with a single time constant. The general state of charge can then be modeled as a function of time as


where SOC0 is the initial state of charge (when the discharge command is given) and τ is the time constant. The time constant is determined using a linear regression of ln (SOC) versus t. The time constant (τ) is, therefore, the negative inverse of the slope of the linearized discharge curve.

The results of three simulations in the typical velocity range are shown along with their respective least-squares fit in Fig. 5. Two things are readily apparent: (1) higher velocities increase the discharge rate of the core, and (2) the one-parameter curve fit is less accurate at lower velocities.

To quantify the error between two discharge curves, the mean square error (mse) is defined as 

where (t0, tf) is the interval of time under consideration and the subscripts a and b denote the two curves being compared. The typical interval used in this study will be from t0 = 0 to the time at which the stove has reached SOC = 0.01.

The mse of the three velocities modeled in Fig. 5 are shown in Table 2, in addition to an identical simulation except using a parabolic velocity profile in the air channel. The time constants are equal to within 0.1%, indicating the shape of the air channel velocity profile has negligible effect on the discharge behavior of the stove. The mse data also show that the flat and parabolic velocity curves at each flow rate are nearly identical to each other. Therefore, a flat velocity profile is used for the remainder of the study.

Thermal Gradients.

As the ETS core charges when excess electric power is available (e.g., high wind event) and discharges due to thermostatic demand, thermal gradients can build up within the storage media over time. The extent to which these gradients affect the subsequent discharge behavior of the stove can be quantified. This relationship will depend on the core material, geometry, air flow rate, and charge-discharge history.

Here, the typical parameter values from Table 1 are used for comparison between two cases at various states of charge: immediate and equilibrated discharge. In the case of immediate discharge, a strong thermal gradient is induced by allowing the core to charge from 0 to SOC0 at the full electric charge rate, followed immediately by discharging until SOC = 0 (Fig. 6). The equilibrated case allows for a short period of time after charging to remove thermal gradients before subsequent discharge. The resulting time constants and mse values for these two different simulations are shown in Table 3. The air velocity was 4 m/s in both simulations.

The mse between immediate and equilibrated discharges for each SOC0 tested are very small, on the order of 10−6, indicating minimal correlation between the presence of a thermal gradient and the ETS discharge behavior. However, there is some small variability (<2%) in the value of the time constant for discharge. These results are valid for charge rates at or below 2660 W/m2.

Parameter Sweep.

Thus far, it is has been shown that the air velocity profile in the core channels and the thermal gradients within the media do not significantly affect the discharge behavior of the core. The next analysis will show that the core material properties and core geometry have a much more substantial effect on SOC(t). Furthermore, the bulk air flow velocity (u) will also be shown to have a significant effect as well. A range of the three fundamental parameter values is tested in a sweep study.

Material Properties.

A close look at the governing equation (3) confirms the three core material properties are density (ρ), specific heat (Cp), and thermal conductivity (k). However, the ability for heat to diffuse through the core over time is best described by the thermal diffusivity (α), which is a function of these three properties 

When these individual properties are varied over the range in Table 4, the thermal diffusivity will span a range between 1/8 and 6 times its nominal value.

Core Geometry.

The nominal 4-in (10.2 cm) height of the simulation region H (Fig. 2) is varied within a range of 3–5 in (7.6–12.7 cm) to simulate the effect of varying core geometry (Table 4).

Bulk Air Velocity.

The same range of velocities shown to affect the stove discharge curve (Sec. 2.2) will be analyzed here, in combination with the other fundamental variables (Table 4).

A complete permutation of these three fundamental parameters results in a 33 study (27 simulation points). The results are presented and discussed in Sec. 4.

Results and Discussion

As shown in Sec. 3, three fundamental parameters drive the behavior of the stove core's time-dependent output, which can be modeled as an exponential discharge (7). The material properties, core geometry, and air flow rate can be combined into a single independent parameter and generally used to describe the output of a range of core designs and operating scenarios. This section explores the shape and accuracy of a generalized one-parameter time constant model and an alternative two-parameter (quadratic) exponential model.

One-Parameter Model.

The three fundamental parameters driving the discharge behavior of the core can be represented in a single independent dimensionless quantity of a modified Graetz number similar to that introduced in Ref. [9]. This modified Graetz number (Gz*) is the product of the Reynolds number (Re), the Prandtl number (Pr), and the dimensionless quantity L/DH 

In Eq. (10), u is the bulk fluid velocity. The fluid dynamic viscosity, ν, cancels and does not appear in Gz*. In the case of channel flow, the hydraulic diameter, DH, would be four times the ratio of channel cross-sectional area over perimeter (4A/P). However here, DH cancels and leaves L in Gz*, where L is the length parallel to the flow direction. In order to retain the effect of the finite thickness of the storage medium, a characteristic length L is defined as the ratio of storage medium volume over the heat transfer area inside the channel (effectively the channel length, and reducing to that exactly for small storage medium volumes, analogous to a thin-walled pipe). The thermal diffusivity, α Eq. (9), contains the material properties relevant to heat transport.

This modified Graetz number is the independent variable of the parameter sweep. The study dependent variable describes the output of the stove core. In Eq. (7), the single parameter that describes this output is the time constant, τ.

The time constants derived from the least squares fit (Sec. 3.1) can be calculated for each simulation in the parameter sweep (Sec. 3.3). τ is a function of the convective heat transfer coefficient, h, and the mass (m), specific heat (cp), and heat transfer surface area (A) of the simulation region (11) 
The heat transfer coefficient, h, is present (along with other known simulation parameters) in the Biot number, Bi (12) 

In addition to h, Bi is also a function of the volume of the ceramic within the simulation region (V), the thermal conductivity of the heat transfer media (k), and the total heat transfer surface area (A). The Biot number represents the ratio of heat transfer at the outside of a body to the heat transfer within it. In this study, the Biot number is the dependent variable plotted against the independent variable, Gz*. The results of the parameter sweep are shown in Fig. 7 spanning several powers of ten.

The direct and nearly linear relationship between log(Bi) and log(Gz*) could be expected as the mode of heat transfer at the surface of the channel is entirely advective. Similarly, for low values of Gz* (convective driven heat transport), heat transfer within the body (low values of Bi) overtakes heat transfer on the outside. A linear fit can be applied to Fig. 7, relating Bi to Gz* as 

Note that it is common for a lumped parameter model to be restricted to Biot numbers no greater than 10−1. This accounts for only 10% of the data in Fig. 7. To resolve the error associated with heat gradients within the core, a second parameter can be used.

Two-Parameter Model.

The single parameter model presented in Sec. 4.1 covers a wide range of foreseeable material properties, core geometries, and air flow rates. However, error is introduced at high Biot numbers. This error can be removed with an added quadratic term in the exponential model. The discharge curve is then 

where a and b are the two parameters that describe the fit.

One such improvement is shown in Fig. 8. The one-parameter general fit is a plot of Eq. (7) using a time constant calculated from Eqs. (10)(13). The two-parameter regressive fit uses a least squares method to find the coefficients of Eq. (14). An 84% improvement in the mse is achieved with the additional parameter.

Analysis of the two-parameter model for all simulations indicates that the b parameter is a strong function of u/H. However, a may range from −10−8 to +5 × 10−8 with little correlation to u, H, or α. It is, therefore, recommended that the single time constant model be used for economic feasibility studies and general system performance analyses when high precision accuracy is not necessary to capture the desired behavior characteristics.

Core Energy Balance.

A comprehensive core energy balance (in rate form) can now be stated. The sensible energy (E) within the core is 

where E has unit of energy (Joules) and depends on the density, specific heat, and temperature within the core volume.

Consider such an ETS core subject to two digital signals: the charge command (cc) and discharge command (dc). When one of these values is zero, no energy transfer occurs in its respective mode. For example, if dc = 0, the core is not discharging because the thermostat does not demand heat. A generalized energy balance of the ETS core presented herein can then be expressed as 

In rate form, the storage term (E˙) is the sum of the charge (positive), discharge (negative), and loss terms. The charge rate (first term) is a function of only the charge command (cc) and the charge power (Q˙chrg). The charging heat flux should be within or around the 2660 W/m2 used in this study in order to observe the lumped capacitance assumption.

The controlled discharge rate (middle term) is the product of the discharge command (dc) and the forced discharge rate (time derivative of Eq. (7)), expressed in W 

The time constant (τ) can be obtained from Eqs. (10)(13). The value of E0 becomes the current sensible energy of the core whenever a discharge event initiates (at t0).

The uncontrolled discharge (loss) term in Eq. (16) is a function of the insulation and stove body around the core 

where A is the heat transfer surface area, Tw is the wall temperature of the core, and hstove is a heat transfer coefficient for the entire stove, determined empirically.

Stove Modeling.

Equation (16) can be used to model the charge and discharge cycles of an ETS core given the digital input signals and time constant. Mass and volume calculations can be made using either a representative stove simulation region (as in this paper) or using the entire core, as long as consistency is maintained. The model is compatible with changing room air temperatures, within reasonable proximity to 20 °C, as would be seen with set-back thermostats. Colder room temperatures would tend to cause higher air mass flow rates in the channels and consequently increase the discharge rate.

Some ETS devices may use a variable speed forced air blower for air circulation throughout the core. This changes the bulk air velocity in the channels and alters the time constant. Such an effect would not change the values of dc and cc, which remain binary.

A final consideration in modeling an entire electric thermal stove is air mixing in the outlet manifold. With some designs, such as the Steffes 2102, a small fraction of the forced air is directed through the core, with the majority being used to moderate the outlet temperature. (It is unsafe to allow 900 K air to be released directly into a residence.) As such, the actual air flow rate through the core must be used in the calculation of the core's time constant.


A numerical model of a typical room size ETS heater core was developed and used to study the effects of material, core geometry, air flow rate (and velocity profile), and thermal gradients on discharge rate. Bulk air velocity, core dimensions, and material properties were shown to have the strongest effect on the discharge curve, which can be approximated by a simple exponential time constant and improved with a second parameter fit. Time constants were calculated for a wide range of parameters and summarized against a single nondimensional modified Graetz number. It was shown that the SOC of a typical ETS heater core is only a function of its initial SOC (SOC0) and the core time constant (τ). Heat gradients and air velocity profiles were shown to have a negligible impact. An energy balance equation was given in terms of the findings within the present study.

The information within this study can be put to practical use in the development of economic case studies for ETS applications in cold climates.


The authors would like to thank Dylan Cutler at the National Renewable Energy Laboratory for his thorough review and helpful feedback on this manuscript. An additional thanks to Dennis Meiners at Intelligent Energy Systems for his support of this ETS research.

This work was supported by the U.S. Department of Energy (DOE), Office of Science, Basic Energy Sciences, under Award No. DE-SC0004903.


  • a =

    first coefficient for two-parameter fit

  • A =

    heat transfer surface area

  • b =

    second coefficient for two-parameter fit

  • Bi =

    Biot number

  • cp =

    specific heat

  • cc =

    charge command (digital)

  • DH =

    hydraulic diameter

  • dc =

    discharge command (digital)

  • E =

    energy content of core

  • Gz* =

    modified Graetz number

  • h =

    convective heat transfer coefficient

  • H =

    height of simulation region

  • k =

    thermal conductivity

  • L =

    characteristic length

  • m =

    mass of simulation region

  • Pr =

    Prandtl number

  • Q˙ =

    heat flux

  • Re =

    Reynolds number

  • SD =

    storage disparity

  • SOC =

    state of charge

  • t =


  • T =


  • Th =

    high temperature limit of core

  • Tmax =

    highest temperature in the core

  • Tmin =

    lowest temperature in the core

  • Tw =

    wall temperature of core

  • T0 =

    lowest potential core temperature

  • T =

    room air temperature

  • u =

    bulk air velocity in channel

  • V =


  • γ =

    ratio of specific heats

  • ν =

    fluid dynamic viscosity

  • ρ =


  • τ =

    time constant


D. S.
, “
Tapping the Energy Storage Potential in Electric Loads to Deliver Load Following and Regulation, With Application to Wind Energy
Energy Convers. Manage.
), pp.
Steffes Corporation
, “
Owner's and Installer's Manual for Room Heating Units
,” Steffes Corporation, Dickinson, ND, accessed Apr. 11, 2017, http://www.steffes.com/wp-content/uploads/2016/09/2100SeriesOwnerInstallersManual.pdf
, “
Meeting Residential Space Heating Demand With Wind-Generated Electricity
Renewable Energy
), pp.
, and
, “
Design and Modelling of Dispatchable Heat Storage in Remote Wind-Diesel Systems
,” Renewable Energy Research Laboratory, Department of Mechanical and Industrial Engineering, University of Massachusetts, Amherst, MA, Report No.
, “
The Technical Potential for Off-Peak Electricity to Serve as Backup in Wind-Electric Thermal Storage Systems
Int. J. Green Energy
), pp.
B. Y.
, Moreau, A., Parent, M., and Labrecque, B.,
, “
Central Electric Thermal Storage (ETS) Feasibility for Residential Applications—Part 1. Numerical and Experimental Study
Int. J. Energy Res.
), pp.
M. E.
, “
Modeling and Analysis of Masonry Electro-Thermal Heating and Storage for Optimal Integration With Remote Stand-Alone Wind-Diesel Systems
M.S. thesis
, University of Alaska Fairbanks, Fairbanks, AK.
M. A.
, Lamarche, L., Karjl, S., and Moreau, A.,
, “
Model of Room Storage Heater and System Identification Using Neural Networks
CLIMA 2000 Conference
, Brussels, Belgium, Aug. 30–Sept. 2, Vol.
, pp. 265–277.
, and
Heat Transfer
University Press, Cambridge, UK, p.