^{1}Finnish Institute of Occupational Health, Developing Indoor Environments, Finland
^{2}Aalto University, School of Engineering, Department of Civil and Structural Engineering, Finland
Received date: May 14, 2015; Accepted date: June 07, 2015; Published date: June 14, 2015
Citation: Holopainen R, Salonen EM (2015) Regenerator Analysis with Finite Elements. J Fundam Renewable Energy Appl 5:173. doi: 10.4173/20904541.1000173
Copyright: © 2015 Holopainen R , et al. This is an openaccess article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original author and source are credited.
Visit for more related articles at Journal of Fundamentals of Renewable Energy and Applications
The classic onedimensional regenerator problem has been solved numerically in a new way. Instead of marching in time, the hot and cold periods are treated as boundary value problems in rectangular domains. To solve the boundary value problems, a spacetime finite element method using the commercial package COMSOL 4.3 with MATLAB with a Galerkin formulation has been applied. This approach minimizes the coding effort. The code solves a problem automatically with three consecutive meshes of increasing densities. We present results of three example cases in some detail.
Regenerator; COMSOL; MATLAB
Latin symbols
A total heat transfer surface area (m^{2})
c specific heat at constant pressure (k_{j}/k_{gk} )
e_{r} relative error estimate ( − )
overall heat transfer coefficient (w/m^{2}k)
L length of storage unit (m )
m mass (kg )
mass rate of fluid flow
P period time length (s )
r effective convergence rate ( − )
dimensionless temperature ( − )
t temperature (°C ,k )
x axial coordinate (m)
Greek symbols
nondimensional time coordinate ( − )
η_{REG }thermal ratio ( − )
reduced length ( − )
λ grid refinement ratio ( − )
nondimensional position coordinate ( − )
reduced period ( − )
time (s)
Subscripts
c coarse
f fluid, fine
i inlet
m medium
s solid (storage) material
o initial
Superscripts
' Hot period
" Cold period
Space average
Heat recovery is a topical issue, as it has a major effect on buildings’ energy efficiency. A typical heat recovery application is a ventilation system, in which intake air is heated by warm outlet air. This article focuses on regenerative heat exchanger computational modelling, in which mass is alternately heated and cooled by changing the direction of the air flow.
Classic onedimensional counter flow regenerator analysis demands normally the use of numerical methods. In this article, the analysis is performed in a new way by applying the finite element method. This approach has certain advantages from the coding effort and accuracy perspectives. We are not aware of the use of the finite element method as a pure boundary value problem in this connection.
We have based our presentation strongly on the theory and arguments given in Schmidt and Willmott [1]. The same theme has been treated, for example, in a rather recent text on regenerative heat transfer in Willmott [2]. However, this latter reference contains no essential changes with respect to the theme of the present article. We have preferred mostly the notations of Schmidt and Willmott [1].
Basic Equation
As mentioned, we will follow the theory given in Schmidt and Willmot [1]. We first consider the basic equations in the single blow case. The regenerator case – to be described in Section 3 – can then be seen as consisting of more or less of consecutive single blow cases.
The setting is briefly illustrated in Figure 1a. The fluid and solid temperatures to be determined t_{f} and t_{s} depend on device axial position x and time τ. As the initial condition, t_{s} is given on boundary τ = 0 and as the boundary condition, t_{f} is given on fluid inflow boundary x = 0.
Schmidt and Willmott [1] describes the governing energy balance equations and the procedures for transforming the equations into dimensionless forms, as indicated in Figure 1b. The dimensionless fluid and solid temperatures are T_{f} and T_{s} and these depend on the dimensionless position coordinate ξ and dimensionless time η. The dimensionless length measures corresponding to L and P in Figure 1a are ٨ and Π. The mathematical treatment of the problem takes place most conveniently in the dimensionless form.
The governing field equations are
(1)
(2)
Notations R refer to field equation residuals. This terminology is common in the finite element method literature – e.g. (Zienkiewicz and Taylor [3] – and here means that certain formulas to follow can be written in a short form.
The initial and boundary conditions become, respectively:
T_{s} (ξ, 0) = 0 (3)
T_{f} (0, η) = 1 (4)
The standard numerical approach used in e.g. (Schmidt and Willmott [1] is based on the finite difference method. Furthermore, the time coordinate is also taken into account in Schmidt and Willmott [1] in the standard way by “marching” in the time direction in a stepwise manner.
The numerical approach employed in this article is based on the finite element method. In addition, the space and time coordinates are treated equally. Considering the very similar roles of ξ and η in equations (1) and (2), this seems appropriate. So no marching in time is performed. The case shown in Figure 1b is considered as just as a boundary value problem, with boundary conditions (3) and (4) on the two boundary parts η = 0 and ξ = 0, respectively. Thus, we now also alternatively call (3) a boundary condition.
As is well known, the finite element method is based on integral type presentations with the method of weighted residuals. The most common versions of the method of weighted residuals are the Galerkin method and the leastsquares method. We experimented with two Galerkin type versions (Formulations I and II) and with a leastsquares
version (Formulation III). The weak forms corresponding to these formulations are listed below
Formulation I the weak form is
(5)
Formulation II The weak form is
(6)
Formulation III The weak form is
(7)
The integrals are over domain Ω= (0,Λ) × (0,Π) in Figure 1b. The Galerkin method is normally defined as one in which the same basis functions used in the approximations of the unknowns are employed as weighting functions. This can also be interpreted so that the weighting functions (before discretization) are variations (δT_{f}, δT_{s} ) of the unknown functions (T_{f},T_{s} ). The two options for choosing the weighting functions are seen in Formulations I and II. The leastsquares version is based on the functional
(8)
Requiring the functional to have a stationary value produces the weak form (7).
The use of the finite element method may be justified by the fact that certain general shelf software – especially COMSOL Multiphysics with the finite element method – allows the applier to feed his/her problem directly as a weak form. We used this option. Obviously, this results in the required coding effort becoming rather small.
Because of the simple rectangular solution domain, uniform meshes consisting of fournode rectangular Lagrangian elements were used for both unknowns T_{f} and T_{s}.
Figure 2 shows crude meshes for two geometries considered below. It should be noted, however, that the elements are not quite standard, as they are actually spacetime elements. Boundary conditions (3) and (4) are satisfied by giving the corresponding nodal values.
In the generated MATLABcode used in connection with COMSOL Multiphysics, care was taken to gain confidence in the accuracy of the numerical results. We more or less followed the recommendations of references Roache [4]. Sinclair et al. [5] in every calculation case, we used three meshes of increasing densities: coarse, medium and fine, and we use the corresponding subscripts c, m and f respectively when referring to the appropriate values. If h_{c} is the typical mesh size of a coarse mesh, the corresponding mesh sizes for the medium and the fine meshes are h_{m} = h_{c}⁄λ and h_{f} = h_{c}⁄λ^{2} , respectively, where λ is the refinement ratio or the scale factor [5]. We employ the value λ=2.
The effective convergence rate r for quantity f is obtained from
(9)
Following (Roache, 2009), the relative error estimate expression e_{r} for result f with the finest mesh is
(10)
where the multiplier 1.25 is a “factor of safety”. Use of the Richardson extrapolation was not considered necessary.
The case sketched in Figure 2a is taken from Example 2.1 in Schmidt and Willmott [1] and the case in Figure 2b is the same as the first case but with the domain measures switched. The crude meshes shown in the figures are selected so that the elements are roughly squareshaped. The meshes used in the first case are thus 4×8, 8×16 and 16×32, and in the second case 8×4, 16×8 and 32×16.
As mentioned above, the single blow case was considered first, as the situation closely resemble a typical regenerator case period, and even more so because an analytical solution is available that making use of the zero order modified Bessel function of the first kind [6], so that here, the actual errors of the numerical solutions could be determined.
We recorded solid temperature results at point A (Figure 2), the fluid and solid temperatures at point B and the average solid temperature on line AB for the three formulations described above. All three formulations worked very well. However, Formulation I generally seemed to provide the most consistent results. The magnitude of its actual maximum relative (percentage) error with the fine meshes, calculated using the analytical solution, was 0.07%. In addition, error estimate e_{r} bounded from above the actual error. Based on the results obtained, only Formulation I was selected for use in the regenerator calculations below.
In a counter flow regenerator, hot and cold periods follow cyclically. As in [1], we used a single prime to signify the hot period, and double primes for the cold period. As the names indicate, inflow temperature t_{fi} ' of the fluid during the hot period is higher than inflow temperature t_{fi} '' of the fluid during the cold period. The periods here had fixed time lengths P' and P'', so that the length of one total cycle was P'+P''. We considered only the simplest case, in which mass flow rates and were constants in time, similarly to t_{fi} ' and t_{fi} ''. We approached the periodic solution by going through enough cycles beginning with an arbitrary given distribution t_{s} (x,0).
In the terminology of Schmidt and Willmott [1], this approach is an open method, as opposed to the closed method. The process began with the hot period taking the flow direction in the positive xaxis direction, so that the flow inflow boundary condition was given at x=0. The temperature distribution in solid material t_{s} ' (x, P’), obtained at the end of the hot period, was used as the initial (or boundary) condition for the cold period. In reality, there is a short filling period between the hot and cold period (as well as between the cold and hot period) until the governing equations are valid once more. We ignored the filling period and started the cold period analysis immediately after τ=P'. Using dimensionless formulations, the computational situation is thus described schematically in Figure 3.
For a hot period, the governing dimensionless field equations are (see (1) and (2))
(11)
(12)
with the boundary condition
(13)
Similarly, for a cold period the field equations are
(14)
(15)
with the boundary condition
(16)
Equation (14) shows a change of signs. Here ξ' and ξ'' are taken to grow in the same direction (Figure 3). However, mass flow rate m_{f} '' is considered positive even when the flow direction is from right to left in the cold period. This is the reason for the change of signs in (14).
Using the dimensionless formulation, the initial (or boundary) conditions at the start of the hot and cold periods took the following forms, respectively,
(17)
(18)
The solution for each period was obtained by Formulation I, naturally taking into account the changes in signs in (14) and the righthand side boundary condition (16).
The code proceeded by solving the equations consecutively in hot and cold periods until error limit δ' Λδ'' <0.001 was reached. The error was checked at each period with previous and and present value using and values using expression
(19)
(20)
The regenerator thermal ratios η’_{ REG} and η'' REG described in [1] are of considerable final importance are of considerable final importance. Schmidt and Willmott [1] have described the thermal ratios as follows: “The effectiveness of regenerator behavior is measured in terms of the thermal ratio η_{REG}. This is defined to be the ratio of the actual heat transfer rate to the thermodynamically limited maximum obtainable heat transfer rate in a counterflow regenerator of infinite heat transfer area. Could this maximum rate be achieved, the temperature of the gas leaving the regenerator in hot/cold period would be equal to the entrance temperature.”
Expression
(21)
(22)
can be derived. Here is the mean dimensionless temperature (with respect to position) of the solid:
(23)
We describe three example cases. For each case, we record similar calculation results.
The first example is a symmetrical case with the measures
(24)
The coarse, medium and fine meshes consisted of 4×8, 8×16 and 16×32 elements, respectively. Three total cycles were needed for “convergence” for each mesh density. The calculation results for the thermal ratios are given in Table 1.
Numerical solution in reference [1]  Numerical solution (4 × 8)  Numerical solution (8 × 16)  Numerical solution (16 × 32) 


η'_{REG}  0.494  0.49268  0.49332  0.49350 
η"_{REG}  0.494  0.49268  0.49332  0.49350 
Table 1: Thermal ratio η'_{REG} and η"_{REG}
It should be noted that the results from Schmidt and Willmott [1] (Table 5) were also obtained numerically. Temperature distributions in Ω are shown in Figure 4 and dimensional temperature distributions at τ=P' in Figure 5.
Effective convergence rate r and relative (percentage) error estimate er calculated for η’_{ REG} and η”_{ REG}are given in Table 2.
Hot period, η'_{REG}  Cold period, η"_{REG}  

r  1.79863  1.79862 
e_{r}  0.01874%  0.01874% 
Table 2: Effective convergence rate for thermal ratios and relative (percentage) error estimates
The second example case is also taken from Schmidt and Willmott [1]; Example 5.5. The dimensionless measures are
(25)
The coarse, medium and fine meshes consisted of 8×4, 16×8 and 32×16 elements, respectively. Twenty five total cycles were needed for each mesh density. The calculation results for the thermal ratios are given in Table 3.
Numerical solution in reference^{a}  Numerical solution (25 × 4)  Numerical solution (50 × 8)  Numerical solution (100 × 16)  

η'_{REG}  0.947  0.96111  0.94956  0.94760 
η"_{REG}  0.635  0.64519  0.63746  0.63616 
^{a}Schmidt and Willmott (1981)
Table 3: Thermal ratio η'_{REG} and η"_{REG}
Again, it should be noted that the results from (Schmidt and Willmott, (1981) were obtained numerically.
Temperature distributions in Ω are shown in Figure 6 and dimensional temperature distributions at τ=P' in Figure 7.
Effective convergence rate r and relative (percentage) error estimate er calculated for η’_{ REG} and η”_{ REG}are given in Table 4.
η'_{REG}  η"_{REG}  

r  2.55493  2.57525 
e_{r}  0.05317%  0.05138% 
Table 4: Effective convergence rate for thermal ratios and relative (percentage) error estimates.
For the third example case, we slightly modified the first case so that the ratio was:
Λ'⁄Λ''= Π'⁄ Π''= 0.83333 and
(26)
The coarse, medium and fine meshes consisted of 4×8, 8×16 and 16×32 elements, respectively. Four total cycles were needed for each mesh density. The calculation results for the thermal ratios are presented in Table 5.
Numerical solution (4 × 8)  Numerical solution (8 × 16)  Numerical solution (16 × 32)  

η'_{REG}  0.36210  0.36271  0.36279 
η"_{REG}  0.36210  0.36271  0.36279 
Table 5: Thermal ratio η'_{REG} and η"_{REG}
Theory demands here that η'_{ REG} = η'' _{REG}, which is seen to be satisfied by the numerical results. Temperature distributions in Ω are shown in Figure 8 and dimensional temperature distributions at τ=P' in Figure 9.
Effective convergence rate r and relative (percentage) error estimate er calculated for η’_{ REG} and η”_{ REG}are given in Table 6.
η'_{REG}  η"_{REG}  

r  2.84466  2.84835 
e_{r}  0.00469%  0.00466% 
Table 6: Effective convergence rate for thermal ratios and relative (percentage) error estimates.
As a general observation concerning Figures 4, 6 and 8, we may see the following. The dimensionless geometries in Figures 4 and 8 are rather similar and so are also the dimensionless temperature distributions. However, in the time direction narrow geometries of Figure 6 the temperature during the hot period “has not enough time” to penetrate far into the regenerator.
The results obtained in the three example cases presented above indicate that very good accuracy can be achieved with the finest meshes. The approach minimizes the effort of coding. It should not be difficult to extend the analysis to cases in which the inflow fluid temperature depends on time. This is achieved by simply changing the boundary conditions in the time direction appropriately. A more detailed presentation concerning the themes of the present article is given in Holopainen and Salonen [7].
The developed model can be applied to regenerators of both fixed and rotary heat storage mass. The model can be used as a tool in the design and optimization of heat exchangers when, for example, studying the effect of different hot and cold period lengths on outlet air temperatures and on thermal ratios.
The authors wish to thank Jorma Kinnunen, MSc (Tech.); Pasi Marttila, MSc (Tech.) and Jukka Tarvo, MSc (Tech.) for their advice on the use of COMSOL Multiphysics with MATLAB
Make the best use of Scientific Research and information from our 700 + peer reviewed, Open Access Journals