An enthalpy based discrete thermal modelling framework for particulate systems with phase change materials

Powder Technology 354 (2019) 505–516

Contents lists available at ScienceDirect

Powder Technology journal homepage: www.elsevier.com/locate/powtec

An enthalpy based discrete thermal modelling framework for particulate systems with phase change materials T. Zhao, Y.T. Feng ⁎ Zienkiewicz Centre for Computational Engineering, Swansea University, UK

a r t i c l e

i n f o

Article history: Received 1 April 2019 Received in revised form 13 June 2019 Accepted 14 June 2019 Available online 19 June 2019 Keywords: Phase change material Discrete thermal element method Effective thermal conductivity Stefan problem

a b s t r a c t The latent thermal energy storage of phase change materials (PCM) is an attractive technique to use renewable energy. Systems with PCM capsules can be found in a wide variety of applications, but PCMs are usually approximated as a continuous phase in previous studies. The current work investigates this problem from the discontinuous point of view. The main objective is to develop an enthalpy based discrete thermal formulation that can take both heat conduction and phase change transition into consideration. The computational aspect of the formulation is fully discussed. The resulting algorithm is simple and effective. Its validity is demonstrated by solving a discrete/particle version of the one-phase Stenfan problem. In addition, the equivalent thermal properties of bulk particle materials with phase change are also derived based on a simple multi-scale modelling scheme. Numerical simulations are conducted to illustrate the effectiveness of the proposed enthalpy based discrete thermal modelling (DTEM) framework. © 2019 Elsevier B.V. All rights reserved.

1. Introduction The use of renewable energy, such as solar energy, has received much attention worldwide with the increase of energy shortage and environmental considerations. To reduce the time discrepancy between the energy production and the consumption, it is necessary to improve thermal energy storage (TES) techniques. Generally, solar heat energy can be stored in different ways including sensible heat, latent heat, or thermo-chemical reaction heat. In latent heat storage systems, the thermal energy is stored when the phase change material (PCM) undergoes a phase change, usually from solid to liquid. With the advantages of having a high density of energy storage and a narrow range of the operational temperature, the latent thermal energy storage is a particularly attractive technique [1–4]. The PCM needs to be contained properly for the successful utilisation in a latent thermal energy storage system. The most common type of PCM containment is macro-encapsulation in which a significant quantity of PCM is packaged in tubes, pouches, spheres, panels or other receptacles [5]. An undesirable property of PCM is its relatively low thermal conductivity. Therefore, it is needed to increase the surface/volume ratio to improve the heat transfer rate. The packed bed TES system is developed by packing a large number of small PCM particles prepared by micro-encapsulation techniques with the heat transfer fluid [6]. This packed bed TES system can be found in a variety of applications, such as

⁎ Corresponding author. E-mail address: [email protected] (Y.T. Feng).

https://doi.org/10.1016/j.powtec.2019.06.028 0032-5910/© 2019 Elsevier B.V. All rights reserved.

heating and cooling systems in buildings [1], solar thermal energy storage [2], solar cooling [3] and compressed air energy storage [4]. Since its wide applications, understanding the nature of a packed bed TES system has long been of a keen research interest. The main research methods used in previous investigations include the experiment [7,8] and numerical simulation. Considering technological and economical barriers of the experimental work, numerous contributions have been made in the field of computational modelling to study the performance of the TES system. These numerical models can mainly be divided into two groups: single-phase model and separate-phase model. In the single-phase model [9], the packed bed is regarded as a quasi-homogeneous medium; while in the separatephase model [10], both solid and fluid phases are considered separately with two energy conservation equations and are coupled by a heat exchange term. An extensive comprehensive review of different numerical models for PCM packed bed systems can be found in the review article [11]. The fixed grid method [12] and the adaptive mesh method [13] are adopted to solve the boundary value problems for partial differential equations where the phase boundary can evolve with time. The solutions of the differential equations are realised by finite-difference or finite element approximations [14]. It can be realised that the solid material (granular phase changing composites obtained by micro or macro-encapsulating PCMs) of the packed bed TES system is approximated as a continuous phase in the previous studies while it is actually the discontinuous particulate system. Therefore, it is necessary to investigate this problem from the discontinuous point of view.

506

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

Nomenclature α β Keff n x ΔT λ v νs ρ θ ε ξ a b Cp, C, c Capp E E∗ Fn H K k L mPCM N Q q R r r∗ r0 T V v

Contact angle Thermal expansion coefficient Effective thermal conductivity, W/(mK) Outward oriented unit boundary normal Relative coordinate of the particle Temperature difference, °C Freezing constant Poisson's ratio Solid volume fraction Density, kg/m3 Heating or cooling rate of the DSC measurement Thermal resistance per unit length, K/(mW) Phase transition parameter Contact radius, m Pipe length, m Specific heat capacity, J/(kgK) Apprent heat capacity, J/(kgK) Young's modulus, Pa Effective Young's modulus, Pa Normal force, N Enthalpy, W Contact conductance, W/K Thermal conductivity, W/(mK) Latent heat, W Mass of the phase change material, kg Number of particles Heat flux, W/s Heat flow, W Thermal resistance, K/W Current radius of particles, m Effective radius of particles, m Initial radius of particles, m Temperature, °C Particle volume, m3m 3 Heat diffusivity, m2/s

Since it was originated, the discrete element method (DEM) [15] has emerged as a reliable and effective numerical technique to model scientific and engineering problems involving particulate materials. Despite its extensive utilisations in the field of mechanical properties of dynamic particulate systems, DEM has also been used to simulate heat transfer in granular assemblies. Hunt [16] first introduced DEM to solve heat transfer problems of granular material flows and determined the effective thermal conductivity and self-diffusivity. Vargas et al. [17] developed the thermal particle dynamics (TPD) method based on DEM to study heat conduction in granular materials and stress effects on the conductivity of particulate beds [18]. Chaudhuri et al. [19] simulated heat transfer in granular flows in rotating vessels, in which granular flows and heat transport properties are taken simultaneously into account to understand their effects on dryer and calcination performance. A simple particle-particle heat transfer model is used in the above researches to account for the thermal effects and may inevitably suffer from a poor solution accuracy. Feng et al. [20] proposed the discrete thermal element method (DTEM) for accurate and effective modelling of heat conduction between individual particles in particulate systems in 2D cases. This method was extended to a simplified version, termed the pipe-network model [21], which significantly simplifies the solution procedure of the original DTEM. These studies and others have demonstrated the capability of the DEM as a powerful tool for investigating

