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

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 𝛼∈(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 _𝑎𝐷^𝛼_𝑡𝑓(𝑡); 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Γ()∫𝑚−𝛼^𝑡𝑎𝑓(𝜏)(𝑡−𝜏)^{𝛼+1−𝑚}𝑑,for𝑚−1<𝛼<𝑚,_𝑚𝑓(𝑡)𝑑𝑡^𝑚,for𝛼=𝑚. (1)

1.3. Shifted Grunwald Approximation [3]

Consider 𝜕^𝛼𝑇(𝑥,𝑡)𝜕𝑥^𝛼=1Γ(−𝛼)lim_𝑀→∞1ℎ^𝛼×^𝑀_𝑘=0Γ(𝑘−𝛼)Γ()𝑇(()),𝜕𝑘+1𝑥−𝑘−1ℎ,𝑡_𝛼𝑇(𝑥,𝑡)𝜕𝑥^𝛼≈1ℎ^𝛼𝑖+1_𝑘=0𝑔_𝑘𝑇^𝑛_{𝑖−𝑘+1},where𝑔_𝑘=1Γ()Γ()−𝛼𝑘−𝛼Γ().𝑘+1 (2)

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 𝑥=0 and 𝑥=𝑙. At 𝑥=0, convecting cooling is applied to the mould at temperature 𝑇_∞, and the surface 𝑥=𝑙 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 𝜌_𝑠𝑐_𝑠𝜕𝑇_𝑠𝜕𝑡=𝑘_𝑠𝜕^𝛼𝑇_𝑠𝜕𝑥^𝛼+𝑞,1<𝛼≤2,0≤𝑥≤𝑆(𝑡). (3) In liquid region 𝜌_𝑙𝑐_𝑙𝜕𝑇_𝑙𝜕𝑡=𝑘_𝑙𝜕^𝛼𝑇_𝑙𝜕𝑥^𝛼+𝑞,1<𝛼≤2,𝑆(𝑡)≤𝑥≤𝑙. (4) The conditions at the moving interface 𝑥=𝑆(𝑡) are given as 𝑘_𝑠𝜕^𝛼−1𝑇_𝑠𝜕𝑥^𝛼−1−𝑘_𝑙𝜕^𝛼−1𝑇_𝑙𝜕𝑥^𝛼−1(𝑡)=𝜌𝐿𝑑𝑆(𝑡),𝑇𝑑𝑡at𝑥=𝑆_𝑠=𝑇_𝑙=𝑇_phat𝑥=𝑆(𝑡), (5)(6) 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 𝐻(𝑇)=^𝑇_𝑇0𝑐𝑑𝑇, (7) where 𝑇₀ is the reference temperature, (3) to (5) become 𝜌𝜕𝐻()𝜕𝜕𝑡=𝐾𝑥,𝑡^𝛼𝑇(𝑥,𝑡)𝜕𝑥^𝛼+𝑞,0≤𝑥≤𝑙. (8) The relation between enthalpy and temperature is given as [17, 18] 𝐻(𝑇)=𝑐_𝑠𝑇−𝑇_𝑚𝑠,for𝑇<𝑇_𝑚𝑠,𝑇−𝑇_𝑚𝑠𝑐_𝑠+𝑐_𝑙2+𝐿𝑇_𝑚𝑙−𝑇_𝑚𝑠for𝑇_𝑚𝑠≤𝑇≤𝑇_𝑚𝑙,1𝐿+2𝑐_𝑠+𝑐_𝑙+𝑐_𝑙𝑇−𝑇_𝑚𝑙,for𝑇>𝑇_𝑚𝑙, (9) and thermal conductivity is given as 𝑘𝐾=^𝑠for𝑇<𝑇_𝑚𝑠,𝑘_𝑠+𝑘_𝑙2for𝑇_𝑚𝑠≤𝑇≤𝑇_𝑚𝑙,𝑘_𝑙,for𝑇>𝑇_𝑚𝑙, (10) where 𝑇_𝑚𝑙 and 𝑇_𝑚𝑠 are liquidus and solidus temperature, respectively.

3.2. Initial and Boundary Conditions

(i) Initially, the system is at 255°C; that is, 𝑇(𝑥,𝑡)=255^∘C,at𝑡=0. (11)

(ii) Convective boundary condition is at 𝑥=0; that is, 𝐾𝜕𝑇(𝑥,𝑡)𝑇()𝜕𝑥=ℎ𝑥,𝑡−𝑇_∞,at𝑥=0, (12) where ℎ is the heat transfer coefficient and 𝑇_∞ is the temperature of the cooling media.

(iii) Adiabatic condition is at 𝑥=𝑙; that is, 𝜕𝑇(𝑥,𝑡)𝜕𝑥=0,at𝑥=𝑙. (13)

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 𝐻_𝑖^𝑛+1=𝐻^𝑛_𝑖+Δ𝑡(Δ𝑥)^𝛼𝜌^𝑛_𝑖𝐾^𝑛_𝑖^𝑖+1_𝑘=0𝑔_𝑘𝑇^𝑛_{𝑖−𝑘+1}+Δ𝑡𝜌^𝑛_𝑖𝑞^𝑛_𝑖,[𝐻]^𝑛+1=[𝐻]^𝑛+Δ𝑡(Δ𝑥)^𝛼[𝐺]^𝑛[𝑇]^𝑛𝑞+Δ𝑡𝜌^𝑛, (14) where 𝐺𝐺=𝐾𝑖,𝑗=0if𝑗≥𝑖+2,𝑖,𝑗^𝑖𝜌_𝑖𝑔₀𝐺=𝐾if𝑗=𝑖+1,𝑖,𝑗^𝑖𝜌_𝑖𝑔₁𝐺=𝐾if𝑗=𝑖,𝑖,𝑗^𝑖𝜌_𝑖𝑔_{𝑖−𝑗+1}if𝑗≤𝑖−1. (15) 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 Δ𝑡/(Δ𝑥)²≤1/𝑑_max2 [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.

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

Table 1 

Thermophysical properties of Sn-5 wt%Pb [21].

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 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 𝑞=500 W/m³ and ℎ=1500 W/m² 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 𝛼=2 at 𝑡=1 second, 𝑡=100 seconds, 𝑡=460 seconds, 𝑡=679 seconds, 𝑡=1000 seconds, 𝑡=2000 seconds, 𝑡=3000 seconds, 𝑡=3779 seconds, 𝑡=3850 seconds, and 𝑡=4265 seconds. Stage 1 ends at 𝑡=1 second when the mushy state starts at 𝑥=0. At 𝑡=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 2 

Temperature profile in region at different time 𝑡 for 𝛼=2.

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 𝑡=1378seconds, 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 3 

Temperature profile in region at different time 𝑡 for 𝛼=1.9.

Figure 4 shows the temperature profile for 𝛼=1.8 at 𝑡=0.3 seconds, 𝑡=34 seconds, 𝑡=100 seconds, 𝑡=274 seconds, 𝑡=1000 seconds, 𝑡=1500 seconds, 𝑡=2084 seconds, 𝑡=2120 seconds, 𝑡=2200 seconds, and 𝑡=2703 seconds. Stage 1 ends at 𝑡=0.3 seconds, when the mushy state starts at 𝑥=0. 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.

Figure 4 

Temperature profile in region at different time 𝑡 for 𝛼=1.8.

Figure 5 shows the temperature profile for 𝛼=1.7 at 𝑡=0.5 seconds, 𝑡=48 seconds, 𝑡=200 seconds, 𝑡=537 seconds, 𝑡=1500 seconds, 𝑡=2500 seconds, 𝑡=3805 seconds, 𝑡=3900 seconds, 𝑡=4000 seconds, and 𝑡=4739 seconds. Stage 1 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 5 

Temperature profile in region at different time 𝑡 for 𝛼=1.7.

Figure 6 shows the temperature profile for 𝛼=1.6 at 𝑡=1second, 𝑡=81 seconds, 𝑡=500 seconds, 𝑡=1000 seconds, 𝑡=1804 seconds, 𝑡=4000 seconds, 𝑡=7000 seconds, 𝑡=8837 seconds, 𝑡=9100 seconds, and 𝑡=10247 seconds. Stage 1 ends at 𝑡=1second, 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. Stage 3 ends at 𝑡=8837 seconds, when the solidification takes place in the entire region. At 𝑡=9100second the system is in stage 4; it ends when it attains its steady state at 𝑡=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 𝛼.

Figure 6 

Temperature profile in region at different time 𝑡 for 𝛼=1.6.

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 7 

Temperature profile for different 𝛼 at 𝑡=1second and 𝑡=36 seconds.

Figure 8 

Temperature profile for different 𝛼 at 𝑡=277 seconds and 𝑡=1378 seconds.

Figure 9 

Temperature profile for different 𝛼 at 𝑡=1812 seconds.

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.

Figure 10 

Liquidus interface for different 𝛼.

Figure 11 

Solidus interface for different 𝛼.

Figures 12 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.

Nomenclature

𝑐∶	Specific heat (J/kg∘C)
𝐻:	Enthalpy of tissue (J/m3)
ℎ:	Heat transfer coefficient (W/(m2°C))
𝑘:	Thermal conductivity of the tissue (W/m∘C)
𝐿:	Latent heat (J/m3)
𝑞:	Volumetric heat generation (W/m3)
𝑆(𝑡):	Position of phase change interface (m)
𝑠:	Solid state
𝑇:	Temperature (°C)
𝑡:	Time (second)
𝑥:	Space coordinate (m)
𝑑max:	Maximum thermal diffusivity (m2/s)
Δ𝑥,Δ𝑡:	Increment in space (m) and time (second)
𝜌:	Density (kg/m3).

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.