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

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.


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 2 International Journal of Engineering Mathematics 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 ∈ (1, 2] 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 .

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. [3]. Consider

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

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 ph are concerned for the solid state, liquid state, and phase change, respectively.
On making the use of enthalpy where 0 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.

Initial and Boundary
Conditions. (i) Initially, the system is at 255 ∘ C; that is, (ii) Convective boundary condition is at = 0; that is, where ℎ is the heat transfer coefficient and ∞ is the temperature of the cooling media.
(iii) Adiabatic condition is at = ; that is,

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 Δ /(Δ ) ≤ 1/ max , where max is the maximum value for thermal diffusivity. For = 2, this stability criterion is reduced to Δ /(Δ ) 2 ≤ 1/ max 2 [18]. On calculating the enthalpy at ( +1)th level using (14), temperature distribution at ( + 1)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.

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 = 0. 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 = 0; during this stage, the entire region may be in the solid-mushy or solidmushy-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 = 500 W/m 3 and ℎ = 1500 W/m 2 and observe the temperature profiles, liquidus interface, solidus interface, and rates of liquidus, and solidus interface for different values of .   = 100 second the system is in stage 2; that is, the mushy and liquid phases coexist. Stage 2 ends at = 679 seconds, while stage 3 starts at = 460 seconds. Here stage 3 starts before the end of stage 2. At = 1000 seconds, there exist solid mushy and liquid states. At = 3000 seconds, the solid and mushy phases coexist; that is, the system is in stage 3. Stage 3 ends at = 3779 seconds, when the entire region is solidified. At = 3850 seconds, the system is now in stage 4, which ends at = 4265 seconds when the steady state is obtained. Figure 3 shows the temperature profile for = 1.9 at = 0.3 seconds, = 36 seconds, = 100 seconds, = 277 seconds, = 800 seconds, = 1100 seconds, = 1378 seconds, = 1420 seconds, = 1500 seconds, and = 1812 seconds. Here, stage 1 ends at = 0.3 seconds, when the mushy state starts at = 0. At = 100 seconds, the system turns in stage 2; that is, the region where mushy and liquid phases coexist. Stage 2 ends at = 277 seconds, and the process of stage 3 starts at = 36 seconds. Here stage 3 starts before the end of stage 2. At = 800 seconds, system is in solid-mushy-liquid states. At = 1100 seconds, the solid and mushy phases coexist; that is, the system enters in stage 3. Stage 3 ends at = 1378 seconds, when the entire region converts into the solid state. At = 1420 seconds and = 1500 seconds, the system is now in stage 4, which ends when the steady state is obtained at = 1812 seconds. Figure 4 shows the temperature profile for = 1.   At = 100 seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at = 274 seconds, while stage 3 starts at = 34 seconds. Here the stage 3 starts before the end of stage 2. At = 1000 seconds, there exist solid mushy and liquid states. At = 1500 seconds, we observe the solid and mushy phases; that is, the system is now at stage 3. Stage 3 ends at = 2084 seconds, when solidification takes place in the entire region. At = 2120 seconds and = 2200 seconds, the system enters in stage 4, which ends when it attains its steady state at = 2703 seconds.      ends at = 0.5 seconds, when the mushy state starts at = 0. At = 200 seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at = 537 seconds, while stage 3 starts at = 48 seconds. Here stage 3 starts before the end of stage 2. At = 1500 seconds, and there exist solid mushy and liquid states. At = 2500 seconds, we observe the solid and mushy phases, that is, the system is in stage 3. Stage 3 ends at = 3805 seconds, when the solidification takes place in the entire region. At = 3900 seconds and = 4000 seconds, the system is in stage 4; it ends when it attains its steady state at = 4739 seconds. Figure 6 shows the temperature profile for = 1.6 at = 1 second, = 81 seconds, = 500 seconds, = 1000 seconds, = 1804 seconds, = 4000 seconds, = 7000 seconds, = 8837 seconds, = 9100 seconds, and = 10247 seconds. Stage 1 ends at = 1 second, when the mushy state starts at = 0. At = 500 seconds, the system is in stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at = 1804 seconds, and process of stage 3 starts at = 81 seconds. Here stage 3 starts before the end of stage 2. At = 1000 seconds, and there exist solid mushy and liquid states. At = 4000 seconds and = 7000 seconds, we observe the solid and mushy phases; that is, the system is in stage 3. 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 = 1 second, = 36 seconds, = 177 seconds, = 1378 seconds, and = 1812 seconds, respectively. On comparing the temperature for different values of at fixed time, we observe that temperature is lower in a system for = 1.9 as compared to = 1.8, 1.7, and 1.6. Figure 10 shows the position of liquidus interface for different values of . We see that interface reaches the end point = 0.05 meter for different values of the , and for = 2, 1.9, 1.8, 1.7, and 1.6 interface takes 679 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 = 0.05 meter for different values of as = 2, 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.9, 1.8, and 1.7 as compared to = 2. Figures 12 and 13 exhibit the velocity of liquidus and solidus interface, respectively. We observe that the velocity of interfaces increases as increases.

Conclusion
In this study, we observe the following.
(i) The system takes less time to attain the steady state for = 1.9 and 1.8 as compared to = 2; that is, normal diffusion.
(ii) Solidification starts earlier for = 1.9 and 1.8 as compared to = 2; that is, normal diffusion.
(iii) Solidus and liquidus interfaces take less time to reach and the end point for = 1.9, 1.8, and 1.7 as compared to = 2. (iv) The velocity of the solidus and liquidus is higher at = 1.9 as compared to = 1.8, 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.9, 1.8, and 1.7.