A Series Solution for Heat Conduction Problem with Phase Change in a Finite Slab

and Applied Analysis 3 Table 1: Fundamental operations of two-dimensional differential transform. Original function y(t, x) Transformed function Y(k, s) y (t, x) = α ⋅ f (t, x) Y (k, s) = α ⋅ F (k, s) y (t, x) = f (t, x) ± p (t, x) Y (k, s) = F (k, s) ± P (k, s) y (t, x) = ∂f (t, x) ∂t Y (k, s) = k + 1 H ⋅ F (k + 1, s) y (t, x) = ∂f (t, x) ∂x Y (k, s) = s + 1 G ⋅ F (k, s + 1) y (t, x) = ∂ V+w f (t, x) ∂tV∂xw Y (k, s) = [(k + V)!/k! ⋅ (s + w)!/s! ⋅ F (k + V, s + w)] HVGw y (t, x) = f (t, x) p (t, x) Y (k, s) = k


Introduction
There is a strong demand for the analyses of heat conduction problems with phase change in a broad range of fields such as ice thermal storage, refrigeration and thawing of foods, freeze dehydration, freeze-drying, cryosurgery, and freeze preservation of living materials.In addition, this type of analysis is important in terms of the quality and productivity estimates of casting products.Unlike normal heat conduction problems, phase change problems are characterised by nonlinearity due to the motion of the change-of-phase front.Thus, exact solutions can be obtained only for a few cases.In particular, when heat flows in both the liquid and solid phases are considered (i.e., the two-phase problem) and/or the object to be analysed is confined to a finite region, the analysis difficulty increases.An excellent textbook and two review articles have been published on the mathematical modelling of phase change problems [1][2][3].
Thus far, only a limited number of researchers have studied the two-phase problem using analytical methods.It should be noted that an important study on a semiinfinite slab by Neumann [1] underlies all the analytical studies.Zubair and Chaudhry [4] presented a closed-form solution to phase change problems with a convective boundary condition.Oliver and Sunderland [5] and Briozzo and Tarzia [6] studied a two-phase problem in which the material properties of the phases were linear and nonlinear functions of the temperature, respectively.Using perturbation methods, Yang et al. [7] analysed the solidification of a finite slab with shrinkage that was convectively cooled at one end.Dursunkaya and Nair investigated the solidification of a finite medium subjected to a constant boundary temperature [8] or a periodic boundary temperature [9].McCue et al. [10] analysed the two-phase problem for spheres by applying the method of matched asymptotic expansions.Das and Rajeev [11] derived an analytical solution of the two-phase problem in a finite domain by using the finite sine transform technique.An and Su [12] developed a lumped parameter model and analysed the melting of a finite slab with volumetric heat generation.Most recently, the phase change problem with location-dependent latent heat was studied [13].
Since the late 1990s, the differential transform method (DTM) has been attracting attention.It is a powerful tool based on Taylor series expansion and is used to explicitly solve not only linear differential equations but also nonlinear ones.This method yields a solution in a simple power-series form.The main advantage of this method is its direct applicability to nonlinear differential equations without requiring linearisation, discretisation, or perturbation.Although there is a criticism that the DTM is purely and solely the traditional Taylor series method [14], applications of the DTM to nonlinear heat conduction analyses are increasingly reported.
A brief review of the relevant literature published before 2011 can be found in our previous paper [15], which analysed the steady temperature field and related thermal stresses in an annular disc of variable thickness with temperaturedependent parameters.Research achievements in the past few years are described briefly as follows.S. Mosayebidorcheh and T. Mosayebidorcheh [16] analysed steady heat conduction in a straight fin with temperature-dependent thermal conductivity, which loses heat by simultaneous convection and radiation.Torabi et al. extended their work to a moving fin [17], a T-shaped fin that requires more nonlinear terms [18], and longitudinal fins of variable thickness [19].Tabatabaei et al. [20] solved heat-like and wave-like equations with variable coefficients to show the efficiency and simplicity of the DTM.Ndlovu and Moitsheki applied the one-dimensional DTM and two-dimensional DTM to steady [21] and transient [22] heat conduction analyses for fins of variable thickness with temperature-dependent parameters, respectively.Recently, the application range of the DTM has been extended to the analysis of heat conduction in nonhomogeneous bodies [23,24].
In the present paper, the two-dimensional differential transform method is applied to solve one-dimensional transient heat conduction problems with phase change.In particular, the two-phase problem for a slab of finite thickness with temperature-dependent material properties is analysed using the apparent specific heat method [25].An approximate analytical (series) solution for the temperature profile of a melting or solidifying slab is derived.Numerical results illustrate the effects of the temperature-dependent parameters on the transient temperature profile of the slab.

