Modelling the behaviour of thermal energy harvesting devices with phase-change materials

This paper presents a new general theoretical model of thermal energy harvesting devices (TEHDs), which utilise phase-change materials (PCMs) for energy storage. The model's major goal is to identify a set of parameters under which these devices perform optimally, that is, attain the largest thermal buffering capacity and exchange heat with the surrounding phase as quickly as possible. For the first time, an expression for the characteristic harvesting time is developed from the constructal theory viewpoint under the optimal performance assumption, and a dimensionless criterion that characterizes PCM performance is provided. Furthermore, a new non-field solution of the energy equation governing the process of heat transfer within TEHDs with PCMs has also been derived. An expression for the effective thermal effusivity is then obtained. Finally, under a given set of boundary conditions and geometrical constraints, a novel simple technique for the optimal choice of PCMs in TEHDs has been established.


Indices
Sparing energy, more efficient transport and conservation of energy, as well as more efficient cooling of systems generating heat during their operation are among the most important tasks to be tackled by modern science and technology. As pointed out in the authors' previous work, phase-change materials (PCMs) are gaining a lot of attention as a sustainable approach to store energy at high ambient temperatures and release it at lower ones, thus buffering undesirable temperature oscillations 1 . Among the recent works, two excellent papers by Lin and coauthors provide the most comprehensive reviews on the state of the art in PCM applications and modelling 2,3 . As can be seen from these reviews, in spite of significant growth of the body of knowledge about PCM performance and applications, some lack of understanding is still present. In particular, no quantitative criteria are known for choosing PCM physical properties to provide their optimal performance. Diurnal temperature differences, peak power usage of electronic components, and solar heat exposure are some examples of transient energy sources that should be dissipated as rapidly as possible, or even better, saved for use when the temperature drops. These applications dictate the choice of the desired properties of PCMs, such as their high thermal buffering capacity and ability to exchange heat quickly with adjacent (matrix) materials. In view of the fact that the thermal buffering capacity cannot be represented by a single quantity, because it is defined by several integral parameters (e.g., latent heat and dimensions of the PCMs), the ability to exchange heat with a matrix material is best described by the thermal effusivity otherwise known as thermal inertia or thermal permeability 4 .
This paper presents a theoretical model of the behaviour of thermal energy harvesting devices (TEHDs), in which PCMs are used for energy storage. The main aim of the model is to establish a set of conditions, under which TEHDs with PCMs operate in an optimal way, that is, achieve the highest thermal buffering capacity and, at the same time, rapidly exchange heat with the adjacent phase (matrix material).
The study presented here originates from a practical need to develop a model of the TEHDs with PCMs in the course of designing a 300 kW heat accumulator for district heating networks. The accumulator is intended to be connected with the secondary system of the heat transfer station. The output medium from such an accumulator is heated water intended for supplying the house heating. The accumulator uses PCM substances to accumulate heat. This heat is transferred to the heat transfer medium by both convection and conduction through a designed metal interface. The input heat transfer medium to the accumulator is the returning water, which has fully transferred its heat within the transfer stations of individual houses. The purpose of heat accumulation in district heating systems is to significantly reduce the load on the heat source and distribution network at the time of maximum consumption and in an effort to use the production and transmission capacity of the entire system at the time of minimum consumption. The result is a significant reduction in the size of these heat sources and networks (e.g. boilers, valves, pipes, etc.), as well as a reduction in heat and hydraulic losses of these production and distribution systems.
The starting point of this study is the constructal theory 5,6 , which defines the frame for establishing a set of criteria for optimal performance of TEHDs with PCMs. In the next section, an expression for the characteristic harvesting time is derived under the optimal performance assumption from the constructal theory viewpoint.
"Problem formulation" of this paper presents the general theoretical formulation of the problem in question. In this section, for the first time. An expression to quantify the volumetric heat sink associated with the presence of the PCM is derived there as well.
In the beginning of "Solution procedure", for the first time, a dimensionless criterion, which characterises the PCM performance is introduced. Further in the same section, a new non-field solution of the energy equation governing the heat transfer process within TEHDs with PCMs has been obtained by the method of  . By this the method has been extended to solving optimization problems. From the solution thus obtained it follows that the thermal effusivity of the matrix material is one of the two parameters, which control the performance of TEHDs with PCMs. An expression for the effective thermal effusivity is then derived, which combines both the thermal effusivity of the matrix material and the dimensionless criterion of the PCM performance introduced earlier.
"Some practically important cases of the boundary conditions" presents the analysis of several important cases of the boundary conditions imposed on TEHDs with PCMs. In particular, the case of periodic or quasiperiodic boundary conditions is considered in more detail, because this case is important if energy harvesting is to be accomplished during diurnal cycles.
From the expression of the dimensionless criterion, which characterises the PCM performance, it follows that the volume fraction of PCMs plays an important role in the energy harvesting process. Hence, "Geometrical www.nature.com/scientificreports/ constraints on the PCM" is devoted to establishing an optimal value of the volume fraction of PCMs within the matrix material. In the same section, for the first time, the issue of the PCM particles distribution within the matrix material is discussed as well as the geometrical constraints related to the size of those particles. In addition, it has been shown how the thermal diffusing capacity of TEHDs with PCMs is related to the effective diffusivity of the constituent materials involved in TEHDs design. The last section of the paper provides the general discussion of the model developed in this study. The main aim of this section is to formulate a new concise procedure for the best choice of PCMs in TEHDs under a given set of the boundary conditions and geometrical constraints.