heat transfer problems in particulate systems [22,23] and for determining the corresponding effective thermal conductivity [24–26]. However, little effort has been made to investigate particulate systems with phase change materials. The main objective of this paper is thus to develop an enthalpy based discrete thermal formulation that can take both thermal conductivity and phase change into consideration in particle systems. The equivalent thermal properties of bulk particle materials with phase change will also be derived based on a simple multiscale modelling scheme. This paper is organised as follows. In Section 2, a brief description of the PCM together with its thermal material characterisation and mathematical modelling is provided. The enthalpy based DTEM framework is developed in Section 3. The computational aspect of the formulations for both heat conduction modelling and phase change modelling are fully discussed and described. In Section 4, the way to determine the effective thermal properties are addressed in detail. The validity of the DTEM method is demonstrated in Section 5 by solving a particle version of the classic one-phase Stenfan melting problem [27]. Section 6 is devoted to conduct numerical simulations which illustrate the effectiveness of the proposed method to model the particulate system of PCM capsules. The conclusion is drawn in Section 7. 2. Phase change material Thermal energy storage is an important technology to balance the supply and demand of energy. It has been widely used in applications such as building, concentrated solar plants and thermal management of batteries [28], to name a few. A TES system can be a sensible, latent or chemical heat storage as shown in Fig. 1 (a). The sensible heat storage applies a temperature gradient to the material to store or release heat. The disadvantage of the sensible heat storage is the low storage capacity which leads to a huge volume of the system. The latent heat storage can provide a high storage density with a small temperature difference between storing and releasing heat. The materials used for the latent heat storage are called phase change materials and are characterised by storing large amount of thermal energy while changing from one phase to other (usually solid-liquid states) at a narrow temperature range which has been regarded as a constant temperature in some research, and presenting high heat of phase change (latent heat) [29]. The enthalpy variation with temperature of both sensible and latent materials are shown in Fig. 1(b). The enthalpy-temperature relation of real materials in nature for all solid, liquid or phase change mixtures is non-linear as shown in the black line, but it is often simplified as a (piecewise) linear relation (blue line) for convenience in engineering applications. PCMs can be classified as organic and inorganic, of which the organic material paraffin is the most popular commercial PCM due to its long-term stability, suitable phase change temperature and acceptable price. 2.1. Encapsulation of PCM Main drawbacks of PCMs include low thermal conductivity, leakage, subcooling and flammability. Appropriate methods need to be developed to ensure the successful utilisation of PCMs. Containment methods can be classified as the bulk storage in tank heat exchangers, the macro or micro encapsulation. In the encapsulation process, individual particles or droplets of solids or liquid PCMs (the core) are surrounded or coated with a film of polymeric materials (the shell). This method has the advantages of providing a large heat transfer area and reducing the reactivity of the materials with the outside environment. The characteristic diameter of such granular phase change composites ranges from a few

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

507

Fig. 1. Energy storage materials.

millimeters to a few centimeters. Granular phase change materials (GPCM) offer the advantage of maintaining their macroscopic solid form during phase change which is the solid component of the latent heat thermal energy storage system (LHTES) [30]. 2.2. Material characterisation The phase change of most PCMs occurs in a temperature range rather than at a constant temperature. From heat flow results of differential scanning calorimeter (DSC) measurements, the apparent heat capacity of granular PCMs is acquired as [12]. C app ðT Þ ¼

qðT Þ þ Cp mPCM θ

ð1Þ

where mPCM is the mass of the PCM used in the DSC, θ is the heating or cooling rate of the DSC measurement and Cp is the specific heat of the solid or liquid PCM. Packed-bed column experiments during charging and discharging modes are usually carried out to analyse the granular phase change material [5]. 2.3. Mathematical models The complex transient nature of granular PCM systems and the high cost of the set-up for an experimental testing makes it necessary to use numerical models to investigate PCM systems. Existing numerical models developed can mainly be divided into two groups: single or two phase model [31]. In the single phase model, the solid phase (PCM spheres) and the fluid phase are considered as one phase. Ismail and Stuginsky [32] demonstrated that this model is useful in analysing fixed beds of both high thermal conductivity and thermal capacity in comparison to the working fluid. In the two phase model, solid and fluid phases are considered separately, and the model can be classified into three (sub-)categories: the concentric dispersion model [3], the continuous solid phase model [33], and the Schumann model [10]. All the four models mentioned above assume the solid phase as a porous material, not as a medium comprised of individual particles. Thus, the accuracy of these models depends on the effective thermal conductivity of the porous material and the total heat transfer coefficient between the fluid and the solid, both of the parameters are usually determined from empirical correlations.

view. Therefore, an enthalpy based discrete thermal modelling framework is developed in the current work. It should be noted that only the solid component (PCM spheres) is considered as the main object of this work is to propose a novel method to treat the heat conduction and phase change transition in the particulate system with PCMs rather than model a specific packed-bed column experiment. 3.1. Heat conduction modelling Different mechanisms exist in heat transfer of particulate systems, including thermal conduction through the solid; thermal conduction through the contact area between two particles; thermal conduction to the interstitial fluid; heat transfer by fluid convection; and radiant heat transfer between the surfaces of particles [34]. According to Batchelor and O'Brien [35], when the interstitial medium is stagnant with a smaller thermal conductivity compared to that of the particles, thermal conduction through the contact area plays the dominate role in heat transfer process. Thus, the current work of heat transfer in particulate systems with phase change material is focused on contact conductance between particles. 3.1.1. Thermal DEM To simulate the heat transfer process in a particle system, temperature is introduced into the DEM as an additional degree of freedom. For the particulate system with a large number of particles, it is extremely time consuming to determine the temperature distribution of each particle. Thus, each particle has been assumed as isothermal in DEM simulations. Heat flow occurs via conduction in the active contact area that connects two particles (i and j) whose (average) temperatures are Ti and Tj respectively. The Fourier's law of heat conduction denotes the amount of heat transported across their mutual (contact) boundary per unit time as
 Q ij ¼ K ij T j −T i

ð2Þ

where Kij, termed the contact conductance, is the amount of heat transported per unit temperature difference per unit time. The evolution of the temperature of particle i is given as dT i Q i ¼ dt Ci

ð3Þ

3. The enthalpy based DTEM framework

where Qi is the total amount of heat transported to particle i from its neighbor particles calculated from Eq. (2) as

