This paper presents a fractional mathematical model of a one-dimensional phase-change problem (Stefan problem) with a variable latent-heat (a power function of position). This model includes space–time fractional derivatives in the Caputo sense and time-dependent surface-heat flux. An approximate solution of this model is obtained by using the optimal homotopy asymptotic method to find the solutions of temperature distribution in the domain 0 ≤x≤s(t) and interface’s tracking or location. The results thus obtained are compared with existing exact solutions for the case of the integer order derivative at some particular values of the governing parameters. The dependency of movement of the interface on certain parameters is also studied.

This paper presents a fractional mathematical model of a one-dimensional phase-change problem (Stefan problem) with a position’s latent-heat power function. This model includes space–time fractional derivatives in the Caputo sense and time-dependent surface-heat flux. An approximate solution of this model is obtained using the optimal homotopy asymptotic method to find approximate solutions of temperature distribution in the domain and using the interface’s tracking or location. The results thus obtained are compared with existing exact solutions for the case of the integer order derivative at some particular values of the governing parameters. A study also is performed of the interface movement’s dependence on certain parameters.

In recent years, many researchers have used fractional derivatives in various mathematical models due to their applicability in different fields of science and engineering [1-4]. It is well known that a fractional derivative is a good tool for taking into account the memory mechanism, particularly in some diffusive processes [5]. Stefan problems (moving boundary problems) with fractional derivatives [6-10] are typical problems from a mathematics point of view because of their nonlinear nature and the presence of a moving interface. Some exact solutions to Stefan problems can be seen in [8], [11], and [12]. Exact solutions to such problems are limited. Therefore, several approximate analytical methods [13-17] have been used to solve the Stefan problems with fractional derivatives. The approximate analytical method used in this literature is the optimal homotopy asymptotic method (OHAM).

The OHAM was developed by Marinca et al. [18], and it has been applied to solve a wide range of nonlinear differential equations [19-23]. Ghoreishi et. al. [24] presented the comparison between the homotopy analysis method and the OHAM for nonlinear age-structured population models. In 2013, Dinarvand and Hosseini [25] also used the OHAM to investigate the temperature distribution equation in a convective straight fin with temperature-dependent thermal conductivity and the convective–radiative cooling of a lumped system with variable specific heat.

This paper presents a mathematical model for a Stefan problem [12] with a space–time fractional derivative. In this model, the OHAM is used to find the expression of the temperature distribution in a given domain and location of a moving interface with the help of the Taylor series [13]. The obtained results are compared with the existing exact solution for a standard case and are in good agreement. An approachable analysis for a fractional case also is discussed.

2. Mathematical formulation

In this section, a mathematical model of a one-dimensional Stefan problem with a position’s latent-heat power function [12] is considered. For this problem, we present a fractional model that involves space–time fractional derivatives, as given in [11]. The governing equations are as follow:

(1)

(2)

(3)

, (4)

where is the temperature distribution, is the moving interface, is thermal conductivity, is the thermal diffusion coefficient, b is a constant ( for melting, for freezing), is the variable latent heat per unit volume, and is an non-negative integer. The operators and are the Caputo fractional derivatives [11,13], which are defined as

(5) , (6)

where is the Gamma function.

In this paper, the following properties of fractional derivatives [13-14] are used:

(a) (7)

where is a constant.

(b) , (8)

where and is the Caputo fractional derivative of .

3. Solution of the problem

First, Eqs. (1)–(4) are written in operator form as follows:

(9) , (10)

where is a linear operator, is a nonlinear operator, and B is a boundary operator.

According to the OHAM [16, 21], we construct an optimal , which satisfies

(11)

(12)

where is an embedding parameter, is an unknown function, and is a nonzero auxiliary function for and . Obviously, if ,

, (13)

and when , then

. (14)

Clearly, as p increases from 0 to 1, the unknown function varies from to the solution .

Now, we choose the auxiliary function in the following form:

, (15)

where are constants to be determined later.

The solution to Eq. (11) is considered in the following series form:

(16)

and

, (17)

where and .

Now, we expand the nonlinear term into the following series form (as given in [24]):

(18)

where .

Now, by substituting Eqs. (16) and (18) into Eq. (11) and equating the coefficients of like powers of , the following problems are obtained:

, (19)

(20)

(21)

and the general equation for is given as

(22)

where .

Substituting Eqs. (16) and (17) into the boundary conditions of (6) and (7) provides the following:

(23)

and

, (24)

where .

The comparison of various powers of can be shown by expending in Taylor series form [13, 14] at a point, , as follows:

(25)

where and .

Eqs. (24) and (25) provide the following:

. (26)

The interface condition (4) becomes

(27)

By considering Eq. (19) and comparing the coefficients of the power of from Eqs. (23), (26), and (27), the following system can be obtained:

(28)

