Numerical modelling of melting and settling of an encapsulated PCM using variable viscosity

Thermal energy storage units using macro-encapsulated PCM in industrial and residential applications are contemporary due to better efficiency during charging and discharging. This article focuses on numerical modelling of the melting process in a macro-encapsulated PCM. Accounting the non-linear enthalpy–temperature relation and ramping down the velocity in solid phase is therefore fundamental. In the present article the variable viscosity method is implemented to ramp down the solid velocity and allow settling of the solid phase. This complete numerical model of melting and settling of PCM in a capsule is implemented in OpenFOAM. The numerical results for different solid viscosities are validated with experiments and show good agreement. The influence of the solid viscosity value and the pressure–velocity convergence is studied. It is observed that the pressure–velocity convergence only plays a greater role in the case where the computation of the exact solid velocity is required.


Introduction
Statistics reveal that the global energy demand is set to grow by 37% until 2040 [1]. This increasing requirement and the decrease in production of fossil fuels challenge industries and scientists to produce clean renewable energy to meet this deficit. Reusing waste heat released as a biproduct during an industrial process may help to meet the phase change phenomena contain three phases namely solid, liquid and mushy i.e. the interface where both liquid and solid coexist. Therefore taking the non-linear enthalpy and temperature relation into account is a challenging task. A few iterative approaches have been employed to solve the temperature-enthalpy coupling [10,12,13]. Apart from modelling the non-linear enthalpy and temperature relationship ramping down velocity in the solid phase is also required. Methods used to ramp down solid velocities are generally called velocity switch-off techniques.
The velocity overwrite approach proposed by Morgan [9] overwrites the computed velocity in the cell with solid PCM [14]. Abrupt change in the phase change at a very small temperature range can produce instabilities in this approach. In the porosity approach, the discharge flux from Darcy law is implemented in the momentum equation as a source term [12]. The source term takes a very large value in cells with full solid phase and zero in cells with full liquid phase and thereby ramps down the velocity in the solid phase. Gartling [11] proposed the variable viscosity method (VVM) on considering viscosity, µ as a function of temperature. Viscosity is generally associated with fluids but in this method the solid phase is proposed to have higher viscosity values. The finite volumes with only solid phase show a very high resistance to the inertial forces due to the consideration of a large viscosity. The suitable equation proposed by Voller et al. [14] is where B is a large constant, µ l is the viscosity of the liquid phase, L is the latent heat value and H l is the enthalpy of the fluid phase in a finite volume. In contrast to the most applications the solid PCM is not fixed within macro-capsules and settles down during melting, resulting in close contact melting. To take account of the settling Assis et al. [6] implemented a low Darcy term value. Therefore, the velocity in the solid PCM is not suppressed completely resulting in a slow settling solid. Due to the low Darcy constant the melting front of the PCM in the capsule is affected, as it influences the melting rate. Later Rösler [7] implemented a Darcy term value with a settling velocity as an additional source term. Whereby, the settling velocity is calculated by an equilibrium force balance of the gravitational force and the pressure on the solid surface earlier proposed by Asako et al. [5] to reduce the computational effort. Later Ghasemi and Molki [15] have studied the influence of the dimensionless numbers influencing the solid settling. But, instabilities can occur due to the forced settling model during a dynamic melting rate influenced by solid settling. Therefore, this article does not only discuss the implementation of VVM with solid settling but also its convergence and influence on the melting rate.

Numerical model
The governing Navier-Stokes equations to solve incompressible flow considering the Bousinessq approximation are [16] Conservation of mass: Conservation of momentum:

Conservation of energy:
The source term S B which is equal to ρgβ(T − T ref ) represents the change in density ρ causing natural convection is implemented through the Bousinessq approximation [17]. Whereas u is the velocity, p is the pressure, t is the time, S i is an additional source term used to represent any added forces in a physical phenomenon and β is the expansion coefficient of liquid relevant to a reference temperature T ref . In the energy equation, T is the temperature, κ the thermal conductivity and H the total enthalpy.

