Numerical Solution of Heat Transfer Process in PCM Storage Using Tau Method

Thermal energy storage units that utilize phase change materials have been widely employed to balance temporary temperature alternations and store energy in many engineering systems. In the present paper, an operational approach is proposed to the Tau method with standard polynomial bases to simulate the phase change problems in latent heat thermal storage systems, that is, the two-dimensional solidification process in rectangular finned storage with a constant end-wall temperature. In order to illustrate the efficiency and accuracy of the present method, the solid-liquid interface location and the temperature distribution of the fin for three test cases with different geometries are obtained and compared to simplified analytical results in the published literature. The results indicate that using a two-dimensional numerical approach can predict the solid-liquid interface location more accurately than the simplified analytical model in all cases, especially at the corners.


Introduction
Storage of thermal energy by latent heat during the melting or solidification process of phase change materials (PCMs) occurs in several practical application areas, from solar water heating and air conditioning in buildings to electronic cooling and the automobile industry.In applications with a small temperature swing, latent heat thermal storage (LHTS) is preferred to sensible heat storage because of the highenergy storage capacity of PCMs and their nearly isothermal storage mechanism.Paraffins, salt hydrates, and fatty acids are the commonly used PCMs in thermal energy storage (TES) systems.Storing or releasing large amounts of energy during phase change is the greatest advantage of PCMs [1].
In the late 19th century, J. Stefan derived the formulation to determine the temperature distribution and moving solidliquid interface history of a solidifying slab of water.Since then, particularly in the past three decades, the problem bearing his name has been extended to include such complex phenomena as TES, the solidification of alloy systems, and supercooling [2].Generally, heat transfer including phase change is referred to as a moving boundary problem (MBP), which is a transient, nonlinear phenomenon.The nonlinearity is the main challenge in MBPs and therefore analytical solutions for these problems have only been generated for specific situations that have a simple geometry and simple boundary conditions [3][4][5].To the knowledge of the authors of this paper, there is no exact analytical solution for the solidliquid interface location in two-dimensional PCM stores with internal fins.Some approximate analytical solutions for this problem published in the literature have simplified the twodimensional heat transfer problem to a one-dimensional one, thereby driving an analytical solution in two separate region [1][2][3][4].The most well-known precise analytical solutions and analytical approximations for one-dimensional MBPs with various boundary conditions have been published by Özis ¸ik [6] and Alexiades and Solomon [7].
PCMs have some disadvantages, such as their low thermal conductivity.In addition, conduction is the main heat transfer mechanism evident during the solidification process.Due to the temperature difference in a liquid PCM, there is natural convection at the solid-liquid interface, but this has a negligible effect on the solid-liquid interface location in comparison with heat conduction in a solid PCM [8,9].The most common way to enhance the internal heat transfer of PCM storage is through the use of fins [10].An exact 2 Mathematical Problems in Engineering analytical solution for the temperature distribution and solidliquid interface location in two-dimensional PCM storage with internal fins has not yet been developed.Lamberg and Sirén [11] simplified the two-dimensional phase change problem to a one-dimensional one and derived an analytical solution that predicts the solid-liquid interface location and temperature distribution of the fin for the melting process.This includes a constant imposed end-wall temperature in a semi-infinite PCM storage with internal fins.The analytical and numerical results showed good agreement for three storage sizes.
In addition to simplified analytical solutions, mathematical models based on the effective heat capacity [2] or the enthalpy formulation have been developed to address the two-dimensional phase change problem, along with the finite difference or finite element procedure.Talati et al. [3] presented a two-dimensional numerical solution based on the enthalpy method for the solidification process in rectangular finned storage with a constant heat flux boundary condition.The numerical results agreed satisfactorily with the approximated analytical ones.Mosaffa et al. [1] solved the same problem with convective cooling boundaries numerically based on the enthalpy method.Approximate analytical and two-dimensional numerical predictions of the solid-liquid interface location for the different test cases were compared to evaluate the accuracy and performance of the proposed solution.Chiba [5] applied a series solution for the heat conduction problem with phase change in a finite slab.To verify the numerical solution, the results were compared with analytical solutions published in the literature.As discussed above and in related papers [12,13], it can be said that the phase change problem has been studied numerically, experimentally, and semianalytically by many authors.
In addition to the mentioned numerical methods, linear and nonlinear heat conduction problems in TES systems can be solved by the operational Tau method, which can be described as a modification of the spectral Galerkin method with arbitrary polynomial bases [14][15][16][17].The main advantage of this technique is using simple operational matrices that reduce the computational costs remarkably, since only nonzero elements of operational matrices are required [17,18].The Tau method was developed by Lanczos [19].Moreover, Ortiz and Samara [20] proposed an approximation technique for solving nonlinear ordinary differential equations using the operational Tau method.In recent years, the operational Tau method has been extended for the numerical solution of different kinds of ordinary differential equations (ODEs) [21,22] and partial differential equations (PDEs) [23,24].In addition, various authors have implemented the Tau method for solving integral and integrodifferential equations [25,26].Hosseini et al. [18] developed the operational Tau method for solving one-dimensional nonlinear transient heat conduction equations.In their paper, some numerical examples were given to clarify the efficiency and accuracy of the proposed method.Talati et al. [27] proposed the numerical solution of one-dimensional transient heat conduction equation with variable thermophysical properties using the Tau method.Finally, some numerical examples have been solved by the proposed method and the results were compared with solutions obtained by the other methods.Saadatmandi and Razzaghi [28] proposed an approximate method for solving the diffusion equation with nonlocal boundary conditions.The method was based upon constructing the double shifted Legendre series to approximate the required solution using Legendre tau method.Finally, numerical examples were included to demonstrate the validity and applicability of the method and a comparison was made with existing analytical results.
In this work, the two-dimensional solidification process is investigated in a rectangular PCM storage with horizontal internal fins (see Figure 1).
The purpose of the present study is to develop a numerical solution for heat transfer problem in TES systems based on the operational Tau method.The remainder of this paper is organized as follows: Section 2 describes the problem configuration, computational domain, governing equations, and assumptions.According to the evaluation of the accuracy of this method, a simplified analytical solution is presented for this problem in Section 3. In Section 4, some preliminaries of the Tau method are introduced.Then, the unknown functions and their derivatives are extended as truncated series to obtain a matrix form of the given differential equations and boundary conditions.Finally, some conclusions are drawn in Section 5.

