International Journal of Engineering Mathematics

Volume 2013 (2013), Article ID 785609, 8 pages

http://dx.doi.org/10.1155/2013/785609

## Numerical Solution of Fractional Diffusion Equation Model for Freezing in Finite Media

Department of Applied Mathematics and Humanities, S. V. National Institute of Technology, Surat-395 007, India

Received 28 March 2013; Accepted 8 August 2013

Academic Editor: C. Nataraj

Copyright © 2013 R. S. Damor et al. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

#### Abstract

Phase change problems play very important role in engineering sciences including casting of nuclear waste materials, vivo freezing of biological tissues, solar collectors and so forth. In present paper, we propose fractional diffusion equation model for alloy solidification. A transient heat transfer analysis is carried out to study the anomalous diffusion. Finite difference method is used to solve the fractional differential equation model. The temperature profiles, the motion of interface, and interface velocity have been evaluated for space fractional diffusion equation.

#### 1. Introduction

Heat diffusion during the change of state, for example, melting and solidification, has important applications in science and technology including casting of nuclear waste materials, vivo freezing of biological tissues, and solar collectors. Problems involving heat conduction in materials in which a phase change occurs are known as moving boundary problems. First time Stefan [1] discussed this problem to study the melting of polar ice, so this is referred as Stefan problem. Solidification is the phase change problem in which phase transformation takes place by descending the energy of the system, that is, heat extraction from the liquid region. The problem related to the heat conduction involving solidification has drawn the attention of many researchers. Solidification process can be grouped into two major categories, that is, solidification of pure substances and alloy. Solidification of pure substances is usually characterized by a sharp solid or liquid interface. On the other hand, a mixed phase region characterizes solidification of an alloy. The mixed phase regions are commonly termed as the mushy region; this is a combination of liquid solute and solid crystals [2].

Fractals and fractional calculus have been used to improve the modelling accuracy of many phenomena in natural sciences. The most important advantage of using fractional differential equations is their nonlocal property. This means that the next state of a system depends not only upon its current state but also depends upon all of its historical states. These are more realistic and one of reasons to make the fractional calculus more and more popular [3, 4]. The fractional model is appropriate for modelling complex dynamics as the ions are undergoing anomalous diffusion. Fractional order models are more general and useful than the previously used integer order model. Fractional derivative enables the description of the memory and hereditary inherent properties in various processes governed by anomalous diffusion. Fractional diffusion equations are studied by many researchers. Meerschaert et al. [5] gave a second-order accurate numerical approximation for the fractional diffusion equation, Murio [6] discussed implicit finite difference approximation for time fractional diffusion equation, Lijuan et al. [7] gave finite difference-approximation for the fractional advection diffusion equation using Riemann-Liouville space fractional derivative. Fractional differential equation models are used to study anomalous diffusion; when a space fractional derivative is replaced by the second-order derivative in a diffusion model, it leads to super-diffusion [8].

Some researchers have studied simple or classic Stefan problem using fractional diffusion equation. Li et al. [9] gave some exact solution to Stefan problem by adopting Caputo form of fractional time derivative and Dirichlet boundary conditions. Voller [10] obtained an exact solution of a Stefan problem by considering Stefan number as zero. Recently, Atkinson [11] wrote a very interesting survey paper on moving boundary problems for time and space fractional diffusion. They studied the classic Stefan problem and discussed analytical or semianalytical solution in closed form for infinite or semi-infinite region.

Alloy solidification problems are distinct from classic or simple Stefan problem as the melting temperature is not known in advance; this depends on the composition of the alloy [12]. Since in alloy phase change occurs on wide range, the mathematical formulations are highly nonlinear as given by Carslaw and Jaeger [13]. A numerical solution can be obtained using finite difference method. The enthalpy method is based on weak formulation, which was used by many researchers (Alexiades and Solomon [14], Shamsunder and Sparrow [15], Voller [16], Bonacina et al. [17], Hoffman and Bischof [18], Jiji and Gaye [19], and Kumar and Katiyar [2, 20]). The essential features of the basic enthalpy methods are that they do not need tracking of the solid-liquid interface; a fixed grid can be used for numerical computations, and the nonlinearity at the moving boundary can also be avoided.

In present study, the enthalpy formulation of space fractional diffusion model is proposed to study the alloy solidification in finite media. Space fractional derivative of order is considered of the Riemann-Liouville form, and it is approximated by shifted Grunwald formula [3]. Finite difference method is used to obtain numerical solution of the model. To study the effect of fractional order derivative , the temperature profiles and motion of freezing interface are obtained for different values of .