Considering the limitations of the continuous based models and the discontinuous nature of granular PCMs, it will be a new direction to simulate the particulate system with PCMs from the discontinuous point of

Qi ¼

N X j¼1

Q ij

ð4Þ

508

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

in which N is the number of particles in contact with particle i; Ci is the total heat capacity of particle i which is defined as C i ¼ ρi c i V i

ð5Þ

where ρi, ci and Vi are the density, specific heat capacity and volume of particle i respectively. Eq. (3) can be discretised by expressing the time derivative using, for instance, the forward finite difference to update the particle temperature as T ihtþΔti ¼ T iht i þ

Qi Δt Ci

  3 F n r  1=3 K ij ¼ 2k ¼ 2ka  4E

ð8Þ

where k is the thermal conductivity, Fn is the normal force acting between the particles, a is the contact radius, r ∗ = ri ∗ rj/(ri + rj) is the effective particle radius, and E ∗ is the effective Young's modulus for the two particles 2

ð6Þ

There are two criteria that need to be satisfied in this thermal DEM method: (1) The temperature of particle i keeps the same value at different contact points with all the neighbor particles. This requires that the resistance to heat transfer inside the particle is significantly smaller than the resistance between the contacting particles, which is generally true because the contact area is often far smaller than the size of particles; (2) The stability criterion of the explicit time integration scheme requires that the temperature change of each particle sufficiently slow so that thermal disturbances do not propagate further than its immediate neighbors during one time-step. This can be satisfied by choosing a sufficiently small time-step. Another aspect considered in modelling of thermal conduction is the thermal expansion and contraction of the particle size. Similar to the previous treatment [36–38], the radius expansion/contraction is considered as r ¼ r 0 ð1 þ βΔT Þ

process in thermal DEM as

ð7Þ

where r0 is the initial radius of particles at the reference temperature T0, r is the current radius, ΔT = T − T0 is the temperature change, and β is the thermal expansion coefficient (assuming constant in the current work). When the thermal expansion of particles is taken into account, granular materials would exhibit the settling behavior under thermal excitations which needs to be simulated by considering kinematic of particles in conjunction with the thermal modelling. This has two consequences. One is the increase of the packing density and the other is the particle rearrangement. The packing density has a major influence on thermal properties of granular phase change materials as can be seen in Section 6.2, while the effect of particle rearrangement on thermal properties is less clear for randomly packed particle configurations. Nevertheless, the effect of particle size change may not be significant in the thermal modelling of phase change materials as the narrow range of the operational temperature during the charging and discharging processes and the randomness of particle arrangement under thermal excitations. Thus most of the thermal simulations conducted in the current work, unless those mentioned specifically, do not consider the size expansion to reduce the computation time because the required timestep in the thermal DEM simulation is often orders of magnitude larger than the corresponding time-step in the mechanical DEM [24]. The thermal induced size change effect will be discussed in detail in Section 6.4. 3.1.2. Thermal contact model To accomplish the above thermal DEM simulation, the contact conductance Kij in Eq. (2) needs to be determined properly. Contact conductance refers to the heat transmit ability of two touching materials. Most work related to this topic can be found in fields of microelectronics, aircraft industry, nuclear industry and nano-technologies [39]. The thermal contact model in thermal DEM modelling borrows some mature results from the early work of contact conductance [40,41]. The approximate analytical solution of the contact conductance between two smooth, elastic spheres is often adopted to consider the heat transfer

1−ν 2i 1−ν j 1 þ  ¼ E Ei Ej

ð9Þ

in which νi and νj are the Poisson's ratios of the two particles. Despite the above relation (8) that is often used in most of thermal DEM simulations [17,18,24,42], a thermal pipe contact model is developed in Particle Flow Code [43] in which heat flows between particles through a one-dimensional pipe where the contact conductance is defined as K ij ¼

1 εb

ð10Þ

in which ε is the thermal resistance per unit length, and b is the pipe length. The determination of the contact conductance in the above models is mainly done in an ad hoc manner rather than based on a rigorous theoretical foundation. Feng et al. [20,21] proposed a 2D pipe-network model for the modelling of heat conduction in particulate systems which provides a more rational and accurate model to represent heat conduction. In this 2D pipe-network model, one particle can be conceptually represented by a simple star-shaped pipe network model in which the centre is connected to each contact zone as shown in Fig. 2. For each individual pipe j of the ith particle, the corresponding thermal resistance Ri, j is given by

Ri; j

1 ¼ πki

2 α 4i; j α 6i; j 3 α i; j −lnα i; j þ þ þ þ 2 36 2700 79380

! ð11Þ

where αi, j is the contact angle of the contact region defined as the ratio of the contact length and the particle radius. Then, the thermal contact conductance Kij between the ith and jth particles is formulated as K ij ¼

1 Ri; j þ R j;i

ð12Þ

A similar approach can also be developed to define the thermal contact conductance for contacting spheres. 3.2. Phase change modelling Phase change materials are viewed as latent heat storage units as a large amount of energy is absorbed or released when materials changes from solid to liquid and vice versa. In the DTEM method, the phase change transition is taken into account by using the enthalpy concept. The change in enthalpy equals to the heat absorbed or released for the enclosed system at a constant pressure. The quantity of thermal energy possibly stored depends on the enthalpy variation in the working temperature range. In the current DTEM framework, a PCM sphere is regarded as a phase change material as shown in Fig. 1, and the phase change process occurs at the temperature Ts and ends at the temperature Tl as shown in Fig. 1 (b), and the enthalpy variation is equal to the latent heat L. The enthalpy

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

509

Fig. 2. Discrete thermal element and pipe-network model for disc [20,21].

of each sphere is Hi and can be obtained by H ihtþΔt i ¼ Hihti þ Q i Δt

4.1. Effective thermal conductivity ð13Þ

For the phase change transition from solid to liquid, the onset enthalpy is Hs = CiTs and the end enthalpy is Hl = Hs + L as shown in Fig. 1(b). The specific capacity ci and thermal conductivity ki of each particle need to be modified according to the enthalpy Hi of that particle as follows 0

Cs; H i ≤Hs H i −Hs B ðC l −C s Þ; H s bHi bHl ci ¼ @ C s þ Hl −H s Cl; H i ≥Hl

ð14Þ

and 0

ks ; H i ≤H s H i −H s B ðkl −ks Þ; H s bH i bH l ki ¼ @ ks þ H l −Hs kl ; H i ≥H l