Optimal performance of thermal energy harvesting devices
In this study, criteria for optimal performance of thermal energy harvesting devices will be derived from the constructal law.
The constructal law is the law of physics that accounts for the phenomenon of evolution (configuration, form, design) throughout nature, inanimate flow systems and animate systems together. In its present form, the law was stated by Bejan in 1997 as follows 5,6 : For a finite-size system to persist in time (to live), it must evolve in such a way that it provides easier access to the imposed currents that flow through it.
Reformulated for the case of thermal energy harvesting devices, the constructal law states that, the optimal performance of these devices is achieved when a given amount of thermal energy is stored within the shortest time span possible. In practice, this means maximising the energy storage rate (power) consumed by the PCMs.
Thus, for a given amount of thermal energy to be harvested, Q hv , the optimal performance of the energy harvesting device is achieved if where t hv,min is the characteristic harvesting time.
The amount of thermal energy, which can be harvested, is given by a simple relation, namely: where m PCM is the total mass of the PCM used for the harvesting purpose and h f is the PCM enthalpy of fusion. From the latter equation, it may look like the task of maximising the amount of thermal energy harvested reduces to a trivial increase of the PCM total mass and choosing the PCM with the highest possible value of the enthalpy of fusion. It will be shown here, however, that geometric restrictions, imposed on the heat transfer process, define a constraint on the size of PCM particles and their relative position with respect to each other. This, in turn, sets the constraint on the total mass of PCM. Furthermore, the characteristic harvesting time is also related to the total mass of PCM.
Indeed, for a thermal energy harvesting device to perform in an optimal way, its operation cycle should exclude any idle periods. In other words, the amount of the thermal energy harvested must equal the amount of heat delivered to the PCM through the matrix material by conduction (thermal diffusion), that is, The latter expression is, of course, nothing else but the conservation of energy principle applied to the case in question.
The amount of heat delivered by thermal diffusion is given by where q ′′ is the heat flux through the surface of the PCM particle and A PCM is the total surface area of PCM available for energy harvesting. The heat flux in the latter equation can be estimated as where k is the thermal conductivity of the matrix material, T is the characteristic temperature difference, and ℓ denotes the average distance between the PCM particles within the matrix material. The characteristic temperature difference is the difference between the phase-change (e.g., melting) temperature and the ambient temperature. For Eq. (3) to remain valid during the entire operation cycle, the average distance between the PCM particles within the matrix material must be as close as possible to the characteristic length of thermal diffusion, that is, where α = k/(ρc p ) is the thermal diffusivity of the matrix material of the density ρ and specific heat capacity c p , respectively.
From combining Eqs. (2)-(6) and after some algebra follows the expression for the characteristic harvesting time  4 , the significance of which will be discussed in the end of this section. As can be easily seen from the latter equation, the characteristic harvesting time in general depends on the ratio of physical properties of the PCM and matrix material as well as the size of the PCM particles. It is also worth noting that the characteristic temperature difference, T , plays a crucial role as well. This temperature difference is defined from the boundary conditions imposed on the thermal energy harvesting device and the phase-change temperature of the PCM.
Observe that the result given by Eq. (7) is, in its mathematical form, practically identical to the result reported by Bejan 12 . Hence, the model developed in this study can, in principle, be easily extended to thermal energy harvesting devices operating in a periodic ("pulsating") mode.
As will be shown in the following sections, the thermal effusivity plays an important role in defining the performance of TEHDs. In fact, it is one of the two most important parameters, which control the behaviour of TEHDs. Hence, the physical meaning of the thermal effusivity is discussed here in more detail.
The thermal effusivity of a material is a measure of its ability to exchange thermal energy with its surroundings. The effusivity (thermal inertia) is defined as the square root of the product of the material's thermal conductivity and its volumetric heat capacity (see Eq. 8).
From Eq. (8), it is easy to notice that the units, in which the effusivity is measured, include s −1/2 . Hence, if the effusivity is to be a parameter in a certain equation, only its product with the square root of time may lead to physically meaningful quantities. This feature is further explored in the following sections.
Furthermore, it is also worth noticing here that, if two semi-infinite bodies initially at temperatures T 1 and T 2 are brought in perfect thermal contact, the temperature at the contact surface will be given by their relative effusivities: Equation (9) is valid for all times for semi-infinite bodies in perfect thermal contact. It is also a good first guess for the initial contact temperature for finite bodies.
From the heat transfer viewpoint, a finite-size body can be treated as semi-infinite (thermally thick) as long as its characteristic length scale, ℓ , remains less than or equal to the corresponding length scale of thermal diffusion, ℓ th , that is, where α denotes the thermal diffusivity of the material and t is time.
Because the thermal diffusivity is α = k/(ρc p ) , the time scale, during which a finite-size body can be treated as semi-infinite is Therefore, under the assumptions given by Eqs. (3) and (6), the heat transfer domain in this study can be treated as semi-infinite-at least, as long as the PCMs within a given TEHD operate near their phase-change temperatures.
To conclude this section, it is worth noting that, although Eq. (4) implicitly imposes a certain constraint on the size of PCM particles through their total surface area, the exact size constraint can be defined only after an assumption about the shape of PCM particles. The detailed analysis of this is presented in "Geometrical constraints on the PCM".