##### 1.1. Used Fractional Derivative

Fractional derivative is denoted as ; the subscripts and denote the two limits related to the operation of fractional differentiation, which is called the terminal of fractional differentiation. The negative values for are denoted for fractional integrals of arbitrary order.

##### 1.2. Riemann-Liouville Fractional Derivative [3]

Consider

##### 1.3. Shifted Grunwald Approximation [3]

Consider

#### 2. Model Description

One-dimensional geometry is considered as shown in Figure 1. A molten material with initial temperature and uniformly distributed volumetric heat generation is confined between two surfaces and . At , convecting cooling is applied to the mould at temperature , and the surface is taken as adiabatic.

#### 3. Mathematical Model

##### 3.1. Governing Equation

We propose the governing fractional diffusion equation model for alloy solidification as follows.

In solid region In liquid region The conditions at the moving interface are given as where , , , , , and represent density, specific heat, thermal conductivity, temperature, time, and distance, respectively. Moreover, , , and are denoted for position of phase change interface, volumetric heat generation, and latent heat, respectively. The subscript , , and are concerned for the solid state, liquid state, and phase change, respectively.

On making the use of enthalpy where is the reference temperature, (3) to (5) become The relation between enthalpy and temperature is given as [17, 18] and thermal conductivity is given as where and are liquidus and solidus temperature, respectively.

##### 3.2. Initial and Boundary Conditions

(i) Initially, the system is at 255°C; that is,

(ii) Convective boundary condition is at ; that is, where is the heat transfer coefficient and is the temperature of the cooling media.

(iii) Adiabatic condition is at ; that is,

#### 4. Numerical Solution

Using forward difference for time derivative and shifted Grunwald approximation [3] for space fractional derivative, finite difference scheme for (8) is given as where Here and are space and time steps, respectively. Space and time increments are denoted by and , respectively. These are adjustable subjects to satisfy the stability criteria, as , where is the maximum value for thermal diffusivity. For , this stability criterion is reduced to [18]. On calculating the enthalpy at th level using (14), temperature distribution at th time level is calculated by reverting (9). Once the new temperature field is obtained from enthalpy, the process repeats until the system reaches steady state. Isotherms at 224.9°C and 183.0°C give the position of liquidus and solidus interface, respectively.

#### 5. Results and Discussion

Thermophysical properties of the Sn-5 wt% Pb alloy used for simulations are enlisted in Table 1 [21]. To study the transient behaviour of the system, the whole process of solidification is divided into four different stages as in [2, 22].

In stage 1, the entire region is in liquid state; this stage ends when the mushy region starts at . During stage 2, there is the coexistence of liquid and mushy regions; this stage ends when the entire region converts into mushy region. Stage 3 begins when solidification takes place from ; during this stage, the entire region may be in the solid-mushy or solid-mushy-liquid state. This stage ends when the entire region is solidified. Only solid phase exists in stage 4, which ends when the system reaches the steady state. In our study, we fix the values of W/m^{3} and W/m^{2} and observe the temperature profiles, liquidus interface, solidus interface, and rates of liquidus, and solidus interface for different values of .

Figure 2 shows the temperature profile for at second, seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, and seconds. Stage 1 ends at second when the mushy state starts at . At second the system is in stage 2; that is, the mushy and liquid phases coexist. Stage 2 ends at seconds, while stage 3 starts at seconds. Here stage 3 starts before the end of stage 2. At seconds, there exist solid mushy and liquid states. At seconds, the solid and mushy phases coexist; that is, the system is in stage 3. Stage 3 ends at seconds, when the entire region is solidified. At seconds, the system is now in stage 4, which ends at seconds when the steady state is obtained.

Figure 3 shows the temperature profile for at seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, and seconds. Here, stage 1 ends at seconds, when the mushy state starts at . At seconds, the system turns in stage 2; that is, the region where mushy and liquid phases coexist. Stage 2 ends at seconds, and the process of stage 3 starts at seconds. Here stage 3 starts before the end of stage 2. At seconds, system is in solid-mushy-liquid states. At seconds, the solid and mushy phases coexist; that is, the system enters in stage 3. Stage 3 ends at seconds, when the entire region converts into the solid state. At seconds and seconds, the system is now in stage 4, which ends when the steady state is obtained at seconds.

Figure 4 shows the temperature profile for at seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, and seconds. Stage 1 ends at seconds, when the mushy state starts at . At seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at seconds, while stage 3 starts at seconds. Here the stage 3 starts before the end of stage 2. At seconds, there exist solid mushy and liquid states. At seconds, we observe the solid and mushy phases; that is, the system is now at stage 3. Stage 3 ends at seconds, when solidification takes place in the entire region. At seconds and seconds, the system enters in stage 4, which ends when it attains its steady state at seconds.