ð15Þ

To track the phase change evolution of a particulate system, the phase change time Tic of each particle is recorded. A phase transition is assumed to occur at the time when a particle has stored or released a certain amount of the latent heat, or the phase transition latent heat Hpt(ξ) determined by a given parameter ξ ∈ [0,1] as H pt ðξÞ ¼ ð1−ξÞH s þ ξH l

ð16Þ

With the phase change time of each particle is known, the position evolution of the phase change front in a particulate system can be obtained in DTEM modelling. 4. Effective thermal properties In design of a LHTES, it is necessary to properly determine thermophysical properties including the effective thermal conductivity and other thermal properties of the system. For the particulate system of phase change material, it is not easy to characterise granular phase change composites due to the nature of the material and the heterogeneity of the sample. Small and non-uniform diameters of spherical capsules also complicate the evolution of latent heat in a heat storage system. The current DTEM method makes it possible to numerically obtain the effective thermal properties of a granular PCM system by applying a homogenisation based technique.

To predict the effective thermal conductivity (ETC), Keff, of granular materials is an important subject for many scientific and industrial applications involving particulate systems. Various experimental tests and analytical models have been proposed to measure or analyse the ETC of granular materials [44–46]. Compared to these methods, numerical simulations can be employed to obtain heat transfer of granular materials, which makes the numerical simulation using DTEM a powerful means to predict the ETC. In theory, the ETC of granular assemblies can be determined by solving the equation of Fourier's law of heat conduction. Therefore, the heat flux q″ and temperature gradient ∇T need to be obtained from DTEM simulations. A simple way is to obtain them from a boundary value problem where the granular sample is considered as a continuum phase. A more accurate method is to obtain the average heat flux 〈q″〉 and temperature gradient 〈∇T〉 of a granular assembly [24], which is analogous to compute the average stress and strain of a particulate system by using the multi-scale method [47–49]. In granular materials, the micro scale is taken to be a scale where individual particle responses are measurable, while the macro scale is where the whole assembly is considered as a continuum phase. The following description of using the average heat flux and temperature gradient to predict the ETC of a particle assembly is mainly adopted from [24]. 4.1.1. Average temperature gradient Consider the steady state heat transfer of a granular system, Ω, i where the approximate temperature T~ of each particle is expressed as i T~ ¼ ∇T  xi

ð17Þ

in which xi is the relative coordinate of the particle, and ∇T is the temperature gradient at the macroscale level where the granular assembly is considered as a continuum phase. Because of the fluctuation of the particle temperature distribution in i the particulate system, the real temperature T^ is not exactly match the i approximate temperature T~ , so

i T^ −∇T  xi ≠0

ð18Þ

To find the specific 〈∇T〉, the square sum of the deviations in Eq. (18) should be minimised. Thus by employing the least-squares approach,

510

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

〈∇T〉 can be obtained as a solution of the following linear equations X

 i  xi T^ −h∇T i  xi ¼ 0

ð19Þ

i∈Ω

4.1.2. Average heat flux In the state of heat equilibrium without heat sources, the temperature of each particle remains constant so the total heat exchange through the boundary of the granular assembly, ∂Ω, is zero: Z ∂Ω

X

q″  ndA ¼

qib ¼ 0

ð20Þ

i∈∂Ω

where q″ is the heat flux flowing through the boundaries, n is the outward oriented unit boundary normal and ∫∂Ω is the integration on the boundary of the particulate system, qib is the total heat flowing through the boundary contacting with particle i. The average heat flux 〈q″〉 of the particulate system is given by  ″ 1 q ¼ V

Z Ω

q″ dV ¼

1 X ib ib q x V

ð21Þ

i∈∂Ω

where V = ∣ Ω∣ is the volume of the particulate system and xib denotes the coordinates of the contact point between particle i and the boundary. Hence, the average heat flux of the particulate system can be computed from the simulation. 4.1.3. Effective thermal conductivity By applying Fourier's law of heat conduction, the effective thermal conductivity Keff should satisfy  Keff h∇T i ¼ − q″

ð22Þ

Keff, as a tensor, is required to be symmetric and positive definite which are automatically satisfied by enforcing the symmetry of Keff [50]. Then Keff is obtained from solving the following overdetermined linear system of equations 21

T ;1

6 6 6 6 62 6 T ;1 6 6 6 6 6 63 6 T ;1 6 6 6 6 6 6 6 4

1

2

3

T ;2

T ;2

T ;2

1

2

3

T ;3

T ;3

T ;3

2

3

T ;1

T ;1

T ;1

1

2

3

T ;2

T ;2

T ;2

1

2

3

T ;3

T ;3

T ;3

−1

1

T ;1

1

T ;2

2

T ;1

2

T ;2

3

T ;1

3

T ;2

−1

1

8 1 ″ 9 q > > >  ″ 1 > > > > 1 > > q 2> > > > >  > > ″ 1 > > > > q > > >  3 > > > 2 ″ > > > q 1> > > > 2 ″ > > > > > > > q > = < 2  ″ 2 > q ¼−  ″ 3 > > 3 > > > > > q 1 > > > > > ″ 3 > > q > > >  2 > > > ″ 3 > > > q 3> > > > > > > > > 0 > > > > > > > > 0 > > ; : 0

1

1

4.2. Effective latent heat The effective latent heat Leff of the particulate system of PCMs is an important property which indicates its heat storage capacity. For a granular system, the effective latent heat Leff can be simply obtained as Leff ¼ Lν s

ð24Þ

where L is the latent heat of the granular particles, and νs is the solid volume fraction, or packing density, of the particulate system. 5. Validation As our current work is to treat granular PCMs from the discontinuous point of view, its validity will be demonstrated by applying it to solve a discrete/particle version of the classical phase change problem – the Stefan problem, where a particulate system with PCMs will be treated as a one-dimension continuous medium. The interface evolution of the particulate system should be analogous to its continuum counterpart. For the continuum phase, the phase transition is described as a boundary value problem of partial differential equations (PDE), where the phase boundary moves with time. In the following subsections, the moving boundary Stefan problem is introduced first. The latent heat is an important parameter which affects the movement of the phase boundary (solid-liquid interface) of the system. The validity of the DTEM method will be demonstrated by comparing the actual effective latent heat Leff of the particulate system with the fitted one Leff‐fit derived from the numerical results. Note that no thermal expansion/contraction is considered in this section. The thermal induced particle size change effect will be discussed in Section 6.4. 5.1. The stefan problem This classical moving boundary problem dates back to the study of the melting of glaciers by the Slovenian physicist Jozef Stefan in 1889 [27]. It aims to describe the temperature distribution in a homogeneous medium undergoing a phase change. The physical constraints of this problem are the conservation of energy and the local velocity of the interface depends on the heat flux discontinuity at the interface. From a mathematical point of view, the phases are regions in which the solutions of the underlying PDE are continuous. The moving boundary are infinitesimally thin and the PDE are not valid at phase change interfaces. Therefore, an additional Stefan Condition is needed to obtain closure. The analytical solutions for this problem is available for onedimensional cases of an infinite region with simple initial and boundary conditions.