Two-Dimensional Differential Transform Method
The basic theory of the two-dimensional (2D) DTM is briefly described here.We consider (, ) to be an analytic function in the domain of interest.The domain of interest is a time-space domain whose upper limits are  for the time coordinate and  for the space coordinate, with lower limits of zero for both the coordinates.The differential transform of the function (, ) is defined as follows [26]: where   /  indicates the th derivative with respect to  and (, ) is the transformed function.The inversion formula for (, ) is given by As seen in ( 2), the original function is expressed by a powerseries function.Substituting (1) into (2) yields which shows that the concept of the differential transform is based on Taylor series expansion (about a point (, ) = (0, 0)).In practice, the function (, ) in ( 2) is approximated with a finite number of terms as where  and  are determined based on the required accuracy of the representation of the original function , which is often a solution to a physical problem.In most practical cases, exact solutions are not available, and thus it is impossible to compute the absolute error |(  ,   ) − ỹ, (  ,   )| for an arbitrary point (  ,   ) in the time-space domain.However, the convergence behaviour of the function ỹ, to the original function  can be shown [27].In order to achieve a convergence of the function/solution with smaller values of  and , the domain split technique developed by Jang et al. [28] is adopted in our numerical calculations (see Section 4).Table 1 lists some of the fundamental operations (similarity, linearity, differentiation, multiplication, etc.) of the 2D differential transform [26].These mathematical operations are easily proved with the aid of ( 1) and (2).In Table 1, (, ), (, ), and (, ) are the transformed functions of original functions (, ), (, ), and (, ), respectively.Using the relations listed in Table 1, the differential transform turns differential equations into recursive polynomial equations, which are much easier to solve.

Heat Conduction Problem with Phase Change
Let us consider a slab with finite thickness  as shown in Figure 1.The surfaces of the slab are subjected to a prescribed temperature   at the right side, which is lower than the fusion temperature   and convective thermal loading with heat transfer coefficient ℎ at the left side.The ambient temperature is given by  ∞ .The temporal variable  indicates the location of the change-of-phase front, at which the temperature is constant (=   ).The slab, which is homogeneous and isotropic with temperature-dependent thermophysical properties, is initially at a uniform temperature  0 ( 0 ≥   ).
In the present paper, the heat conduction in the slab is analysed on the basis of the apparent specific heat method [25].This formulation allows for a continuous treatment of a system undergoing phase transition because the latent heat effect of the phase change is included in the apparent specific heat of the substance.Moreover, the phase change of a substance in which the latent heat is evolved over a range of temperatures can be treated as well as that with a sharp phase change temperature.The former includes liquid solutions, biological tissues, and alloys, whereas the latter is represented  by water.The limitation of the apparent specific heat method in terms of analysis accuracy was discussed by Civan and Sliepcevich [29].
If the volume change during the phase change is neglected (i.e., the density is assumed to be constant), then the heat conduction problem for the slab is formulated as follows: where  is the space variable,  is the time,  is the density,  is the apparent specific heat, and  is the thermal conductivity.
The apparent specific heat of the slab is considered to exhibit the following dependence on the temperature: where  1 ,  2 ,  ∞ ,  (≤ 0), and  1 are constants.Figure 2(a) graphically represents (6), in which the normal distribution function is used for a smooth variation in the apparent specific heat [25].The relationship between the latent heat and the parameters in ( 6) is described in Appendix A. The thermal conductivity of the slab is assumed to satisfy the relationship where  ∞ ,  2 , and  3 are constants.Figure 2(b) describes this temperature dependency.Substituting ( 6) and ( 7) into (5) gives the following equations in dimensionless form:

Thermal conductivity of liquid
Thermal conductivity of solid Equation ( 7)

𝜃 = 𝜃
where Here, it should be noted that the initial condition of ( 9) and boundary condition of ( 11) are discontinuous at (, ) = (0, 1) [28].However, continuity is required at any points in the domain of interest to use the DTM.Thus, the initial condition of ( 9) is modified in the form of a polynomial function, as proposed by Jang et al. [28], where N is a positive integer that is identical with that appearing in (4).Equation ( 8) can be expanded as in which Taking the 2D differential transform of each term in (13) and rearranging the terms yield the following recurrence relation: where  and  denote the dimensionless time and space intervals of interest, respectively, and Θ(, ) is the transformed function of (, ).Functions Ψ(, ) and Φ  (, ) for  = 1,2 are the transformed functions of () and   () defined by ( 14) and ( 15), which are obtained using the procedure developed by S. Chang and I. Chang [30] as follows: where () represents the Kronecker delta, which is equal to 1 if  = 0 and 0 otherwise.Additionally, applying the 2D differential transform to (9  ), (10), and (11) yields the initial and boundary conditions as follows: where the following relation must be fulfilled because of the need for continuity between the initial condition and the boundary condition of ( 10) at (, ) = (0, 0): Finally, we can obtain the following solution by the differential inverse transform of Θ(, ): It should be noted that this series solution is valid for  ∈ [0, ] and  ∈ [0, ].Functions Θ(, ), Ψ(, ), and Φ  (, ), where  = 1,2, are readily obtained through a simple recursive procedure using ( 16)-( 17), which is initiated using Θ(0, ) for  = 0, 1, 2, . . ., .