Problem formulation
The heat transfer domain composed of a matrix material, in which PCM is imbedded, is governed by where q ′′′ (r, t) models the effective volumetric heat sink associated with the presence of the PCM. The minus sign in front of the last term in the right side of Eq. (12) shows that heat is depleted by the PCM.
Note that Eq. (12) was obtained under the assumption that the Fourier law of heat conduction is valid 13 , namely that is, no ultra-fast heat transfer processes are involved and, hence, no thermal waves are present 14 . The domain of interest is initially in the state of thermal equilibrium, that is, T(r, 0) = T 0 = const . The two boundary conditions, which are to be imposed on Eq. (12), are not provided intentionally. This will become clear from the following sections.

Modelling the volumetric heat sink associated with the presence of the PCM. The total amount of heat, consumed by the PCM, is
where the amount of energy necessary to heat the PCM from its initial temperature, T 0 , to the fusion temperature, T f , equals whereas the amount of heat necessary to change the PCM phase equals where h f , is the enthalpy of fusion of the PCM. In Eqs. (15) and (16), m PCM = ρ PCM V PCM denotes the total mass of the PCM within the matrix material; ρ PCM and V PCM are the density of the PCM and its total volume, respectively.
Hence, the strength of the volumetric heat sink, associated with the presence of the PCM, can be quantified as where ϕ = V PCM /V denotes the volume fraction of the PCM. Finally, because the characteristic time of the heat transfer process is t hv = ℓ 2 /α from Eq. (6), the volumetric power of the PCM heat sink is In the following section a non-field solution of Eq. (12) for the boundary values of temperature and heat flux is to be derived 11 .
It is worth noting here that non-field solutions were not used in previous studies for analysis of optimization problems. Hence, the material, presented in the following section, is of methodological importance.
The initial condition, imposed on Eq. It is worth to remind here that the assumption of the domain's being semi-infinite remains valid as long as the PCM can be treated as an infinite heat sink-the latter is true if the PCM operates in an optimal way, that is, near its phase-change temperature.
For a semi-infinite domain, A(s) ≡ 0 , so that Eq. (22) reduces to (14)  Non-field solutions are known for that they relate local values of an intensive property-temperature in this case-and the corresponding flux. These solutions are valid everywhere within the domain including the boundaries. In view of the latter, non-field solutions proved to be extremely useful to transform boundary conditions of one kind into boundary conditions of another kind (e.g., the Dirichlet boundary conditions into the Neumann boundary conditions and vice versa) 11 .
In particular, the non-field solution, given by Eq. (27), can be inverted following the procedure adopted by Kulish and Lage 7 , that is, The non-field solutions, given by Eqs. (27)  Notice also that, choosing q where the second term in the right side is constant.