3 1

Cartesian coordinate system. This equation is solved as a linear leastsquares problem.

78 eq 9 7> k > 7> 11 > T ;3 7> > > keq > 7> > > 12 > 7> > > eq > > 7> > k > 13 > > 7> > > keq > 7> = < 21 > 2 eq T ;3 7 7 k 7> 22 > eq > 7> > > k23 > 7> > > 7> eq > > k31 > > > 7 > > 3 > eq > > T ;3 7 > > k32 > 7> > > 7> ; 7: keq 33 5 1

−1

ð23Þ

where all blank coefficients of the matrix are 0 and the superscripts 1, 2 and 3 represent the boundary problems in the three directions of the

5.1.1. Formulations For a semi-infinite region, the initial temperature is T0. At time t N 0, the temperature of the boundary at x = 0 is suddenly kept at Tw. The governing equations for the temperature distribution T(x, t) in the region are formulated as follows 0

  ∂T ∂ ∂T @ C l ∂t ¼ ∂x kl ∂x T l ð0; t Þ ¼ T w

0bxbsðt Þ;

tN0

ð25Þ

tN0

0

  ∂T ∂ ∂T sðt Þbxb þ ∞; @ C s ∂t ¼ ∂x ks ∂x T s ðþ∞; t Þ ¼ T 0

tN0 tN0

ð26Þ

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

latent heat Leff‐fit is derived based on the simulation results as

subject to the initial condition ð27Þ

T ðx; 0Þ ¼ T 0

where s(t) is the interface position, and the subscripts s and l indicate the solid and liquid phases. The moving rate of the interfaces is controlled by the latent heat lost or absorbed at the boundary. The following equation, known as the Stefan Condition, describes this process. kl

∂T l ∂T s ds −ks ¼ −Λ dt ∂x ∂x

at x ¼ sðt Þ

ð28Þ

where Λ = ρL is the heat of phase change per unit volume (note that L is the latent heat coefficient) and ds/dt is the velocity of this interface. Moreover, the temperature verifies: Tl ¼ Ts ¼ Tm

at x ¼ sðt Þ

ð29Þ

where Tm is the melting temperature. This problem can be formulated in non-dimensional variables for a finite sheet 0 ≤ x ≤ l (l is a standard length) by assuming constant thermal values and using a simple scaling [51]. The non-dimensional form makes it convenient to obtain the analytical solution and save computational costs of the numerical solutions. 5.1.2. Analytical solution The exact solution of the one-dimensional Stefan problem, available in the literature [52], can be derived by using the similarity variable to transform the governing equation from partial derivatives to an ordinary differential equation. The solution of the temperature distribution is T ðx; t Þ ¼ T 0 −

 T0 x erf pffiffi erf ðλÞ 2 t

ð30Þ

where λ is the freezing constant given by the root of the following transcendental equation pffiffiffi 2 λeλ erf ðλÞ ¼ Ste= π

ð31Þ

in which erf is the error function and Ste is the Stefan number defined as Ste ¼

CΔT L

511

ð32Þ

Leff‐fit ¼

CΔT pffiffiffi 2 λeλ erf ðλÞ π

ð34Þ

If the proposed enthalpy based DTEM method for simulations of problems with PCMs is valid, the Leff‐fit obtained by Eq. (34) should be close to Leff calculated by Eq. (24). Numerical simulations using the enthalpy based DTEM method are conducted to derive the fitted effective latent heat Leff‐fit of a particulate system with granular PCMs. Three cuboid shaped boxes with 1 m in both length and width and 5 m in height filled with PCM capsules have been randomly generated with some initial contacts between neighbouring particles. These samples have different packing densities νs = 0.6,0.65 and 0.7 respectively. The diameters of the capsules have a uniform distribution between 0.05 m to 0.06 m. The top and bottom of each box have a constant temperature of 100 °C and 60 °C respectively. Heat is conducted between the top and bottom of the box with the contacting particles in the sample. The side surfaces of the boxes are adiabatic or insulated. The initial temperature of all the PCM capsules is 60 °C. It is assumed that the particles are ideal phase change materials and the melting process occurs at a temperature point Tm = 60 °C rather than a temperature range. The phase transition latent heat Hpt is set by choosing the parameter ξ = 0.5. A particle is assumed melted when the enthalpy of the particle H ≥ Hpt. The thermal properties of the PCM capsules are set as Cp = 430 J ⋅ kg−1 ⋅ K−1, ks = 4 W ⋅ m−1 ⋅ K−1, kl = 8 W ⋅ m−1 ⋅ K−1, and L = 10 kJ/kg. The thermal simulation is performed until a steady-state heat transfer is achieved for the particulate system. The particle system is then treated as a one-dimensional continuous medium along the vertical (z) direction. The corresponding equivalent thermal properties of the system are obtained. Fig. 3 shows the temperature distribution of Sample 1 (νs = 0.6). During the initial heating period, the capsules near the top are heated up while the remaining part keep the initial temperature. As time elapses, the temperature rises gradually in the system. The phase change times and positions of all the particles in Sample 1 are depicted in Fig. 4, and their relation is curve fitted using Eq. (33), from which the freezing constant λ is obtained as 0.0812. Then, by applying Eq. (24) and Eq. (34), the effective latent heat Leff and the fitted one Leff‐fit can be computed. Both values of Leff and Leff‐fit for all the three samples are listed in Table 1. It can be seen that the fitted values of the effective latent heat Leff match with the actual values of the

where C is the specific heat of solid phase in the freezing process, but is the specific heat of the liquid phase in the melting process; ΔT is the temperature difference between the two phases. The moving front position is expressed as pffiffiffiffiffi sðt Þ ¼ 2λ vt

ð33Þ

