*5.1. Modelling a Heat Accumulator Using an Explicit Finite-Difference Method*

The system of Equations (5) and (8) with boundary condition (10) and initial conditions (11) and (12) was solved using the explicit finite-difference method. An accumulator of length *Lr* was divided into *N* finite volumes of length Δ*x* (Figure 8). By introducing a dimensionless coordinate, the dimensionless length of the finite volume is given by the formula:

$$
\Delta \mathbf{x}^+ = \frac{1}{\mathbf{N}'} \tag{53}
$$

**Figure 8.** Finite difference grid used in air (*Tg*) and packing (*Tw*) temperature calculations.

The derivative of the air temperature *Tg* after the *x*-coordinate in Equation (5) was approximated by the forward difference quotient at the old time step *n*, and the derivative after time by the forward difference quotient for node *i* + 1, which was located at the outlet of the control region contained between nodes *i* and *i* + 1. The remaining terms in Equation (5) were calculated at the old time step *tn* = *n*Δ*t*, *n* = 0, 1, ... . The approximate form of Equation (5); i.e., the difference equation, is as follows:

$$T\_{\mathbf{g},i}^n \frac{T\_{\mathbf{g},i+1}^{n+1} - T\_{\mathbf{g},i+1}^n}{\Delta t} + \frac{1}{N\_{\mathbf{g},i}^n} \frac{T\_{\mathbf{g},i+1}^n - T\_{\mathbf{g},i}^n}{\Delta x^+} = T\_{w,i}^n - \frac{T\_{\mathbf{g},i}^n + T\_{\mathbf{g},i+1}^n}{2} \tag{54}$$

The finite-difference grid is defined as follows:

$$\mathbf{x}\_{i} = (i - 1)\Delta \mathbf{x}\_{\prime} \ \Delta \mathbf{x} = \mathbf{L}\_{\tau} / N\_{\prime} \ \mathbf{x}\_{i}^{+} = \mathbf{x}\_{i} / L\_{\tau} \ \mathbf{i} = 1, \ldots, (N + 1) \ t\_{n} = n\Delta t \ \mathbf{i} \ n = 0, 1, \ldots \tag{55}$$

The differential Equation (8) was approximated in a similar way:

$$\tau\_{w,i}^{n}\frac{T\_{w,i}^{n+1} - T\_{w,i}^{n}}{\Delta t} = -\left(T\_{w,i}^{n} - \frac{T\_{\mathbf{g},i}^{n} + T\_{\mathbf{g},i+1}^{n}}{2}\right)i = 1, \ldots, (N+1)\ n = 0, 1, \ldots \tag{56}$$

Solving Equation (54) for *Tn*+<sup>1</sup> *<sup>g</sup>*,*i*+<sup>1</sup> gives:

$$T\_{\mathbf{g},i+1}^{n+1} = T\_{\mathbf{g},i+1}^{n} + \frac{\Delta t}{\mathbf{r}\_{\mathbf{g},i}^{n}} \left[ \left( T\_{\mathbf{w},i}^{n} - \frac{T\_{\mathbf{g},i}^{n} + T\_{\mathbf{g},i+1}^{n}}{2} \right) - \frac{1}{N\_{\mathbf{g},i}^{n}} \frac{T\_{\mathbf{g},i+1}^{n} - T\_{\mathbf{g},i}^{n}}{\Delta \mathbf{x}^{+}} \right], i = 1, \ldots, N \text{ } n = 0, \ldots \tag{57}$$

The wall temperature *Tn*+<sup>1</sup> *<sup>w</sup>*,*<sup>i</sup>* was determined using Equation (56):

$$T\_{w,i}^{n+1} = T\_{w,i}^n - \frac{\Delta t}{t\_{w,i}^n} \left( T\_{w,i}^n - \frac{T\_{\emptyset,i}^n + T\_{\emptyset,i+1}^n}{2} \right), i = 1, \dots, N \; n = 0, 1, \dots \tag{58}$$

The system of difference Equations (57) and (58) was solved with the boundary condition (10), which takes the form:

$$T\_{\mathbb{S}^1} = \left(T\_{\mathbb{S}}'\right)^n \tag{59}$$

The initial conditions for gas and packing are as follows:

$$T^0\_{\mathcal{G},i} = T\_0 \; \; i = 1, \dots, N+1 \tag{60}$$

$$T\_{w,i}^0 = T\_0 \; i = 1, \dots, N \tag{61}$$

From Equations (57) and (58), the air and packing temperatures were determined as a function of location and time while taking into account the boundary condition (59) and the initial conditions (60) and (61). The time step Δ*t* should satisfy the Courant–Friedrichs– Lewy condition [36]:

$$\frac{\Delta t}{N\_{\ $}r\_{\$ }\Delta \alpha^{+}} \le 1\tag{62}$$

Due to the small value of the product *ρgcpg*, the heat accumulated in the air is minimal. Therefore, air temperature changes over time occur very quickly. Disregarding the heat accumulation in the air; i.e., assuming *τ<sup>g</sup>* = 0, a simplified form of the differential Equation (5) is as follows:

$$\frac{1}{N\_{\mathcal{S}}} \frac{\partial T\_{\mathcal{S}}}{\partial \mathbf{x}^{+}} = T\_{\mathcal{W}} - T\_{\mathcal{S}} \tag{63}$$

Replacing the derivatives in Equation (63) with the difference quotient yields:

$$\frac{1}{N\_{\mathbb{S}^{\bar{i}}}^{n}} \frac{T\_{\mathbb{g},i+1}^{n+1} - T\_{\mathbb{g},i}^{n+1}}{\Delta x^{+}} = T\_{w,i}^{n+1} - \frac{T\_{\mathbb{g},i}^{n+1} + T\_{\mathbb{g},i+1}^{n+1}}{2} \tag{64}$$

The solution of Equation (64) for *Tn*+<sup>1</sup> *<sup>g</sup>*,*i*+<sup>1</sup> has the following form:

$$T\_{\mathbf{g},i+1}^{n+1} = \frac{1}{1 + \frac{\Delta \mathbf{x}^{+}}{2N\_{\mathbf{g},i}^{n}}} \left[ \frac{\Delta \mathbf{x}^{+}}{N\_{\mathbf{g},i}^{n}} T\_{\mathbf{w},i}^{n+1} + T\_{\mathbf{g},i}^{n+1} \left( 1 - \frac{\Delta \mathbf{x}^{+}}{2N\_{\mathbf{g},i}^{n}} \right) \right], \ i = 1, \ldots, N \; n = 0, 1, \ldots \tag{65}$$

The packing temperature is calculated using the relationship (58) with boundary condition (59) and initial conditions (60) and (61).

Assuming zero air heat capacity, the computer simulation time is greatly reduced, as there is no Courant–Friedrichs–Lewy stability condition on time-step length Δ*t* (62). At low air velocities, the air temperature rise per control volume is large, and the approximation of the average temperature by the arithmetic mean of the inlet and outlet temperatures, as in Equation (54) or Equation (64), is insufficiently accurate [14]. The accumulator length must be divided into a large number of control volumes; e.g., assume *N* = 50 or more to ensure high calculation accuracy at low air velocities [14].

#### *5.2. Modelling a Heat Accumulator Using the Crank–Nicolson Method*

In the Crank–Nicolson method [36], the derivative after time in the differential equations for air and packing was approximated by the forward difference quotient. The arithmetic mean of the values approximated the remainder of the equation at the beginning and end of a given time step. Equation (5) can be written in the following form:

$$\pi\_{\mathcal{S}} \frac{\partial T\_{\mathcal{S}}}{\partial t} = -\frac{1}{N\_{\mathcal{S}}} \frac{\partial T\_{\mathcal{S}}}{\partial x^{+}} + \left(T\_{w} - T\_{\mathcal{S}}\right) \tag{66}$$

Equation (66) was approximated by the following difference scheme using the Crank– Nicolson method:

$$\begin{split} \frac{1}{2} \left( \boldsymbol{\tau}\_{\boldsymbol{g},i}^{n} + \boldsymbol{\tau}\_{\boldsymbol{g},i}^{n+1} \right) \frac{T\_{\boldsymbol{g},i+1}^{n+1} - T\_{\boldsymbol{g},i+1}^{n}}{\Delta t} &= \frac{1}{2} \left\{ \left[ -\frac{1}{N\_{\boldsymbol{g},i}^{n}} \frac{T\_{\boldsymbol{g},i+1}^{n} - T\_{\boldsymbol{g},i}^{n}}{\Delta x^{+}} + \left( T\_{\boldsymbol{w},i}^{n} - \frac{T\_{\boldsymbol{g},i}^{n} + T\_{\boldsymbol{g},i+1}^{n}}{2} \right) \right] + \\ + \left[ -\frac{1}{N\_{\boldsymbol{g},i}^{n+1}} \frac{T\_{\boldsymbol{g},i+1}^{n+1} - T\_{\boldsymbol{g},i}^{n+1}}{\Delta x^{+}} + \left( T\_{\boldsymbol{w},i}^{n+1} - \frac{T\_{\boldsymbol{g},i}^{n+1} + T\_{\boldsymbol{g},i+1}^{n+1}}{2} \right) \right] \right\} i = 1, \ldots, N, \; n = 0, 1, \ldots \end{split} \tag{67}$$

Equation (8) for the accumulator packing was transformed in a similar manner:

$$\frac{1}{2} \left( \mathbf{r}\_{w,i}^{n} + \mathbf{r}\_{w,i}^{n+1} \right) \frac{T\_{w,i}^{n+1} - T\_{w,i}^{n}}{\Delta t} = -\frac{1}{2} \left[ \left( T\_{w,i}^{n} - \frac{T\_{g,i}^{n} + T\_{g,i+1}^{n}}{2} \right) + \left( T\_{w,i}^{n+1} - \frac{T\_{g,i}^{n+1} + T\_{g,i+1}^{n+1}}{2} \right) \right] \tag{68}$$
 
$$i = 1, \dots, N, \ n = 0, 1, \dots$$

The system of Equations (67) and (68) was solved for *Tn*+<sup>1</sup> *<sup>g</sup>*,*i*+<sup>1</sup> and *<sup>T</sup>n*+<sup>1</sup> *<sup>w</sup>*,*<sup>i</sup>* . To solve the system of Equations (67) and (68), the Thomas method [37], also known as the tri-diagonal matrix algorithm, was used. The Thomas algorithm is a particular case of the Gauss elimination method that does not require inversion of the coefficient matrix. The unknowns in the system of equations were determined using simple analytical formulas. Therefore, the time needed for numerical calculation was very short. The set of Equations (67) and (68) can also be solved with the iterative method of Gauss–Seidel. The temperatures *T<sup>n</sup> <sup>g</sup>*,*i*+<sup>1</sup> and *Tn <sup>w</sup>*,*<sup>i</sup>* were taken as initial values in the iterative process. Due to the insignificant differences between *Tn*+<sup>1</sup> *<sup>g</sup>*,*i*+<sup>1</sup> and *<sup>T</sup><sup>n</sup> <sup>g</sup>*,*i*+<sup>1</sup> and between *<sup>T</sup>n*+<sup>1</sup> *<sup>w</sup>*,*<sup>i</sup>* and *<sup>T</sup><sup>n</sup> <sup>w</sup>*,*<sup>i</sup>* the number of iterative steps required to obtain the solution was small.