Enthalpy-temperature
The non-linear relation between the total enthalpy and temperature is approximated by a piecewise linear function as depicted in Fig. 1. The total enthalpy in Eq. 4 can be segregated as sensible and latent heat, Voller and Prakash [10] proposed rearranging them as where c is the specific heat capacity of the PCM, which is equal to c s · α s + c l · α l . Here, c s and c l are the specific heat capacities of solid and liquid PCM respectively. Within the source based fictitious method discussed in this article, see Voller and Prakash [10] the latent heat is expressed as Here α l is the volume fraction of liquid PCM and the solid volume fraction α s is equal to 1 − α l . The lower and upper limits of the melting temperature are T s and T l respectively. The rearranged energy equation reads [10] (2) ρ ∂cT ∂t + ∇ · (ucT ) = ∇ · (κ∇T ) − ∂ρH l ∂t + ∇ · (ρuH l ) .
In the source based method, the term ∂ρH l ∂t + ∇ · (ρuH l ) is the latent heat source term computed and corrected using an iterative approach [13]. The discretized form of Eq. 7 considering Eq. 5 is In this discretized equation, P is the node of the control volume calculated, NB represents the neighbouring cells of the node P, ψ NB is equal to a NB T n+1 NB which is the product of a NB , the coefficients in the discretized energy equation at neighbouring nodes to node P and T n+1 NB , the temperature neighbour to the node P, a P is the coefficients in the discretized energy equation at node P, n + 1 is the present time step, n is the old time step and V is the volume of the computing cell. After every iteration the liquid volume fraction is updated through the calculated temperature [13,18]. The liquid fraction of the fluid is updated by In the Eq. 9, URF is the under relaxation factor to relax the iteration steps. It should be noted here that the choice of the URF influences the stability and accuracy of the solution and its value lie between 0 and 1. Equations 8 and 9 are iterated in a loop till convergence [18]. The correction bounds the volume fraction of the liquid between 0 and 1 [13]. In this way the total enthalpy is conserved between its latent and sensible parts maintaining the non-linearity with temperature.

Modelling solid settling
Modelling solid-liquid phase change also demands ramping down solid velocity i.e., the velocity in the solid phase.
In general, viscosity is defined as the ability of a fluid to resist the inertial forces due to gradual deformation by shear and tensile stresses [19]. In the VVM, the ability of the solid to resist the inertial forces is defined as a very large value i.e., a very high viscosity. The cell volumes containing a solid phase with very high viscosity values will resist all possible inertial forces in the cell volume due to which the solid movement is obstructed. This method is a stable approach due to its mean value implicitly computed from the cell faces [20]. The discretized form of the viscous term in the transport equation is as below [20] and pictorially represented in Fig. 2: where φ is a conserved quantity, S is the surface area of the finite volume cell, f are cell faces i.e., the surfaces enclosing the control volume and also the surfaces in contact with the neighbouring cells and Γ the diffusive scalar quantity. The relevant discretization scheme can be used in order to discretize the viscous term in the equation. A more appropriate form to compute the viscosity would be the harmonic mean of the values in cells [14]. However, the stability is also due to the interpolation schemes of the viscous term and cannot be generalized. The stability of the process is also influenced by the mean viscosity which is computed from each neighbouring cell and face in contact with the P as shown in Eq. 10 where Γ is µ in the momentum equation: where P is the computing cell and NB and the computing cell P.
The relation between temperature and viscosity can be represented using various functions. The fluxes on the phase change interface i.e., the mushy zone, influence the nature of the melting. The piecewise linear approach considered in this article defines the total viscosity as fractions of solid and liquid viscosity. The viscosity µ in the momentum Eq. 3 reads  Finite volume discretization and approximation of viscous term in a finite volume cell [21] In Eq. 12, viscosity is considered as a very large value in the full solid phase and experience a linear declination in the fluid phase as shown in Fig. 3.
In the VVM the shear velocity in the flow is restricted by a large viscosity value in the mushy zone and solid phase. The shear velocity is ramped down due to the viscous term of the momentum equation, The source term F i in Eq. 13 is a source term implemented to enhance the natural convection and settling of solid due to the density difference: with and where ρ s and ρ l are the densities of the solid and liquid PCM respectively. Total density in a cell ρ i in Eq. 14 is variable whereas ρ in Eq. 13 is constant. The source term F i is implemented in the pressure-velocity loop implicitly. Settling of solid in a melting process can be accounted by correcting and iterating the pressure-velocity loop [5,22]. Asako et al. [5] defined the settling force of the solid in a melting phenomena with an equilibrium equation which is the sum of inertia, gravity, pressure and shear forces. In this article implementation of the settling force from equilibrium equation see [5], is refrained. But the settling of the solid is enhanced by iterating the pressure-velocity loop in the momentum equation.
The complete numerical model is implemented in an open source C++ finite volume library OpenFOAM-2.2.2. A combination of semi-implicit method for pressure linked equations (SIMPLE) [16] and pressure implicit with splitting of operator (PISO) [23] namely PIMPLE algorithm is used to solve the segregated pressure and velocity equations.