where v = ks/Csρ is the heat diffusivity. 5.2. Effective latent heat of a particle system based on DTEM As mentioned in Section 3.2, the phase change time of each particle can be recorded during the simulation. The time history of the phase change times and positions of all the particles in the system will depict how the melting interface of the system evolves with time. If the particulate system is regarded as an equivalent continuum, it should have a similar relation as described by Eq. (33). The phase change time - position points of the particles can be curve fitted by this relation, and then the resulting freezing constant λ can be used to derive the effective latent heat based on both Eq. (31) and Eq. (32). Thus, the fitted effective

Fig. 3. Temperature distributions and time evolutions of the particulate system with PCM capsules.

512

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

6.1. Temperature evolution

Fig. 4. Determination of the effective latent heat (Sample 1).

Table 1 Effective latent heat of the particulate system. Sample

1

2

3

Packing density νs Freezing parameter λ Leff (kJ/Kg) - Eq. (24) Leff‐fit (kJ/Kg) - Eq. (34)

0.6 0.0812 6.0 6.027

0.65 0.0789 6.5 6.404

0.7 0.0761 7.0 6.880

Another thermal simulation is conducted using Sample 1 where the granular capsules are considered as sensible material by setting the latent heat L to be 0. Fig. 5(a) shows the temperature profile of Sample 1 at different heights for various time instances. A similar profile for the phase change material used in the previous simulations is also displayed in Fig. 5(b) for comparison. The top (z = 0 m) and bottom (z = 5 m) boundaries of Sample 1 have a constant temperature of 100 °C and 60 °C respectively. The initial temperature of all the capsules in the sample is 60 °C. The heat flux flows from the top to the bottom through the granular capsules. The temperature of the capsules increases during the heat-absorbing process. Heat conduction between two capsules occurs because of their temperature difference. The phase change material has the ability to absorb large amount of heat without changing its temperature. Thus it takes longer for heat to transfer from the top to the bottom of the phase change materials in the sample. Fig. 5 (a) shows that at the time 600 s, the heat transfer process is almost completed in the whole sample for the sensible materials; while Fig. 5(b) indicates that only just over 1/10 of the sample o phase change materials is affected by the heat transfer at the same time. Also, the sample with sensible material reaches the steady-state temperature distribution in about 800 s, while for the system with phase change materials, due to the notable heat storage capacity, the same process takes more than 18,000 s, which is more than 20 times that of the sensible one. This obvious significance can also be observed from the temperature evolution of the particles at different heights as shown in Fig. 6. Capsules at five different heights are selected and their temperatures are depicted against time. For the capsules of sensible material, the temperature increases linearly and reaches to the final temperature very quickly. 6.2. Effect of capsule size and packing density

effective latent heat very well, thus demonstrating the validity of the DTEM method.

6. Further illustrations Further simulations are performed to show the effectiveness of the DTEM method in modelling heat transfer in particulate systems involving phase change materials.

Additional two samples with PCM capsules with different radius ranges but a constant packing density νs = 0.6 are generated. The resulting three samples are termed as small (0.04–0.05 m), middle (0.05–0.06 m) and large (0.06–0.07 m). Fig. 7(a) show the total melt fraction of the particulate system with different capsule sizes. It can be seen that the samples with smaller capsules taking less time for charging which indicates the same phenomenon as pointed in the previous research [53], but the effect is weak here. This difference is caused by

Fig. 5. Temperature profiles and evolutions of the particulate system with two different materials at different time instances.

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

513

Fig. 6. Temperature evolutions of the particulate system with two different materials at different heights.

the dissimilarity of the packing configuration of the samples with different particle sizes. The detailed explanation can be found in the work of [54,55]. Fig. 7(b) shows the evolutions of the total melt fraction of the three samples with different packing densities. The difference between these cases is obvious. The sample with the largest packing density attains the complete melting position quickest. The time needed for the loosest sample to reach the final state is more than three times than that of the densest sample. This is because that a denser sample has a larger effective thermal conductivity which can be seen in Table 2 below. 6.3. Determination of effective thermal conductivity As described in Section 4, the effective thermal properties can be calculated from the numerical results. Simulations are conducted on five different samples with different particle size distributions or packing

densities as shown in Table 2. In the current simulations, the temperatures are specified at the top and bottom boundaries of the samples, and all the other boundaries are isolated. So only the thermal conductivity along the z direction is determined. Take Sample 1 as an example. The heat flux at the state of heat equilibrium is obtained when the heat flux at the top and bottom of the sample are equal. Thus, the average heat flux 〈q″〉 = 374.44 Wm−2 is shown in Fig. 8(b). The average temperature gradient is determined by fitting the temperature distribution using a linear relation which leads to 〈∇T〉 = 6.92 K/m as shown in Fig. 8(a). Then, by using Eq. (22), the effective thermal conductivity is computed as shown in Table 2. The effective thermal conductivity of the other samples are calculated by the same method and listed in Table 2. It can be observed that the effective thermal conductivity is affected significantly by the packing density of a sample while the particle size distribution has a minor effect on the conductivity.

Fig. 7. Melt fraction evolutions of the particulate system with different particle sizes and packing densities.

514

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

Table 2 Effective thermal conductivity of the particulate system. Sample

1

2

3

4

5

Packing density νs Radius r/m Keff (W/mK)

0.6 0.05–0.06 54.08

0.65 0.05–0.06 137.95

0.7 0.05–0.06 236.66

0.6 0.04–0.05 55.85

0.6 0.06–0.07 48.61

6.4. Effects of thermal induced particle size and packing configuration changes The size change of particles due to thermal expansion will cause the increase of packing density and the particle rearrangement of the granular system. The consequence of particle rearrangement induced by the

thermal excitation may be similar to that caused by difference packing configurations. Three particle packings are generated using the same parameters as in Sample 2. The only difference between these packings is that they are generated using different random seeds. The same thermal simulations are conducted on these three packings with or without considering the thermal induced size expansion. The particle size expansion is considered following Eq. (7). The thermal expansion coefficient is set to be β = 10−4K−1 based on [56]. The curve fitted finial steady-state temperature distributions and phase change times of these six cases are shown in Fig. 9(a) and (b) respectively. Comparing the phase change times of the same packing with or without size change in Fig. 9(b), it can be seen that the heat transfer is

Fig. 8. Steady-state particle temperature distribution (a) and heat flux evolutions at the boundaries (b) for Sample 1.

