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

Authors: R. S. Damor; Sushil Kumar; A. K. Shukla
Journal: International Journal of Engineering Mathematics
 (2013)

---

## 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.


---

## Body

## 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α∈(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 α.


### 1.1. Used Fractional Derivative

Fractional derivative is denoted asaDtαf(t); the subscripts a and t 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)Dtαaf(t)={dmdtm[1Γ(m-α)∫at‍f(τ)(t-τ)α+1-m],form-1<α<m,dm(f(t))dtm,forα=m.


### 1.3. Shifted Grunwald Approximation [3]

Consider(2)∂αT(x,t)∂xα=1Γ(-α)limM→∞1hα×∑k=0M‍Γ(k-α)Γ(k+1)T(x-(k-1)h,t),∂αT(x,t)∂xα≈1hα∑k=0i+1‍(gk)Ti-k+1n,wheregk=1Γ(-α)Γ(k-α)Γ(k+1).


## 1.1. Used Fractional Derivative

Fractional derivative is denoted asaDtαf(t); the subscripts a and t 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)Dtαaf(t)={dmdtm[1Γ(m-α)∫at‍f(τ)(t-τ)α+1-m],form-1<α<m,dm(f(t))dtm,forα=m.


## 1.3. Shifted Grunwald Approximation [3]

Consider(2)∂αT(x,t)∂xα=1Γ(-α)limM→∞1hα×∑k=0M‍Γ(k-α)Γ(k+1)T(x-(k-1)h,t),∂αT(x,t)∂xα≈1hα∑k=0i+1‍(gk)Ti-k+1n,wheregk=1Γ(-α)Γ(k-α)Γ(k+1).


## 2. Model Description

One-dimensional geometry is considered as shown in Figure1. A molten material with initial temperature T0 and uniformly distributed volumetric heat generation q is confined between two surfaces x=0 and x=l. At x=0, convecting cooling is applied to the mould at temperature T∞, and the surface x=l is taken as adiabatic.Figure 1
Schematic of the physical problem.


## 3. Mathematical Model

### 3.1. Governing Equation

We propose the governing fractional diffusion equation model for alloy solidification as follows.In solid region(3)ρscs∂Ts∂t=ks∂αTs∂xα+q,1<α≤2,0≤x≤S(t).
In liquid region
(4)ρlcl∂Tl∂t=kl∂αTl∂xα+q,1<α≤2,S(t)≤x≤l.
The conditions at the moving interface x=S(t) are given as
(5)ks∂α-1Ts∂xα-1-kl∂α-1Tl∂xα-1=ρLdS(t)dtatx=S(t),(6)Ts=Tl=Tphatx=S(t),
where ρ, c, k, T, t, and x represent density, specific heat, thermal conductivity, temperature, time, and distance, respectively. Moreover, S(t), q, and L are denoted for position of phase change interface, volumetric heat generation, and latent heat, respectively. The subscript s, l, and ph are concerned for the solid state, liquid state, and phase change, respectively.On making the use of enthalpy(7)H(T)=∫T0T‍cdT,
where T0 is the reference temperature, (3) to (5) become
(8)ρ∂H∂t=K(x,t)∂αT(x,t)∂xα+q,0≤x≤l.
The relation between enthalpy and temperature is given as [17, 18]
(9)H(T)={cs(T-Tms),forT<Tms,(T-Tms)(cs+cl2+LTml-Tms)forTms≤T≤Tml,L+12(cs+cl)+cl(T-Tml),forT>Tml,
and thermal conductivity is given as
(10)K={ksforT<Tms,ks+kl2forTms≤T≤Tml,kl,forT>Tml,
where Tml and Tms are liquidus and solidus temperature, respectively.


### 3.2. Initial and Boundary Conditions

(i) Initially, the system is at 255°C; that is,(11)T(x,t)=255°C,att=0.(ii) Convective boundary condition is atx=0; that is,
(12)K∂T(x,t)∂x=h(T(x,t)-T∞),atx=0,
where h is the heat transfer coefficient and T∞ is the temperature of the cooling media.(iii) Adiabatic condition is atx=l; that is,
(13)∂T(x,t)∂x=0,atx=l.


## 3.1. Governing Equation

We propose the governing fractional diffusion equation model for alloy solidification as follows.In solid region(3)ρscs∂Ts∂t=ks∂αTs∂xα+q,1<α≤2,0≤x≤S(t).
In liquid region
(4)ρlcl∂Tl∂t=kl∂αTl∂xα+q,1<α≤2,S(t)≤x≤l.
The conditions at the moving interface x=S(t) are given as
(5)ks∂α-1Ts∂xα-1-kl∂α-1Tl∂xα-1=ρLdS(t)dtatx=S(t),(6)Ts=Tl=Tphatx=S(t),
where ρ, c, k, T, t, and x represent density, specific heat, thermal conductivity, temperature, time, and distance, respectively. Moreover, S(t), q, and L are denoted for position of phase change interface, volumetric heat generation, and latent heat, respectively. The subscript s, l, and ph are concerned for the solid state, liquid state, and phase change, respectively.On making the use of enthalpy(7)H(T)=∫T0T‍cdT,
where T0 is the reference temperature, (3) to (5) become
(8)ρ∂H∂t=K(x,t)∂αT(x,t)∂xα+q,0≤x≤l.
The relation between enthalpy and temperature is given as [17, 18]
(9)H(T)={cs(T-Tms),forT<Tms,(T-Tms)(cs+cl2+LTml-Tms)forTms≤T≤Tml,L+12(cs+cl)+cl(T-Tml),forT>Tml,
and thermal conductivity is given as
(10)K={ksforT<Tms,ks+kl2forTms≤T≤Tml,kl,forT>Tml,
where Tml and Tms are liquidus and solidus temperature, respectively.


## 3.2. Initial and Boundary Conditions

(i) Initially, the system is at 255°C; that is,(11)T(x,t)=255°C,att=0.(ii) Convective boundary condition is atx=0; that is,
(12)K∂T(x,t)∂x=h(T(x,t)-T∞),atx=0,
where h is the heat transfer coefficient and T∞ is the temperature of the cooling media.(iii) Adiabatic condition is atx=l; that is,
(13)∂T(x,t)∂x=0,atx=l.


## 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
(14)Hin+1=Hin+Δt(Δx)αρinKin∑k=0i+1‍(gk)Ti-k+1n+Δtρinqin,[H]n+1=[H]n+Δt(Δx)α[G]n[T]n+Δt[qρ]n,
where
(15)G(i,j)=0ifj≥i+2,G(i,j)=Kiρig0ifj=i+1,G(i,j)=Kiρig1ifj=i,G(i,j)=Kiρigi-j+1ifj≤i-1.
Here i and n are space and time steps, respectively. Space and time increments are denoted by Δx and Δt, respectively. These are adjustable subjects to satisfy the stability criteria, as Δt/(Δx)α≤1/dmaxα, where dmax is the maximum value for thermal diffusivity. For α=2, this stability criterion is reduced to Δt/(Δx)2≤1/dmax2 [18]. On calculating the enthalpy at (n+1)th level using (14), temperature distribution at (n+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.


## 5.  Results and Discussion

Thermophysical properties of the Sn-5 wt% Pb alloy used for simulations are enlisted in Table1 [21]. To study the transient behaviour of the system, the whole process of solidification is divided into four different stages as in [2, 22].Table 1
Thermophysical properties of Sn-5 wt%Pb [21].
Properties
Value
Thermal conductivity of solid (W/m°C)
65.6
Thermal conductivity of liquid (W/m°C)
32.8
Density of solid (kg/m3)
7475
Density of liquid (kg/m3)
7181
Specific heat of solid (J/kg°C)
217
Specific heat of liquid (J/kg°C)
253
Latent heat (KJ/kg)
592.14
Solidus temperature (°C)
183
Liquidus temperature (°C)
224.9In stage 1, the entire region is in liquid state; this stage ends when the mushy region starts atx=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 x=0; 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 q=500 W/m3 and h=1500 W/m2 and observe the temperature profiles, liquidus interface, solidus interface, and rates of liquidus, and solidus interface for different values of α.Figure2 shows the temperature profile for α=2 at t=1 second, t=100 seconds, t=460 seconds, t=679 seconds, t=1000 seconds, t=2000 seconds, t=3000 seconds, t=3779 seconds, t=3850 seconds, and t=4265 seconds. Stage 1 ends at t=1 second when the mushy state starts at x=0. At t=100 second the system is in stage 2; that is, the mushy and liquid phases coexist. Stage 2 ends at t=679 seconds, while stage 3 starts at t=460 seconds. Here stage 3 starts before the end of stage 2. At t=1000 seconds, there exist solid mushy and liquid states. At t=3000 seconds, the solid and mushy phases coexist; that is, the system is in stage 3. Stage 3 ends at t=3779 seconds, when the entire region is solidified. At t=3850 seconds, the system is now in stage 4, which ends at t=4265 seconds when the steady state is obtained.Figure 2
Temperature profile in region at different timet for α=2.Figure3 shows the temperature profile for α=1.9 at t=0.3 seconds, t=36 seconds, t=100 seconds, t=277 seconds, t=800 seconds, t=1100 seconds, t=1378 seconds, t=1420 seconds, t=1500 seconds, and t=1812 seconds. Here, stage 1 ends at t=0.3 seconds, when the mushy state starts at x=0. At t=100 seconds, the system turns in stage 2; that is, the region where mushy and liquid phases coexist. Stage 2 ends at t=277 seconds, and the process of stage 3 starts at t=36 seconds. Here stage 3 starts before the end of stage 2. At t=800 seconds, system is in solid-mushy-liquid states. At t=1100 seconds, the solid and mushy phases coexist; that is, the system enters in stage 3. Stage 3 ends at t=1378seconds, when the entire region converts into the solid state. At t=1420 seconds and t=1500 seconds, the system is now in stage 4, which ends when the steady state is obtained at t=1812 seconds.Figure 3
Temperature profile in region at different timet for α=1.9.Figure4 shows the temperature profile for α=1.8 at t=0.3 seconds, t=34 seconds, t=100 seconds, t=274 seconds, t=1000 seconds, t=1500 seconds, t=2084 seconds, t=2120 seconds, t=2200 seconds, and t=2703 seconds. Stage 1 ends at t=0.3 seconds, when the mushy state starts at x=0. At t=100 seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at t=274 seconds, while stage 3 starts at t=34 seconds. Here the stage 3 starts before the end of stage 2. At t=1000 seconds, there exist solid mushy and liquid states. At t=1500 seconds, we observe the solid and mushy phases; that is, the system is now at stage 3. Stage 3 ends at t=2084 seconds, when solidification takes place in the entire region. At t=2120 seconds and t=2200 seconds, the system enters in stage 4, which ends when it attains its steady state at t=2703 seconds.Figure 4
Temperature profile in region at different timet for α=1.8.Figure5 shows the temperature profile for α=1.7 at t=0.5 seconds, t=48 seconds, t=200 seconds, t=537 seconds, t=1500 seconds, t=2500 seconds, t=3805 seconds, t=3900 seconds, t=4000 seconds, and t=4739 seconds. Stage 1 ends at t=0.5 seconds, when the mushy state starts at x=0. At t=200 seconds, the system enters into stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at t=537 seconds, while stage 3 starts at t=48 seconds. Here stage 3 starts before the end of stage 2. At t=1500 seconds, and there exist solid mushy and liquid states. At t=2500 seconds, we observe the solid and mushy phases, that is, the system is in stage 3. Stage 3 ends at t=3805 seconds, when the solidification takes place in the entire region. At t=3900 seconds and t=4000 seconds, the system is in stage 4; it ends when it attains its steady state at t=4739 seconds.Figure 5
Temperature profile in region at different timet for α=1.7.Figure6 shows the temperature profile for α=1.6 at t=1second, t=81 seconds, t=500 seconds, t=1000 seconds, t=1804 seconds, t=4000 seconds, t=7000 seconds, t=8837 seconds, t=9100 seconds, and t=10247 seconds. Stage 1 ends at t=1second, when the mushy state starts at x=0. At t=500 seconds, the system is in stage 2; that is, the coexistence of mushy and liquid phases. Stage 2 ends at t=1804 seconds, and process of stage 3 starts at t=81 seconds. Here stage 3 starts before the end of stage 2. At t=1000 seconds, and there exist solid mushy and liquid states. At t=4000 seconds and t=7000 seconds, we observe the solid and mushy phases; that is, the system is in stage 3. Stage 3 ends at t=8837 seconds, when the solidification takes place in the entire region. At t=9100second the system is in stage 4; it ends when it attains its steady state at t=10247 seconds. Table 2 represents the comparison of the ends of different stages for different values of α.Table 2
Time taken to the end of the stages for differentα.
Alpha
Stage 1
Stage 2
Stage 3
Stage 4
α
=
2
t
=
1 second
t
=
679 seconds
t
=
3779 seconds
t
=
3850 seconds
α
=
1.9
t
=
0.3 seconds
t
=
277 seconds
t
=
1378 seconds
t
=
1812 seconds
α
=
1.8
t
=
0.3 seconds
t
=
274 seconds
t
=
2084 seconds
t
=
2703 seconds
α
=
1.7
t
=
0.5 seconds
t
=
537 seconds
t
=
3805 seconds
t
=
4739 seconds
α
=
1.6
t
=
1 second
t
=
1804 seconds
t
=
8837 seconds
t
=
10247 secondsFigure 6
Temperature profile in region at different timet for α=1.6.Figures7, 8, and 9 exhibit the temperature profile for time t=1 second, t=36 seconds, t=177 seconds, t=1378 seconds, and t=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 7
Temperature profile for differentα at t=1second and t=36 seconds.Figure 8
Temperature profile for differentα at t=277 seconds and t=1378 seconds.Figure 9
Temperature profile for differentα at t=1812 seconds.Figure10 shows the position of liquidus interface for different values of α. We see that interface reaches the end point x=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 x=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.Figure 10
Liquidus interface for differentα.Figure 11
Solidus interface for differentα.Figures12 and 13 exhibit the velocity of liquidus and solidus interface, respectively. We observe that the velocity of interfaces increases as α increases.Figure 12
Velocity of liquidus interface at interior point for differentα.Figure 13
Velocity of solidus interface at interior point for differentα.


## 6. 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.



---

*Source: 785609-2013-09-10.xml*