Numerical setup and problem definition
A relevant computational grid of the geometry is necessary to simulate the developed model. Different computational grids have been tested to achieve a grid independent solution. The cubical cavity 40 mm × 40 mm × 40 mm with a geometrical grid consisting of 10,000 cell volumes assigned 100 × 100 × 1 cells is selected for the present work. In the simulations the front and back boundaries of the cubical cavity are considered symmetric. Due to the symmetric front and back boundaries the geometry is a virtual 2D geometry contributing 1 cell in z-direction. The PCM RT35HC from Rubitherm [24] is considered in both, numerical simulations and experiments. The measurement of the material properties as seen in Table 1, as well the experiments were performed in a previous work, for more details see Rösler [7].
Two different thermo-fluidic cases have been simulated using the numerical model and the results are compared with the experimental results. The melting rate of the PCM depends on the heat flow into the material. Therefore a more exact mixed heat transfer boundary condition or Dirichlet-Neumann boundary condition is implemented i.e., the heat transfer around the boundary using where n is the normal direction to the heat transfer, HTC is the heat transfer coefficient equalling to 2500 W/(m 2 K) represents the heat transfer through the copper plates in the experiments. T i the boundary element temperature, T ∞ the ambient temperature at the boundary and κ the thermal conductivity of the PCM.

Heated left boundary
The idea of validating the melt front by heating a wall is empirical in the field of moving boundary problems. The melt front and melting rate obtained through numerical simulations are validated with experimental results. In the heated left boundary case the left wall boundary of the geometry is heated with a temperature of T h = 315.15 K whereas the right wall is heated with temperature of T c = 303.15 K. The top and bottom boundaries have been defined with an adiabatic temperature T a in Table 2. On these boundaries, a zero heat flux, Neumann boundary condition is implemented where the temperature gradient on the boundaries is zero (Fig. 4). The solid viscosity is concluded in this case by validating the different melt fronts obtained through different solid viscosities with experimental results.

Close contact melting
Close contact melting refers to the detachment of the solid from the boundary surface and the movement of solid into the liquid phase resulting in a small melting column [25]. Close contact melting enhances the melting rate of the PCM in the geometry. The considered case in this article represents the realistic melting and settling of PCM in a capsule when all boundaries except for the front and the back are heated with a constant temperature (Fig. 5).
All boundaries are heated with an equal temperature of T h = 315.15 K, as shown above in Table 3. In this section, the melting front and also the settling form of the solid PCM in a capsule is analysed and validated with experimental results.

Experiments
The melt front in the capsule is captured using a single lens reflex camera. The solid-liquid interface has been marked and the melted volume has been calculated with time. The direction of lighting in the capsule is parallel to the gravity in the case discussed in Sect. 3.1 and perpendicular to the gravity in the case discussed in Sect. 3.2. In the experiments performed for the Sect. 3.1, a rubber foam insulation with a low thermal conductivity of 0.04 W/m K has been placed on the top and bottom boundaries to support the no heat flux condition in the simulations [7]. This insulation with a thickness of 40 mm has been placed on the top and bottom boundaries resulting in a neglectable low heat transfer. Here, the liquid fraction of PCM portion at different times is calculated by marking the phase front between liquid and solid PCM on the captured pictures obtained from the experimental set-up (Fig. 6).

Results and discussion
The computational grid is generated to simulate the melting phenomena with the developed numerical model. It has been observed that the melt front and melting rate of the PCM are largely influenced by the solid viscosity value and are compared with the experimental results.

Heated left boundary
Solid viscosities ranging from 10 3 to 10 6 Pa s are considered for the set-up described in Sect. 3.1. It is observed that the interface curvature between the solid and liquid phase depends on the solid viscosity value which is also relatively proportional to the melting rate of the PCM. The solid-liquid interfaces obtained through different viscosities are compared with the melt interface from experiments as shown in Fig. 7.
The results for a solid viscosity of µ s = 10 5 Pa s have shown satisfactory results. The melting rate of the PCM using this solid viscosity is quantitatively validated against experimental results in Fig. 8. The numerical results show a good agreement with the experiments.

Close contact melting
The close contact melting case introduced in Sect. 3.2 is simulated considering the defined boundary conditions. Two cases, one with a large viscosity value obtained from the validated melt front in Fig. 7 and another with a lower viscosity to enhance the melt rate have been simulated in this section.