Figure 5 shows the temperature profile for at seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, and seconds. Stage 1 ends at seconds, when the mushy state starts at . At seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at seconds, while stage 3 starts at seconds. Here stage 3 starts before the end of stage 2. At seconds, and there exist solid mushy and liquid states. At seconds, we observe the solid and mushy phases, that is, the system is in stage 3. Stage 3 ends at seconds, when the solidification takes place in the entire region. At seconds and seconds, the system is in stage 4; it ends when it attains its steady state at seconds.

Figure 6 shows the temperature profile for at second, seconds, seconds, seconds, seconds, seconds, seconds, seconds, seconds, and seconds. Stage 1 ends at second, when the mushy state starts at . At seconds, the system is in stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at seconds, and process of stage 3 starts at seconds. Here stage 3 starts before the end of stage 2. At seconds, and there exist solid mushy and liquid states. At seconds and seconds, we observe the solid and mushy phases; that is, the system is in stage 3. Stage 3 ends at seconds, when the solidification takes place in the entire region. At second the system is in stage 4; it ends when it attains its steady state at seconds. Table 2 represents the comparison of the ends of different stages for different values of .

Figures 7, 8, and 9 exhibit the temperature profile for time second, seconds, seconds, seconds, and seconds, respectively. On comparing the temperature for different values of at fixed time, we observe that temperature is lower in a system for as compared to , and .

Figure 10 shows the position of liquidus interface for different values of . We see that interface reaches the end point meter for different values of the , and for , and interface takes seconds, 277 seconds, 274 seconds, 537 seconds, and 1804 seconds, respectively. Figure 11 shows the position of solidus interface for different values of . This takes 3779 seconds, 1378 seconds, 2084 seconds, 3805 seconds, and 8837 seconds, respectively, to reach interface at the end point meter for different values of as , 1.9, 1.8, 1.7, and 1.6. It has been observed that interfaces take less time to reach the end point for , 1.8, and 1.7 as compared to .

Figures 12 and 13 exhibit the velocity of liquidus and solidus interface, respectively. We observe that the velocity of interfaces increases as increases.

#### 6. Conclusion

In this study, we observe the following.(i)The system takes less time to attain the steady state for and 1.8 as compared to ; that is, normal diffusion.(ii)Solidification starts earlier for and 1.8 as compared to ; that is, normal diffusion.(iii)Solidus and liquidus interfaces take less time to reach and the end point for , 1.8, and 1.7 as compared to .(iv)The velocity of the solidus and liquidus is higher at as compared to , 1.7, and 1.6.

Looking at the above observations, we can say that the normal diffusion takes more time to attain the steady state position whereas liquidus and solidus interfaces also take more time to reach the end point as compared to fractional diffusion model for , 1.8, and 1.7.

#### Nomenclature

Specific heat () | |

: | Enthalpy of tissue () |

: | Heat transfer coefficient (W/(m^{2°}C)) |

: | Thermal conductivity of the tissue () |

: | Latent heat () |

: | Volumetric heat generation () |

: | Position of phase change interface (m) |

: | Solid state |

: | Temperature (°C) |

: | Time (second) |

: | Space coordinate (m) |

: | Maximum thermal diffusivity () |

: | Increment in space (m) and time (second) |

: | Density (). |

#### Acknowledgments

The authors are thankful to the reviewers for their valuable suggestions for improvement of the paper. The first author is thankful to Technical Education Department, Government of Gujarat, India, for providing the opportunity to carry out the research work under quality improvement programme at SVNIT, Surat, India.

#### References