Problem Statement
Consider PCM storage with internal straight metal fins filled with a phase change material.Initially, the storage exhibits a uniform temperature   .Because of the symmetrical structure of storage and the imposed boundary conditions, the analysis is reduced to a single symmetrical cell, as shown in Figure 2, where the planes of symmetry are in the middle of the fin and midway between the two adjacent fins.The heat equation for the solid PCM during solidification process with boundary and initial conditions can be written as follows [2,11]: (, , ) =   (, , ) .
Due to the symmetry, we have Similarly, the heat equation for a fin with boundary and initial conditions can be written as (, , ) =   (, , ) .
Due to the symmetry, we have The energy balance for the solid-liquid interface with the initial conditions is as follows [6,29]: where Here,  is the temperature,  is time,  is specific heat,  is density,  is heat conductivity,  is the latent heat of fusion,  is the location of the solid-liquid interface, and   is initial temperature.The subscripts , , and  denote the solid PCM, fin, and wall, respectively.Due to the nonlinear and unsteady nature of the problem, several assumptions have been made to simplify the twodimensional heat-transfer problem.The assumptions are the following: (1) The solidification process is assumed to be isothermal.
(2) Initially the storage is filled with liquid PCM at the solidification temperature (  =   ), and the fin also is initially at   .Hence, the heat conduction and free convection in the liquid phase are assumed to be negligible [3].
(3) The thermophysical properties of the PCM and fin are assumed to be constant.
The problem is a classical Stefan-type problem in which the nonlinearity is the source of difficulties; it destroys the validity of generally used solution methods such as the superposition and separation of variables approaches [7].