Larger solid viscosity
The solid viscosity µ s = 10 5 Pa s, obtained by validating the melt fronts from Sect. 4.1 is implemented in this section to simulate the melting of PCM and settling of solid (Figs. 9, 10).
The qualitative validation of the results has shown some decent agreement with the experiments. A small deviation from the experiments could be due to the lack of the exact viscosity selection in the mushy zone. Modelling the right viscosity in the mushy zone using a linear correlation of total viscosity is difficult. Therefore this article refrains from discussing the influence of mushy viscosity on the  melting rate. The movement of the solid is restricted by obstructing the shear velocity but the settling velocity is not restricted as shown in the Fig. 11.
The settling velocity of solid during melting is variable and changes with the melting rate of the PCM. Due to the implementation of a large solid viscosity value tracking the solid-liquid interface during settling demands a large number of iterations. The iteration count depends upon the velocity differences between the solid and liquid phases. Pressure-velocity loop iterations between 25 and 100 are found to be optimal to introduce the right settling of solid PCM as shown in Fig. 12. The change in results after 25 pressure-velocity iterations is found to be very small and can be neglected here.
The melting rate of the simulations show some good agreement with experiments after t = 10 min as shown in Fig. 13. In the first 10 min the different settling form in the experiments is due to the fixed solid at the beginning which starts to settle after about 4 min. Therefore, the simulation results diverge from experiments in this time duration. An asymmetrical melting and settling of solid in the experiments in Figs. 9 and 10 is due to the liquid drag force affecting the free degree of motion of solid PCM. These liquid drag effects on solid PCM are neglected in the simulations. A much larger value of solid viscosity of PCM may resist non-physical deformation during settling of solid PCM. However, the choice of a larger solid viscosity demands a large computational power due to iterations required to converge the pressure-velocity in tracking the interface between solid and liquid PCM. Therefore, a small viscosity value of µ s = 5 × 10 3 Pa s has been chosen in the next section to influence the convergence rate. Low solid viscosity also proportionally decreases the value of viscosity in the mushy zone.

Reduced solid viscosity
The implementation of a low solid viscosity value influences the melting rate and melting form due to the low friction offered by the solid. The choice of a low solid viscosity is not in agreement with the validated solid viscosity obtained in Fig. 7 but due to the implementation of reduced solid viscosity, only a very few pressure-velocity iterations are required for a converged solution (Fig. 14).
The results from the numerical model using a low viscosity value show a good acceptance with the qualitative experimental results until t = 6 min. It has been observed from the simulations that the solid PCM settles at the bottom of the capsule geometry at this time. Afterwards the solid PCM experiences a deformation due to its own weight and implementation of the low solid viscosity value. In this article, the effective viscosity is defined as piecewise linear as seen in Fig. 3. In this case, the reduction of the viscosity value in the solid phase also proportionally decreases the viscosity value in the mushy zone. Due to which that the material in mushy zone close to the melting temperature tends to flow resulting in the deformation observed in Fig. 15.
Quantitative validation of the results using the low solid viscosity produce satisfactory results as seen in Fig. 16 apart from the deformation observed in the solid in Fig. 15. This material flow deformation of the solid PCM occurred post contact of the solid at the bottom of capsule due to its own weight.

Conclusions
This article focuses on the implementation of a variable viscosity to model solid-liquid phase change problems  including the settling of the solid phase. Different test cases were simulated to validate the numerical model and also to investigate the influence of the solid viscosity value on the melting process. As the choice of these solid viscosities was random and a generalization could not be made, the developed numerical model has shown a decent agreement with the experimental results. Validating the numerical model by comparing melt fronts and the melting rate is empirical, see Voller and Prakash [10]. The choice of solid viscosity in the present article is also concluded in the same way. The influence of the solid viscosity value on the melting phenomena is found to be different based on the application. A higher solid viscosity value has shown a good agreement with experiments in Sects. 4.1 and 4.2, whereas lower solid viscosity has only shown a good agreement with experiments in Sect. 4.2. Within the developed model compromise has to be made between the computational time and the accuracy of the results as the large viscosity value demands a large computational time to model settling. Whereas, the low viscosity value requires a low computational time, but the solid PCM tends to deform due to small viscosity influencing the outflow of solid as shown in Fig. 15.
Concerning this numerical model, it can be further extended to different geometrical applications, viscosity relations and force implementations. It can be strongly concluded that a variable viscosity method (VVM) is a stable approach but the error made due to the linear relationship of the viscosities could be optimized by choosing a relevant function depending on the problem definition. So, there is a strong demand for a detailed investigation on the nature of these ramp down techniques.