Abstract
A two-dimensional differential transform method is applied to solve one-dimensional phase change problems in a slab of finite thickness, which is subjected to convective thermal loading at one surface and a constant prescribed temperature at the other. In the problems, the initial temperature of the slab does not necessarily have to be the same as the fusion temperature. A series solution is derived for the temperature profile in the melting or solidifying slab with temperature-dependent thermal conductivity and volumetric heat capacity. The latent heat effect of the phase change is incorporated into the temperature-dependent heat capacity. Numerical results demonstrate the effects of the temperature-dependent parameters on the transient temperature profile of the slab.
1. 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–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 semi-infinite 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 temperature-dependent 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.
2. 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 power-series function. Substituting (1) into (2) yields which shows that the concept of the differential transform is based on Taylor series expansion (about a point ). 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.
3. 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 ().
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 , , , (), and 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 , , and are constants. Figure 2(b) describes this temperature dependency.
(a)
(b)
Substituting (6) and (7) into (5) gives the following equations in dimensionless form: where Here, it should be noted that the initial condition of (9) and boundary condition of (11) are discontinuous at [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 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 and 0 otherwise.
Additionally, applying the 2D differential transform to , (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 : Finally, we can obtain the following solution by the differential inverse transform of : It should be noted that this series solution is valid for and . Functions , , and , where , are readily obtained through a simple recursive procedure using (16)–(17), which is initiated using for .
4. Numerical Results and Discussion
4.1. 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 . 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 , , , and . The series solution with the -domain division into three subdomains (i.e., ) is sufficiently converged for (Table 2). It is observed that the converged series solution is in good agreement with the analytical solution.
(a)
(b)
4.2. Parametric Studies
All the numerical results shown in this subsection are obtained using the following parameter set: and . 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 (see Figure 2(a)), a higher value reduces the rate of decline of the temperature. When the specific heat is independent of the temperature, that is, , 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 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 raises the levels of the thermal conductivity for both the liquid and solid phases (see Figure 2(b)). Hence, a higher 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 heat required for the phase change completion from the initial state, , is measured here by the Stefan number . In the case of , (6) yields . It is observed in Figure 6 that an increase in the Stefan number accelerates the temperature reduction. Because there is an inverse relationship between the latent heat and Stefan number (see the definition above), the latent heat flux is dominant over the sensible heat flux in the case of small Stefan numbers and consequently the temperature variation is suppressed. This means that a larger Stefan number results in a higher solid phase development rate. From Figures 4–6, one can see that the series solution derived provides reasonable results.
5. 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 temperature-dependent 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.
Appendices
A. Amount of Latent Heat
The relationship between the latent heat and the apparent specific heat (in particular, parameters included in 6) is described here. The latent heat is given by the integral of the normal distribution function in 6 over the entire range of temperature 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 . For times , there is heat dissipation by convection from its boundary surface at into a surrounding at temperature . The heat transfer coefficient is denoted by . The other boundary surface at is kept at temperature . An analytical solution for the present boundary-value problem is obtained as where and are the thermal diffusivity and thermal conductivity, respectively, and Eigenvalues are the positive roots of the transcendental equation
Conflict of Interests
The author declares that there is no conflict of interests regarding the publication of this paper.