Fig. 9. Final temperature distribution (a) and phase change time (b) for different packings with or without size change effect.

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

slightly quicker when the size expansion is taken into account in the simulation, as expected. However, the difference is not significant and comparable with the difference caused by different random packing configurations with the same packing density. Therefore, to reduce the computation time and improve the simulation efficiency, it is a reasonable choice to conduct the thermal modelling without considering the size expansion of particles. 7. Conclusion This work has developed an enthalpy based discrete thermal modelling framework for particulate systems with phase change materials which can consider both the heat conduction process and the phase change transition. The proposed algorithm is numerically very simple and effective. In addition, the equivalent thermal properties of bulk particle materials with phase change have also been derived based on a simple multi-scale modelling scheme. The proposed methodology is assessed by solving a particle version of the classic one-phase Stenfan melting problem. Additional numerical simulations have also been conducted to illustrate the effectiveness of this modelling framework. The effect of capsule size and packing density can be evaluated directly in the current DTEM method. A particle system of phase change materials with smaller capsules or larger packing densities takes lesser time for charging or discharging. The effective thermal conductivity is affected significantly by the packing density while the particle size distribution has a minor effect. The effect of thermal induced particle size change is not predominant when the temperature involved is low and within a narrow range. It is noted, however, that the current method can only consider the solid component of the heat storage system without considering heat transfer through fluid phases or by convection. In addition, PCM capsules are treated as spheres with a homogeneous property, thus a possible detailed structure of real PCM capsules is not the focus of the current work. References [1] Frédéric Kuznik, Damien David, Kevyn Johannes, Jean-Jacques Roux, A review on phase change materials integrated in building walls, Renew. Sust. Energ. Rev. 15 (1) (2011) 379–391. [2] M.K. Anuar Sharif, A.A. Al-Abidi, Sohif Mat, Kamaruzzaman Sopian, Mohd Hafidz Ruslan, M.Y. Sulaiman, M.A.M. Rosli, Review of the application of phase change material for heating and domestic hot water systems, Renew. Sust. Energ. Rev. 42 (2015) 557–568. [3] Xiwen Cheng, Xiaoqiang Zhai, Ruzhu Wang, Thermal performance analysis of a packed bed cold storage unit using composite pcm capsules for high temperature solar cooling application, Appl. Therm. Eng. 100 (2016) 247–255. [4] Hao Peng, Rui Li, Xiang Ling, Huihua Dong, Modeling on heat storage performance of compressed air in a packed bed system, Appl. Energy 160 (2015) 1–9. [5] María Asunción Izquierdo Barrientos, Heat Transfer and Thermal Storage in Fixed and Fluidized Beds of Phase Change Materials, 2014. [6] T. Saitoh, K. Hirose, High-performance phase-change thermal energy storage using spherical capsules, Chem. Eng. Commun. 41 (1–6) (1986) 39–58. [7] Keumnam Cho, S.H. Choi, Thermal characteristics of paraffin in a spherical capsule during freezing and melting processes, Int. J. Heat Mass Transf. 43 (17) (2000) 3183–3196. [8] Eduard Oró, Albert Castell, Justin Chiu, Viktoria Martin, Luisa F. Cabeza, Stratification analysis in packed bed thermal energy storage systems, Appl. Energy 109 (2013) 476–487. [9] D. Vortmeyer, R.J. Schaefer, Equivalence of one-and two-phase models for heat transfer processes in packed beds: one dimensional theory, Chem. Eng. Sci. 29 (2) (1974) 485–491. [10] To EW Schumann, Heat transfer: a liquid flowing through a porous prism, J. Franklin Ins. 208 (3) (1929) 405–416. [11] Alvaro de Gracia, Luisa F. Cabeza, Numerical simulation of a pcm packed bed system: a review, Renew. Sust. Energ. Rev. 69 (2017) 1055–1063. [12] Ciril Arkar, Sašo Medved, Influence of accuracy of thermal property data of a phase change material on the result of a numerical model of a packed bed latent heat storage with spheres, Thermochim. Acta 438 (1–2) (2005) 192–201. [13] K.A.R. Ismail, J.R. Henriquez, Numerical and experimental study of spherical capsules packed bed latent heat storage system, Appl. Therm. Eng. 22 (15) (2002) 1705–1716. [14] Yvan Dutil, Daniel R. Rousse, Nizar Ben Salah, Stéphane Lassue, Laurent Zalewski, A review on phase-change materials: mathematical modeling and simulations, Renew. Sust. Energ. Rev. 15 (1) (2011) 112–130.

515