Verification of Present Solution.
To verify the presented series solution, the transient heat conduction in a finite slab with temperature-independent material properties is considered, for which an analytical solution exists [31].A detailed description of the analytical solution can be found in Appendix B. To improve the accuracy and rate of convergence of the series solution, we divide the full - domain into subdomains, as Jang et al. did [28].The overall  domain is split equally into several subdomains, depending on the numerical examples treated, whereas the time domain is split into many subdomains with interval  = 0.001.The number of power-series terms in ( 20) is determined according to the convergence of the solution.
The temperature distributions calculated from the presented series solution and the analytical solution [31] are shown in Figure 3 for different Biot numbers.The values of the parameters used for this computation are  1 =  2 =  2 = 0,  3 = 0.5,  0 = 1, and   = −1.The series solution with the -domain division into three subdomains (i.e.,  = 1/3) is sufficiently converged for (, ) = (40, 6) (Table 2).It is observed that the converged series solution θ40,6 (, ) is in good agreement with the analytical solution.

Parametric Studies.
All the numerical results shown in this subsection are obtained using the following parameter set:  = 0 and   = −1.Figure 4 shows transient temperature distributions in the slab with different temperature dependence behaviours of the specific heat.Since the specific heat of the liquid phase increases with an increase in the value of  2 (see Figure 2(a)), a higher  2 value reduces the rate of decline of the temperature.When the specific heat is independent of the temperature, that is,  2 = 0, a significant decrease in the temperature is observed due to the specific heat remaining low.It is also observed that the gap between the temperature profiles for different  2 values grows exponentially with time.
Figure 5 illustrates the temperature distributions in the slab with different dependences of the thermal conductivity on the temperature.An increase in the value of  3 raises the levels of the thermal conductivity for both the liquid and solid phases (see Figure 2(b)).Hence, a higher  3 value leads to higher rates of temperature reduction.
Figure 6 shows the temperature distributions in the slab with different latent heat quantities for a constant thermal conductivity and specific heat, that is, the same thermophysical properties between the solid and liquid phases.The latent

Concluding Remarks
The two-dimensional differential transform method has been employed to solve the nonlinear heat conduction problem with phase change in a finite-thickness slab.The slab was subjected to convective thermal loading at one boundary surface and a constant prescribed temperature at the other boundary surface.In addition, the slab had temperaturedependent thermophysical properties: the thermal conductivity and specific heat (or volumetric heat capacity).The treatment of the phase change was based on the apparent specific heat method.The presented analytical method gives an analytical solution in the form of a power series with easily computable components.Numerical results illustrated that the DTM is useful as a new analytical method for solving the phase change problem in a slab with temperature-dependent parameters.Our future plans are (i) to apply the present analytical method to the cases of other geometries (e.g., inward and outward solidification in a cylindrical or spherical geometry) and (ii) to adopt a novel treatment for multiple complex nonlinear terms, as given by ( 14) and ( 15), for the DTM [32], which will improve the computational efficiency.

A. Amount of Latent Heat
The relationship between the latent heat  and the apparent specific heat (in particular, parameters included in ( 6)) is where this integral is equivalent to the area of domain bounded by two smooth solid curves in Figure 2(a).Considering that parameter  is a nonpositive constant, the above integral (i.e., the Gaussian integral) can be evaluated as follows:

B. Analytical Solution for Linear Case
Consider a slab of finite thickness  which is initially at temperature  0 .For times  > 0, there is heat dissipation by convection from its boundary surface at  = 0 into a surrounding at temperature  ∞ .The heat transfer coefficient is denoted by ℎ.

Figure 1 :
Figure 1: Analytical model of one-dimensional phase change process for a finite-thickness slab.

Figure 2 :
Figure 2: Temperature dependencies of (a) apparent specific heat and (b) thermal conductivity.

Table 1 :
Fundamental operations of two-dimensional differential transform.