Taking Eq. (20) and comparing the coefficients of power for from Eqs. (23), (26), and (27) provides the following:

(29)

Similarly, other systems can be found by comparing various powers of .

The solutions of the zeroth-order problem (28) are calculated as the following:

, (30)

and

, (31)

Where , .

Substituting and into the first-order problem (29) and using the above process obtains the following expressions of :

(32)

where ,

,

and .

The expression of can be calculated as:

, (33)

where ,

and .

The approximate solution of the temperature distribution can be determined as

, (34)

and an approximate solution of is given as

. (35)

In order to get the constants involved in Eq. (34) for the expression of , the least square method is used [24]. For this purpose, residual is defined as:

(36)

If , then will be the exact solution. Generally, the OHAM gives an approximate solution. Therefore, in such a case, , but the function can be minimized as

, (37)

where R is the residual. The constants can be obtained optimally from the following conditions:

(38)

4. Numerical results and discussion

In this section, numerical results for interface position are obtained with the help of Wolfram Research (8.0.0) software by considering only , and the results are presented through tables and figures.

Table 1. Comparison between exact and approximate solution of at n = 0.

b

t

Exact value of S(t)

Approximate value of S(t) by OHAM

Error (%)

0.1

0

5

10

15

20

25

0.0

0.4428497

0.6262841

0.7668507

0.8856994

0.9902421

0.0

0.4449844

0.6293030

0.7707735

0.8899689

0.9950155

0.0

0.4820

0.4820

0.5112

0.4820

0.4820

0.25

0

1

2

3

4

5

0.0

0.4728215

0.6686706

0.8189509

0.9456430

1.0030059

0.0

0.4847701

0.6854703

0.8395262

0.9694014

1.0282054

0.0

2.5271

2.5124

2.5124

2.5124

2.5124

0.5

0

0.25

0.50

0.75

1.00

1.30

0.0

0.4193648

0.5930714

0.7263611

0.8387296

0.9562989

0.0

0.4439988

0.6279091

0.7690284

0.8879975

1.0124729

0.0

5.8741

5.8741

5.8741

5.8741

5.8741

Table 2. Comparison between exact and approximate solution of at n = 1.

b

t

Exact value of S(t)

Approximate value of S(t) by OHAM

Error (%)

0.1

0

1

2

3

4

5

0.0

0.4275253

0.6046120

0.7404955

0.8550505

0.9559756

0.0

0.4332152

0.6126588

0.7503507

0.8664304

0.9686986

0.0

1.3309

1.3309

1.3309

1.3309

1.3308

0.25

0

0.5

1.0

1.5

2.0

2.5

0.0

0.4527556

0.6402931

0.7841957

0.9055112

1.0123923

0.0

0.4719994

0.6675079

0.8175269

0.9439988

1.0554227

0.0

4.2504

4.2504

4.2504

4.2504

4.2504

0.5

0.0

0.25

0.5

0.75

1.0

1.25

0.0

0.4222513

0.5971534

0.7313607

0.8445026

0.9441826

0.0

0.4536318

0.6415322

0.7857133

0.9072635

1.0143515

0.0

7.4317

7.4317

7.4317

7.4317

7.4317

Tables 1–2 represent comparisons between the exact and approximate values of the phase front ’s positions at particular times for (standard motion). The tables clearly show that the approximate results are sufficiently accurate and in agreement with the existing exact solution [12] for standard motion.

Figs. 3 and 4 also depict the dependence of phase front ’s path on the thermal diffusion coefficient for n = 1at and , respectively. Figs. 1–4 portray that the interface’s movement increases with an increase in the value of the thermal diffusion coefficient for fractional cases (nonclassical or non-Fickian), which is similar to the case of standard motion [12].

b =0.15

b=0.1

b=0.05

Fig. 5 Plot of vs. at and n = 0.

b =0.15

b=0.10

b=0.05

Fig. 6 Plot of vs. at , and n = 0.

Figs. 5–6 show a variation in ’s path for a different value of b for a non-classical or non-Fickian case. From these figures, it is clear that the phase front’s movement increases with an increase in the value of the constant b; that is, the melting (or freezing) process becomes fast as the value of the constant b increases.

5. Conclusion

In this work, we considered a mathematical model that contains space–time fractional derivatives and time-dependent surface-heat flux. An approximate solution of the model was obtained by the OHAM. It was observed that the interface movement increases with an increase in the value of the thermal diffusion coefficient v as well as the constant b for a nonclassical or non-Fickian case. Moreover, it can be seen that the proposed technique is sufficiently accurate and efficient for solving Stefan problems. It also was observed that it is convenient for controlling and adjusting the convergence of the series solution through the control parameters in the OHAM.

Acknowledgements

The authors express their sincere thanks to the anonymous referees for their valuable suggestions for the improvement of the paper.