[15] Peter A. Cundall, Otto D.L. Strack, A discrete numerical model for granular assemblies, Geotechnique 29 (1) (1979) 47–65. [16] M.L. Hunt, Discrete element simulations for granular material flows: effective thermal conductivity and self-diffusivity, Int. J. Heat Mass Transf. 40 (13) (1997) 3059–3068. [17] Watson L. Vargas, J.J. McCarthy, Heat conduction in granular materials, AICHE J. 47 (5) (2001) 1052–1059. [18] Watson L Vargas, J.J. McCarthy, Stress effects on the conductivity of particulate beds, Chem. Eng. Sci. 57 (15) (2002) 3119–3131. [19] Bodhisattwa Chaudhuri, Fernando J. Muzzio, M. Silvina Tomassone, Modeling of heat transfer in granular flow in rotating vessels, Chem. Eng. Sci. 61 (19) (2006) 6348–6360. [20] Y.T. Feng, K. Han, C.F. Li, D.R.J. Owen, Discrete thermal element modelling of heat conduction in particle systems: basic formulations, J. Comput. Phys. 227 (10) (2008) 5072–5089. [21] Y.T. Feng, K. Han, D.R.J. Owen, Discrete thermal element modelling of heat conduction in particle systems: pipe-network model and transient analysis, Powder Technol. 193 (3) (2009) 248–256. [22] John C. Steuben, Athanasios P. Iliopoulos, John G. Michopoulos, Discrete element modeling of particle-based additive manufacturing processes, Comput. Methods Appl. Mech. Eng. 305 (2016) 537–561. [23] S. Haeri, Y. Wang, O. Ghita, J. Sun, Discrete element simulation and experimental study of powder spreading process in additive manufacturing, Powder Technol. 306 (2017) 45–54. [24] H.W. Zhang, Q. Zhou, H.L. Xing, H. Muhlhaus, A dem study on the effective thermal conductivity of granular assemblies, Powder Technol. 205 (1–3) (2011) 172–183. [25] Yuanbo Liang, Haiying Niu, Yafei Lou, Expression for etc of the solid phase of randomly packed granular materials, Appl. Therm. Eng. 109 (2016) 44–52. [26] Akhil Reddy Peeketi, Marigrazia Moscardini, Akhil Vijayan, Yixiang Gan, Marc Kamlah, Ratna Kumar Annabattula, Effective thermal conductivity of a compacted pebble bed in a stagnant gaseous environment: an analytical approach together with dem, Fusion Eng. Design 130 (2018) 80–88. [27] J. Stefan, Uber einige probleme der theorie der warmeletung, Sitzer. Wien. Akad. Math. Naturw. 98 (1889) 473–484. [28] Yanio E. Milian, Andrea Gutierrez, Mario Grageda, Svetlana Ushak, A review on encapsulation techniques for inorganic phase change materials and the influence on their thermophysical properties, Renew. Sust. Energ. Rev. 73 (2017) 983–999. [29] A. Jamekhorshid, S.M. Sadrameli, Mohammed Farid, A review of microencapsulation methods of phase change materials (pcms) as a thermal energy storage (tes) medium, Renew. Sust. Energ. Rev. 31 (2014) 531–542. [30] Mohamed Rady, Granular phase change materials for thermal energy storage: experiments and numerical simulations, Appl. Therm. Eng. 29 (14–15) (2009) 3149–3159. [31] L. Xia, P. Zhang, R.Z. Wang, Numerical heat transfer analysis of the packed bed latent heat storage system based on an effective packed bed model, Energy 35 (5) (2010) 2022–2032. [32] K.A.R. Ismail, R. Stuginsky Jr, A parametric study on possible fixed bed models for pcm and sensible heat storage, Appl. Therm. Eng. 19 (7) (1999) 757–788. [33] Selvan Bellan, Tanvir E. Alam, José González-Aguilar, Manuel Romero, Muhammad M. Rahman, D. Yogi Goswami, Elias K. Stefanakos, Numerical and experimental studies on heat transfer characteristics of thermal energy storage system packed with molten salt pcm capsules, Appl. Therm. Eng. 90 (2015) 970–979. [34] J.G.H. Borkink, K.R. Westerterp, Influence of tube and particle diameter on heat transport in packed beds, AICHE J. 38 (5) (1992) 703–715. [35] George Keith Batchelor, R.W. O'brien, Thermal or electrical conduction through a granular material, Proc. R. Soc. Lond. A 355 (1682) (1977) 313–333. [36] Watson L. Vargas, Joseph J. McCarthy, Thermal expansion effects and heat conduction in granular materials, Phys. Rev. E 76 (4) (2007), 041301. [37] Xiao-Liang Wang, Dong-Yun Bai, Thermal expansion and thermal fluctuation effects in a binary granular mixture, Int. J. Heat Mass Transf. 116 (2018) 84–92. [38] Xiao-Liang Wang, Dong-Yun Bai, Thermal cycling leads grains to more homogeneous force networks and energy repartition, Powder Technol. 339 (2018) 111–118. [39] M. Michael Yovanovich, Four decades of research on thermal contact, gap, and joint resistance inmicroelectronics, IEEE Trans. Comp.Packag. Technol. 28 (2) (2005) 182–206. [40] M.G. Cooper, B.B. Mikic, M.M. Yovanovich, Thermal contact conductance, Int. J. Heat Mass Transf. 12 (3) (1969) 279–300. [41] B.B. Mikić, Thermal contact conductance; theoretical considerations, Int. J. Heat Mass Transf. 17 (2) (1974) 205–214. [42] Yuanbo Liang, Xikui Li, A new model for heat transfer through the contact network of randomly packed granular material, Appl. Therm. Eng. 73 (1) (2014) 984–992. [43] PFC3D Itasca, Particle Flow Code in 3 Dimensions, User's Guide, 2008. [44] Shazim Ali Memon, Phase change materials integrated in building walls: a state of the art review, Renew. Sust. Energ. Rev. 31 (2014) 870–906. [45] S. Tavman, I.H. Tavman, Measurement of effective thermal conductivity of wheat as a function of moisture content, Int. Comm. Heat Mass Trans. 25 (5) (1998) 733–741. [46] Zhou Zhao, K.M. Feng, Y.J. Feng, Theoretical calculation and analysis modeling for the effective thermal conductivity of li4sio4 pebble bed, Fusion Eng. Design 85 (10 −12) (2010) 1975–1980. [47] Katalin Bagi, Stress and strain in granular assemblies, Mech. Mater. 22 (3) (1996) 165–177. [48] W. Ehlers, E. Ramm, S. Diebels, G.A. dAddetta, From particle ensembles to cosserat continua: homogenization of contact forces towards stresses and couple stresses, Int. J. Solids Struct. 40 (24) (2003) 6681–6702. [49] Ching S. Chang, Matthew R. Kuhn, On virtual work and stress in granular media, Int. J. Solids Struct. 42 (13) (2005) 3773–3793. [50] X.H. Wen, L.J. Durlofsky, M.G. Edwards, Use of border regions for improved permeability upscaling, Math. Geol. 35 (5) (2003) 521–547.

516

T. Zhao, Y.T. Feng / Powder Technology 354 (2019) 505–516

[51] Marek Błasik, Numerical scheme for the one-phase 1d stefan problem using curvilinear coordinates, Sci. Res. Ins. Math. Comp. Sci. 11 (3) (2012) 9–14. [52] H.S. Carslaw, J.C. Jaeger, Conduction of Heat in Solids, Oxford Science Publications, Oxford, England, 1959. [53] A. Felix Regin, S.C. Solanki, J.S. Saini, An analysis of a packed bed latent heat thermal energy storage system using pcm capsules: numerical investigation, Renew. Energy 34 (7) (2009) 1765–1773.

[54] Y.T. Feng, Kuanjin Han, D.R.J. Owen, J. Loughran, On upscaling of discrete element models: similarity principles, Eng. Comput. 26 (6) (2009) 599–609. [55] Y.T. Feng, D.R.J. Owen, Discrete element modelling of large scale particle systemsi: exact scaling laws, Comp. Particle Mech. 1 (2) (2014) 159–168. [56] R.G. Craig, J.D. Eick, F.A. Peyton, Properties of natural waxes used in dentistry, J. Dent. Res. 44 (6) (1965) 1308–1316.