- J. Stefan, “Ueber die theorie der eisbildung, insbesondere ber die eisbil—dung im polarmeere,”
*Annalen Der Physik Und Chemie*, vol. 42, pp. 269–286, 1891. View at Google Scholar - S. Kumar and V. K. Katiyar, “Transient analysis of alloy freezing in finite media with energy generation and convective cooling,”
*International Journal of Applied Mechanics and Engineering*, vol. 15, no. 4, pp. 1155–1168, 2010. View at Google Scholar - I. Podlubny,
*Fractional Differential Equations*, Academic Press, New York, NY, USA, 1999. - F. Mainardi, “Fractional calculus: some basic problems in continuum and statistical mechanics,” in
*Fractals and Fractional Calculus in Continuum Mechanics*, A. Carpinteri and F. Mainardi, Eds., pp. 291–348, Springer, New York, NY, USA, 1997. View at Google Scholar - M. M. Meerschaert, H. P. Scheffler, and C. Tadjeran, “Finite difference methods for two-dimensional fractional dispersion equation,”
*Journal of Computational Physics*, vol. 211, no. 1, pp. 249–261, 2006. View at Publisher · View at Google Scholar · View at Scopus - D. A. Murio, “Implicit finite difference approximation for time fractional diffusion equations,”
*Computers and Mathematics with Applications*, vol. 56, no. 4, pp. 1138–1145, 2008. View at Publisher · View at Google Scholar · View at Scopus - S. Lijuan, W. Wenqia, and Y. Zhoxima, “Finite difference approximations for the fractional advection-diffusion equation,”
*Physics Letters A*, vol. 373, no. 48, pp. 4405–4408, 2009. View at Publisher · View at Google Scholar · View at Scopus - M. M. Meerschaert and C. Tadjeran, “Finite difference approximations for fractional advection-dispersion flow equations,”
*Journal of Computational and Applied Mathematics*, vol. 172, no. 1, pp. 65–77, 2004. View at Publisher · View at Google Scholar · View at Scopus - X. Li, M. Xu, and X. Jiang, “Homotopy perturbation method to time-fractional diffusion equation with a moving boundary condition,”
*Applied Mathematics and Computation*, vol. 208, no. 2, pp. 434–439, 2009. View at Publisher · View at Google Scholar · View at Scopus - V. R. Voller, “An exact solution of a limit case Stefan problem governed by a fractional diffusion equation,”
*International Journal of Heat and Mass Transfer*, vol. 53, no. 23-24, pp. 5622–5625, 2010. View at Publisher · View at Google Scholar · View at Scopus - C. Atkinson, “Moving boundary problems for time fractional and composition dependent diffusion,”
*Fractional Calculus and Applied Analysis*, vol. 15, pp. 207–221, 2012. View at Google Scholar - J. Crank,
*Free and Moving Boundary Problems*, Claredon Press, Oxford, UK, 1984. - H. S. Carslaw and J. C. Jaeger,
*Conduction of Heat in Solids*, Oxford University Press, London, UK, 2nd edition, 1959. - V. Alexiades and A. D. Solomon,
*Mathematical Modelling of Melting and Freezing Process*, Hemisphere, Washington, DC, USA, 1993. - N. Shamsundar and E. M. Sparrow, “Analysis of multidimensional conduction phase change via the enthalpy model,”
*ASME Journal of Heat Transfer*, vol. 97, no. 3, pp. 333–340, 1975. View at Google Scholar · View at Scopus - V. R. Voller, “An implicit enthalpy solution for phase change problems: with application to a binary alloy solidification,”
*Applied Mathematical Modelling*, vol. 11, no. 2, pp. 110–116, 1987. View at Google Scholar · View at Scopus - C. Bonacina, G. Comini, A. Fasano, and M. Primicerio, “Numerical solution of phase-change problems,”
*International Journal of Heat and Mass Transfer*, vol. 16, no. 10, pp. 1825–1832, 1973. View at Google Scholar · View at Scopus - N. E. Hoffmann and J. C. Bischof, “Cryosurgery of normal and tumor tissue in the dorsal skin flap chamber: part I—thermal response,”
*ASME Journal of Biomechanical Engineering*, vol. 123, no. 4, pp. 301–309, 2001. View at Publisher · View at Google Scholar · View at Scopus - L. M. Jiji and S. Gaye, “Analysis of solidification and melting of PCM with energy generation,”
*Applied Thermal Engineering*, vol. 26, no. 5-6, pp. 568–575, 2006. View at Publisher · View at Google Scholar · View at Scopus - S. Kumar and V. K. Katiyar, “Numerical study on phase change heat transfer during combined hyperthermia and cryosurgical treatment of lung cancer,”
*International Journal of Applied Mathematics and Mechanics.*, vol. 3, no. 3, pp. 1–17, 2007. View at Google Scholar - I. L. Ferreira, J. E. Spinelli, J. C. Pires, and A. Garcia, “The effect of melt temperature profile on the transient metal/mold heat transfer coefficient during solidification,”
*Materials Science and Engineering A*, vol. 408, no. 1-2, pp. 317–325, 2005. View at Publisher · View at Google Scholar · View at Scopus - V. K. Katiyar and B. Mohanty, “Transient heat transfer analysis for moving boundary transport problems in finite media,”
*International Journal of Heat and Fluid Flow*, vol. 10, no. 1, pp. 28–31, 1989. View at Google Scholar · View at Scopus