Some practically important cases of the boundary conditions
One of the simplest and yet practically important cases of the boundary conditions to be imposed on TEHDs is the case of a constant heat flux through the TEHD's boundary/surface, that is, q ′′ s = q ′′ 0 = const . In this case, Eq. (27) renders a straightforward result for the temperature variation on the boundary, namely In this case, however, the harvesting power of THEDs is constant with time and is fully defined by the geometry of the PCMs used, that is, and Observe that, not only both the auxiliary Fresnel integrals converge to zero quite rapidly ( z = 5 for practical purposes), but also their values are of the same order of magnitude. In view of this, because the difference between the auxiliary Fresnel integrals defines the last term in Eq. (37), the contribution of these integrals can be neglected for most practical applications. To illustrate the statement made here, Fig. 1 shows the difference fres(z) − gres(z).
Moreover, for practical purposes, it is sufficient to consider only major temperature variations during the operational cycle, such that [a n cos(ω n t) + b n sin(ω n t)] (36a) a n = 2 �t �t 0 T s (t)cos(ω n t)dt √ ω n a n cos ω n t − π 4 + b n sin ω n t − π 4 + www.nature.com/scientificreports/ where (�T) max is the maximal temperature variation in the course of a single operational cycle, the duration of which is t and ω = 2π/�t is the characteristic frequency of the cycle. The corresponding heat flux in this case is given by The total amount of thermal energy, harvested during a single cycle, is then As can be seen from the latter equation, for a given duration of the harvesting cycle, the total amount of thermal energy, which can be harvested during that cycle, solely depends on the physical properties of the matrix material and PCM as well as the geometric configuration of the PCM. Hence, it becomes important to define geometrical constraints, which define the PCM performance.

Geometrical constraints on the PCM
From Eq. (20) it follows that the PCM performance is directly proportional to the volume fraction of the PCM within the matrix material.
Consider a representative elementary volume (REV), which contains a single spherical PCM particle of diameter d-it is assumed here that the PCM particles are spherical. This particle is surrounded by the matrix material so that the total volume of the REV is where ℓ is the average distance between the PCM particles within the matrix material. If the ratio ℓ/d = β , then the volume fraction of the PCM is As follows from Eq. (43), the most compact packing of spherical PCM particles is achieved at β = 0 yielding ϕ max = π/6 ≈ 0.5236 . However, such an arrangement does not meet either heat transfer or geometrical constraints.
As has been pointed out in "Optimal performance of thermal energy harvesting devices", the optimal performance of the PCM is achieved if [see Eq. (6)] Combining the latter equation with Eq. (7) yields the relation between ℓ and d , namely: From the latter equation follows that the parameter β to be used in Eq. (43), in order to guarantee the optimal performance of the PCM, must be However, as have been demonstrated by Kulish and Lage 16 , if the PCM particles are located too close to each other, part of their surface area becomes shielded (blocked) by the neighbouring particles and, hence, excluded from the energy harvesting process. It has been found out that this happens when the volume fraction of the particles exceeds 16 per cent. Therefore, the maximal performance is achieved at φ = 0.16. The latter value has been determined from a numerical experiment with the relative error no greater than 0.1 per cent. At the same time, if φ < 0.16, part of the thermal energy, which could be stored within the PCM, would be conducted freely through the matrix material.

General discussion
The main purpose of this concluding section is to formulate a concise procedure for the best choice of PCMs in TEHDs under a given set of the boundary conditions and geometrical constraints. As pointed out in the preceding sections, the maximal amount of thermal energy is harvested if the entire surface area, which separates the PCM particles from the matrix material, is available for heat transfer. To achieve this, the volume fraction of the PCM is to be no more than 16 per cent. Then, it follows from Eq. (47) that the PCM particles are to be uniformly distributed within the matrix material, so that the distance between them equals almost one of their radius. Equation (46) provides a guideline for the choice of the PCM material with respect to the physical properties of the matrix material, namely: where T 0 is a certain reference temperature (e.g., the ambient temperature).
Once the PCM material is chosen, the PCM performance can be quantified by Eq. (20), after which the heat transfer process within TEHDs with PCMs can be fully modelled by Eqs. (27) and (28), respectively. Obviously, an appropriate set of boundary conditions, under which the TEHD with PCMs operates, is to be imposed beforehand.

Conclusions
In spite of significant growth of the body of knowledge about PCM performance and applications, some lack of understanding is still present 1-3 . In particular, no quantitative criteria are known for choosing PCM physical properties to provide their optimal performance.
In this study, for the first time: an expression for the characteristic harvesting time is derived under the optimal performance assumption from the constructal theory viewpoint and a dimensionless criterion, which characterises the performance of PCMs is introduced.
Furthermore, a non-field solution of the energy equation governing the heat transfer process within TEHDs with PCMs has been obtained. It is worth noting that non-field solutions were not used in previous studies for analysis of optimization problems. Hence, the material, presented in "Solution procedure" of this study, is of methodological importance-it extends the method of Kulish 10 to solving optimization problems with energy sinks/sources. Using the non-field solution, a new expression for the effective thermal effusivity is then derived.
To illustrate the usefulness of the non-field solution, several important cases of the boundary conditions imposed on TEHDs with PCMs have been provided. In particular, the case of periodic or quasi-periodic boundary conditions is considered in more detail, because this case is important if energy harvesting is to be accomplished during diurnal cycles.
Finally, a new concise procedure for the best choice of PCMs in TEHDs under a given set of the boundary conditions and geometrical constraints has been formulated.