Approximate Analytical Approach
Lamberg divided storage into two separate regions.The first region is far from the fin and only the heat sink has been assumed as the constant temperature end wall, so heat transfer only occurs in the -direction.The other region is near to the fin and both the wall and the fin transfer the heat from the PCM to the environment; hence, in this region the solidliquid interface moves one-dimensionally in the -direction.
In the first region, the distance of the solid-liquid interface from the end wall is as follows [2]: Meanwhile, in the second region, the distance of the solidliquid interface from the fin is as follows [2]: where
Since the proofs of all lemmas are similar, we only prove (26) using the definition in (19).Thus, we have the following: where   is the ( + 1)th column of the matrix .
To prove (32), we first find the matrix form of (/) using ( 24) and (31).Thus, we have Equation ( 32) results from substituting / for  in (33) and by the definition of =   .By considering some definitions and lemmas in ( 22)-(33), we recall the Tau method to obtain the matrix representation of the governing equations and boundary conditions, (1)- (12).The matrix form of (1) can be produced by using ( 16), ( 22), (26), and ( 27): With the same procedure, by using ( 16), ( 23), (26), and ( 27), the matrix form of ( 6) is expressed as follows: Finally, the matrix form of the energy balance for the solidliquid interface results from substituting ( 17), ( 22), ( 24), ( 29), and ( 32) into ( 11): The matrix forms of the boundary conditions are obtained by the same procedure.The unknown coefficient matrices,   ,   , and   , at  = 0 are obtained by the initial conditions.For  > 0, the unknown coefficient matrices are calculated by solving a system of algebraic equations at each time step.Thus, the solid-liquid interface location and the temperature distributions of the solid PCM and fin are determined through the time.We can find the unknown coefficient matrices by solving a nonlinear system of equations.Note that the relations contain 2( + 1) 2 + 9( + 1) equations, whereas we need only 2( + 1) 2 + ( + 1) equations to find the unknown coefficients.Thus, we first select algebraic equations obtained from the boundary conditions; the remaining equations have been selected from algebraic equations obtained from the differential equations.

Results and Discussion
The PCM considered in this study is paraffin and the fins are made of aluminum.Table 1 shows the properties of the PCM and the fin material.
In order to test the accuracy and performance of the twodimensional numerical solutions by the Tau method, three test cases with different geometries have been selected.For convenience, the length of the fin divided by half of the height of the cell is introduced as the cell aspect ratio,  =   /  .The  cell aspect ratio  is variable in the different cases, taking on values of 0.2, 3, and 5.The dimensions for the storage used in the simulation are shown in Table 2.In all cases, the half thickness of the fin, , has a constant value of 0.5 mm and the initial temperatures of the PCM and fin are equal to 25 ∘ C. The wall temperature is set at 10 ∘ C. The two-dimensional numerical predictions of the solidliquid interface location and the temperature distribution of the fin obtained from the Tau method ( = 12) were compared with the simplified analytical solution.The computations associated with these experiments were performed in Maple 15.The time step was 1 s and the convergence was checked at each time step, with a convergence criterion of 10 −4 for temperature.
Figures 3-5 show the solid-liquid interface locations for the three test cases with different aspect ratios ( = 0.2, 3, and 5).As mentioned in Section 3, the analytical solution for this problem is one-dimensional in two separate regions inside the solving domain; thus, as seen in the figures, this method predicts the solid-liquid interface locations imprecisely.Meanwhile, because of the two-dimensional nature of the heat transfer process in the PCM storage, the two-dimensional numerical approach gives more accurate results compared to the analytical solutions.
The solid-liquid interface location in case 1 ( = 0.2) is obtained using the analytical and numerical methods at  = 200 s and  = 400 s.In the region far from the fin, the numerical results exhibit good agreement with the analytical ones because the dominant heat transfer is in the -direction in this region.However in the region near to the fin, the numerical solution gives more accurate results compared to the analytical solution because of the two-dimensional heat transfer effects.The solid-liquid interface location at  = 200 s and  = 400 s in case 2 ( = 3) is shown in Figure 4.The interface in the two-dimensional numerical approach moves more quickly than approximate analytical model, except in the region near to the fin.This is because, in reality heat transfer in the storage is two-dimensional and accelerates the solidification process.In this case, the effect of two-dimensional heat transfer is more obvious than in  case 1 ( = 0.2), and heat transfer is significant in both direction and -direction.Thus, the numerical solution gives a rounder shape to the interface than the analytical solution.
In Figure 5, the solid-liquid interface locations in case 3 are obtained using the analytical and numerical methods at  = 300 s and  = 900 s.By increasing  up to 5, the main heat transfer is in the -direction, and so in the region near to the fin, the numerical results exhibit good agreement with analytical ones.However, in the region far from the fin and by increasing the heat transfer effects in both directions, the numerical solution gives more accurate results compared to the analytical solution.As shown in Figures 3-5, the dimensions of the storage have a significant effect on the location and speed of the solidification interface.
When the cell aspect ratio  is much smaller than unity ( ≪ 1), the heat transfer occurs mainly in the -direction.When  is much greater than unity ( ≫ 1), the heat transfer is dominant in the -direction.The heat transfer is remarkable in both directions when the cell aspect ratio  tends to unity.Tables 3 and 4 show the solid-liquid interface distances from the wall and fin, respectively.Distances in the numerical solution are reported on the horizontal symmetry line in Table 3 and on the vertical symmetry line in Table 4, where  1 = 200 s and  2 = 400 s in cases 1 and 2 and  1 = 300 s and  2 = 900 s in case 3.
The temperature distribution of the fin is shown in Figures 6-8 for the test cases.As shown in all cases, the high thermal conductivity of the fin causes the temperature variation along the fin to be low.When the fin is short, it achieves the end-wall temperatures quickly; thus, in case 1, the temperature distribution along the fin is approximately fixed.By increasing the fin length, the temperature difference along the fin is increased (see cases 2 and 3).When the fin is long, the solidification process takes more time; thus, the chosen time steps are higher in case 3 than in cases 1 and 2. In case 1, where  is small, the heat transfer is nearly one-dimensional in the -direction and the proposed numerical model for the temperature of the fin exhibits a good agreement with analytical results.In this case, the maximum difference between the analytical and numerical solutions is approximately less than 0.04 ∘ C at  = 200 s and 0.06 ∘ C at  = 400 s.The differences between the analytical and numerical solution are small for the temperature distribution of the fin in all cases.
Tables 5-7 show the fin temperature at specific times and locations for all three cases.As shown in the tables, the largest error occurs for case 3 when  = 300 s and is about 0.5 ∘ C.

Conclusions
In this paper, a numerical solution based on the operational Tau method was developed for the solidification process of a PCM encapsulated in a rectangular container with horizontal internal fins.The main findings were as follows: (1) The dimensions of the storage have a significant effect on the location and shape of solidification interface.
(2) When the cell aspect ratio is much smaller or greater than unity, the heat transfer occurs mainly in one direction, but the heat transfer is significant in both directions when  tends to unity.Therefore, a two-dimensional approach can be used to predict the solid-liquid interface location more accurately than the simplified analytical model, especially at the corners in all cases.
(3) The differences between the analytical and numerical solutions are small for the temperature distribution of the fin in all cases.The largest error occurs for case 3 when  = 300 s, and this is about 0.5 ∘ C.
The numerical results indicate that the Tau method can be regarded as a structurally simple approach that exhibits better performance that is conventionally applicable to the numerical solution of linear and nonlinear multidimensional heat conduction equations with moving boundaries.

Figure 1 :
Figure 1: Schematic of PCM storage with internal fins.

Figure 3 :
Figure 3: Position of the solid-liquid interface for the approximate analytical model and the two-dimensional numerical method in case 1 ( = 0.2).

Figure 4 :Figure 5 :
Figure 4: Position of solid-liquid interface for the approximate analytical model and the two-dimensional numerical method in case 2 ( = 3).

Figure 6 :
Figure 6: The temperature distribution of the fin in the PCM storage in case 1 ( = 0.2).

Figure 7 :Figure 8 :
Figure 7: The temperature distribution of the fin in the PCM storage in case 2 ( = 3).

Table 2 :
Dimensions of the storage considered.

Table 3 :
Solid-liquid interface distance from the wall.

Table 4 :
Solid-liquid interface distance from the fin.