International Journal of Numerical Methods for Heat & Fluid Flow

Numerical solution of Stefan problem with variable space grid method based on
mixed finite element/finite difference approach
Milos Ivanovic, Marina Svicevic, Svetislav Savovic,
Article information:
To cite this document:
Milos Ivanovic, Marina Svicevic, Svetislav Savovic, (2017) "Numerical solution of Stefan problem with
variable space grid method based on mixed finite element/finite difference approach", International
Journal of Numerical Methods for Heat & Fluid Flow, Vol. 27 Issue: 12, pp.2682-2695, https://
doi.org/10.1108/HFF-11-2016-0443
Permanent link to this document:
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

https://doi.org/10.1108/HFF-11-2016-0443
Downloaded on: 28 November 2017, At: 09:51 (PT)
References: this document contains references to 34 other documents.
To copy this document: permissions@emeraldinsight.com
The fulltext of this document has been downloaded 10 times since 2017*
Access to this document was granted through an Emerald subscription provided by
Token:Eprints:H4GRK56TZI4444QQZ2IU:
For Authors
If you would like to write for this, or any other Emerald publication, then please use our Emerald
for Authors service information about how to choose which publication to write for and submission
guidelines are available for all. Please visit www.emeraldinsight.com/authors for more information.
About Emerald www.emeraldinsight.com
Emerald is a global publisher linking research and practice to the benefit of society. The company
manages a portfolio of more than 290 journals and over 2,350 books and book series volumes, as
well as providing an extensive range of online products and additional customer resources and
services.
Emerald is both COUNTER 4 and TRANSFER compliant. The organization is a partner of the
Committee on Publication Ethics (COPE) and also works with Portico and the LOCKSS initiative for
digital archive preservation.

*Related content and download information correct at time of download.

 The current issue and full text archive of this journal is available on Emerald Insight at:
www.emeraldinsight.com/0961-5539.htm

HFF
27,12 Numerical solution of Stefan
problem with variable space grid
method based on mixed finite
2682 element/finite difference approach
Received 8 November 2016 Milos Ivanovic and Marina Svicevic
Revised 15 January 2017
17 February 2017
Department of Mathematics and Informatics, University of Kragujevac,
Accepted 3 March 2017 Kragujevac, Serbia, and
Svetislav Savovic
Department of Mathematics and Informatics, University of Kragujevac,
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

Kragujevac, Serbia and Department of Mathematics,

City University of Hong Kong, Kowloon, Hong Kong

Abstract
Purpose – The purpose of this paper is to improve the accuracy and stability of the existing solutions to 1D
Stefan problem with time-dependent Dirichlet boundary conditions. The accuracy improvement should come
with respect to both temperature distribution and moving boundary location.
Design/methodology/approach – The variable space grid method based on mixed finite element/finite
difference approach is applied on 1D Stefan problem with time-dependent Dirichlet boundary conditions
describing melting process. The authors obtain the position of the moving boundary between two phases
using finite differences, whereas finite element method is used to determine temperature distribution. In each
time step, the positions of finite element nodes are updated according to the moving boundary, whereas the
authors map the nodal temperatures with respect to the new mesh using interpolation techniques.
Findings – The authors found that computational results obtained by proposed approach exhibit good
agreement with the exact solution. Moreover, the results for temperature distribution, moving boundary
location and moving boundary speed are more accurate than those obtained by variable space grid method
based on pure finite differences.
Originality/value – The authors’ approach clearly differs from the previous solutions in terms of
methodology. While pure finite difference variable space grid method produces stable solution, the mixed
finite element/finite difference variable space grid scheme is significantly more accurate, especially in case of
high alpha. Slightly modified scheme has a potential to be applied to 2D and 3D Stefan problems.
Keywords Finite element method, Stefan problem, Finite difference variable space grid method,
Front-tracking
Paper type Research paper

1. Introduction
Stefan problems involving heat conduction in conjunction with the change of phase are of
great interest in numerous important areas of science, engineering and industry. Such a
process covers a wide range of applications in which the phase changes from liquid, solid or
International Journal of Numerical
Methods for Heat & Fluid Flow
vapor states. The moving boundary problems occur in many areas such as the metal, glass,
Vol. 27 No. 12, 2017
pp. 2682-2695
plastic and oil industries, space vehicle design, preservation of foodstuffs, chemical and
© Emerald Publishing Limited diffusion processes. The material is assumed to undergo a phase change with a moving
0961-5539
DOI 10.1108/HFF-11-2016-0443 boundary whose position is unknown and has to be determined as a part of the analysis.
 Across the phase boundary, the heat flux is not continuous, and the heat equation is Numerical
replaced by a flux condition that relates the velocity of the phase boundary and the jump of solution of
heat flux across the phase front.
Because moving boundary problems require solving the heat equation in an unknown
Stefan problem
region that has also to be determined as a part of the solution, they are inherently the non-
linear. Hill (1987) reports that due to non-linearity of the moving boundary problems, they
can be solved analytically for only a limited number of special cases. Due to difficulties in
obtaining analytical solutions, numerical techniques, first proposed by Goodman (1958), are 2683
far more common. Numerical techniques are specially known to have difficulties with time-
dependent boundary conditions, and very small time steps are often needed for accurate
solutions. Different numerical methods and their use for different types of problems
are systematized by Minkowycz et al. (2006). Esen and Kutluay (2004), Mennig and Özis ik
(1985), Furzeland (1980), Rizwan-uddin (1999), Savovi
c and Caldwell (2003), Singh et al.
(2011), Vitorino et al. (2010), Irawan et al. (2013) and Reutskiy (2011) report various solutions
of Stefan problems and include linear, exponential and periodical variation of the surface
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

temperature or the flux with time.

There are three main approaches to the solution of the Stefan problem. The first is the
front-tracking method, where the position of the phase boundary is continuously tracked. An
example is the heat balance integral method presented by Goodman (1958), which explicitly
tracks the motion of isotherms (the phase boundary being one of them). Variable grid
methods (variable space grid and variable time step) provide a way to track the phase front
explicitly, as in the work of Marshall (1986). Various solutions for heat conduction problems
with phase change use finite elements (Comini et al., 1974; Lewis et al., 1996, 2004; Nithiarasu
et al., 2016), discontinuous Galerkin method in Nguyen et al. (2010), finite volume method in
the work of Gloth et al. (2003), etc.
The second approach is fixed-domain formulation. An example is the isotherm migration
method, which uses the temperature as the independent variable, reported by Churchill and
Gupta (1977). A more common method is the enthalpy method, which uses an enthalpy
function together with the temperature as dependent variable, reported by Esen and
Kutluay (2004), Voller and Cross (1985), Caldwell and Chan (2000) and Kosec et al. (2011).
Level set method introduces a level set function, capturing the interface at its zero-level set.
The examples of this approach can be found in Vrankar et al. (2010), Reutskiy (2011) and
Papac et al. (2013).
The last approach is the front-fixing formulation. Using a suitable coordinate
transformation, one may immobilize the moving front at the expense of solving a more
complicated problem by a numerical scheme. Examples include finite difference schemes as
Kutluay et al. (1997); Mitchell and Vynnycky (2012) and Vynnycky and Mitchell (2015) and
finite elements in the work of Holmes and Yang (2008).
Comparison of various numerical methods has been carried out by Furzeland (1980),
Caldwell and Kwan (2004) and Javierre et al. (2006). Our approach can be classified as a
front-tracking method.
In this paper, we consider 1D Stefan problem with time-dependent boundary conditions
describing the melting process, identical to that solved by Savovi
c and Caldwell (2009).
Their pure finite difference explicit scheme, although efficient, requires special attention to
the choice of the time step to retain stability. The explicit time-stepping scheme involving
small time-step values often leads to the reduced computational efficiency. Also, pure finite
difference scheme is not as flexible as finite elements at handling natural boundary
conditions. The aim of our paper is to improve the scheme described by Savovi
c and
Caldwell (2009) by using mixed finite element/finite difference scheme including implicit
 HFF time stepping. The position of the moving boundary between two phases is calculated using
27,12 finite differences, whereas we use finite element method for the determination of the
temperature distribution. At each time step, the positions of finite element nodes are updated
according to the moving boundary, whereas we map the nodal temperatures with respect to
the new mesh using interpolation techniques.
The paper is organized as follows. We formulate the problem in Section 2, whereas
2684 Section 3 defines our mixed finite element/finite difference variable space grid (MFE/
FDVSG) scheme. Section 4 gives a brief overview of the FDVSG method used by Savovi
c
and Caldwell (2009). Section 5 presents numerical results and the last section is reserved for
conclusions.

2. Formulation of the problem

Stefan problem describes the phenomenon of the phase change. One-dimensional phase
change problem could be demonstrated by a semi-infinite solid, like a thin block of ice
occupying 0 # x < 1, on solidification temperature. At the fixed boundary of the thin block
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

of ice (x = 0), various types of flux functions could act. In this paper, we used the same
boundary condition as Savovi
c and Caldwell (2009), so the temperature at x = 0 increases
exponentially with time. We also assume that the temperature of the entire solid phase
remains at the melting point. Therefore, the problem is to determine temperature
distribution in the liquid phase at time t0, where x < s(t0) and the location of the free
boundary s(t0). At later time t1 > t0, the moving boundary s(t) moves to the right and
occupies the position s(t1) > s(t0) = s0 (Figure 1). The part of the thin block of ice from
position s(t0) to position s(t1) has been melted over the time interval (t0, t1).
The temperature distribution u(x, t) in the liquid phase region 0 # x # s(t) is given by the
heat equation:

@u @2u
¼a 2 (1)
@t @x

under the following boundary conditions:

uðx; tÞ ¼ eat ; x ¼ 0; t>0

(2)
uðx; tÞ ¼ 1; x ¼ sðtÞ; t>0

where a is a physical parameter combining density, specific heat and thermal conductivity.
The location of the moving boundary is given by equation known as Stefan condition:

1 ds @u
¼ ; x ¼ sðtÞ; t>0 (3)
a dt @x

In the general case, initial condition is given by:

Figure 1.
Ice melting
 sð0Þ ¼ 0 (4) Numerical
The exact solution for this specific problem is known:
solution of
Stefan problem
uðx; tÞ ¼ eatx
(5)
sðtÞ ¼ at

As stated above, to determine s(t), we used FDVSG scheme, and for obtaining u(x, t), we used 2685
the finite element method. It should be noted that we consider a generic solid melting to
liquid (depending on the value of a), not just the ice-water case. The following section
describes our approach in more detail.

3. Mixed finite element/finite difference variable space grid method

The finite element method, as well as the FDVSG proposed by Savovi
c and Caldwell (2009),
can be used to estimate temperature distribution in the fixed number of nodes in the liquid
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

phase. Our novel procedure consists of solving equation (1) using finite element method
while keeping the FDVSG given by Savovi
c and Caldwell (2009) as a solution to equation (3).
We designate our combined approach proposed in this paper as MFE/FDVSG method.

3.1 Finite element scheme for the temperature distribution

In this subsection, the derivation of the continuous and discrete equations for a single
element is briefly presented. Equation (1) applies Galerkin method in the way that we
multiply both sides of the equation by weighting function v, integrate throughout the
volume and apply Gauss theorem, giving:
ð ð þ
@u @v @u du
a dV ¼  vdV þ a vdS (6)
V @x @x V @t S dx

This is the weak formulation of equation (1). The problem is to determine the function u
such that equation (6) holds for every v. The value of the function can
Pbe written as a sum of
the nodal values and the interpolation functions at given point u ¼ f j uj . Furthermore, in
place of weighting function introduced in equation (6), the interpolation function f i can be
adopted, where i ¼ 1; K, K being the number of finite element nodes. We get:

@uj
Ci;j þ Ki;j uj ¼ Fi (7)
@t

where:
ð
Ci;j ¼ f i f j dV (8)
V

ð
@f i @f j
Ki;j ¼ a dV (9)
V @x @x

þ
du
Fi ¼ a f dS (10)
S dx i
 HFF In the case that the given element is located in the domain interior, due to the orientation of the
27,12 surface vector, the expression Fi is annulled by the corresponding integral of the adjacent
elements. For the boundary elements, according to the conditions, the value of the given
integral is also zero. As the model implies fixed values on both ends given by equation (2), we
pose these Dirichlet boundary conditions by applying standard penalty method.

2686 3.2 Boundary update and remeshing

Besides solving system equation (7) in each time step, it is necessary to update the position
of the moving liquid phase boundary. Also, while boundary is moving, the mesh has to be
regenerated to keep constant number of elements with equally distributed nodes throughout
the simulation. The temperature gradient at the moving boundary x = s(t) = ND x is given
by following finite difference higher three point backward scheme described by Savovi
c and
Caldwell (2009):

@u  3uN  4uN1 þ uN2
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

 ¼ þ OðDx2 Þ (11)
@x x¼sðtÞ 2Dx

According to equation (11), the Stefan condition equation (3) at x = s(t) can be written as:

Dta  
smþ1 ¼ sm  3uðsm Þ  4u ðsm  hm Þ þ u ðsm  2hm Þ ;
2hm (12)
m ¼ 0; 1; 2; . . .

where m denotes time-step count, D t time-step value and hm ð DxÞ ¼ sNðtÞ nominal space
step. The quantities designated as uðsm Þ; u ðsm  hm Þ and u ðsm  2hm Þ are temperatures at
sm, sm–hm and sm–2hm, in time step m, respectively. It should be noted that these
temperatures usually do not correspond to the nodal values (hence, bar notation). The reason
for that is the addition of an extra node in the new mesh (in time step m þ 1), exactly at the
position of the previous liquid-ice boundary (at sm). We expect the inclusion of this
additional node to improve accuracy. However, an extra node requires interpolation for
determining temperatures u ðsm  hm Þ and u ðsm  2hm Þ. The interpolation is also required
for the transition from the old mesh with nominal step hm to the new mesh with nominal
step hmþ1 ¼ smþ1N . For both purposes we used Akima scheme proposed by Franke (1982),
which provided a slight improvement over classic cubic spline interpolation defined by
McKinley and Levine (1998).
The procedure of the transition between time steps is shown in Figure 2. It should be noted
that mesh consists of N nodes only in the first time step; in all subsequent steps, there are
actually N þ 1 nodes. Of course, if the projection of the node to the previous time step falls
within ice phase, we take temperature as 1, as shown at the far right side of Figure 2. Using this
approach, we have all necessary building blocks to form equation (7) for the time step m þ 1.

4. Finite difference variable space grid method

In this section, we briefly describe the pure finite difference approach taken by Savovi
c and
Caldwell (2009). Taking into account continuous change in the nodal positions due to the
boundary movement, the following expression applies at the i-th grid point:
   
@u  @u  @x  @u 
¼ þ (13)
@t i @x t @t i @t x
 whereas the node xi moves according to: Numerical
dx xi ds
solution of
¼  (14) Stefan problem
dt sðtÞ dt

in which suffixes t, i, and x have to be kept constant during differentiation process and omitted
for clarity in further text. Upon substituting equations (1) and (14) into equation (13), we obtain:
2687
@u xi ds @u @ u 2
¼ þa 2; 0 # x # sðtÞ; t > 0 (15)
@t s dt @x @x

Having temperature gradient at the boundary determined by equation (11), forward

difference approximation for time derivative and central difference approximation for the
space derivative, discretization of equation (15) can be expressed as:
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

kxi;m s_ m ka  
ui;mþ1 ¼ ui;m þ ðuiþ1;m  ui1;m Þ þ 2 uiþ1;m  2ui;m þ ui1;m (16)
2hm sm hm

where ui;m  uðxi;m Þ; s_ m ¼ ðsmþ1  sm Þ=Dt; xi;m ¼ ihm ; tm ¼ t0 þ mk, hm is the space grid
size D x at mth time step, k : Dt is the time step, and t0 is the time at which the numerical
process is initialized.
We obtain the position of the moving boundary using identical formula (12), as in the
case of MFE/FDVSG method, whereas updated grid size becomes hmþ1 = smþ1/N. Further
details on FDVSG approach can be found in Savovi
c and Caldwell (2009).

5. Results and discussion

In this section, we present the results obtained using our MFE/FDVSG scheme, applied to
one dimensional Stefan problem of ice melting. Our tests aimed at obtaining the following:
 accuracy of the temperature distribution;
 time step value sensitivity;
 accuracy of the location and speed of the moving boundary; and
 sensitivity regarding mesh density and time step, only for MFE/FDVSG method.

Figure 2.
Procedure of
updating solution
from the previous
time step including
remeshing
 HFF First of all, we compare our computational results for heat distribution with the exact
27,12 solution for a = 2 and a = 10 and also with those obtained by Savovi
c and Caldwell (2009),
who used pure FDVSG. The parameter values are identical as in Savovi
c and Caldwell
(2009). Nominal mesh density is fixed to ten nodes for both MFE/FDVSG and FDVSG
schemes. The initial time is t0 = 0.01, which, according to equation (5), corresponds to the
initial position of the moving boundary. In the case of a = 2, the boundary starts at s(t0) =
2688 0.02, whereas in the case of a = 10, the boundary starts at s(t0) = 0.1. The position of the
boundary between ice and water, together with the temperature distribution, is observed in

t x/s MFE/FDVSG error FDVSG error

0.1 0 0 0
0.1 0.00002496 0.00003222
0.2 0.00003333 0.00003587
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

0.3 0.00003602 0.00001741

0.4 0.00004023 0.00001816
0.5 0.00005058 0.00006696
0.6 0.00007064 0.00012614
0.7 0.00010322 0.00019363
0.8 0.00015076 0.00026803
0.9 0.00023918 0.00034846
1 0.00043449 0.00043449
0.3 0 0 0
0.1 0.00472279 0.02209676
0.2 0.00180897 0.02962343
0.3 0.00279347 0.02938231
0.4 0.00655247 0.02537378
0.5 0.00869978 0.01985605
0.6 0.00910542 0.01402515
0.7 0.00792414 0.00844575
0.8 0.00539308 0.00332003
0.9 0.00119794 0.00134782
1 0.00565342 0.00565343
0.5 0 0 0
0.1 0.24347081 0.67711050
0.2 0.15592026 0.75866250
Table I. 0.3 0.02961707 0.63241963
Temperature 0.4 0.05352178 0.46327442
0.5 0.09062042 0.31238079
distribution u(x, t) in
0.6 0.09428687 0.19557417
nodal positions 0.7 0.07782808 0.11096057
obtained by FDVSG 0.8 0.05096659 0.05123940
and the MFE/ 0.9 0.01751825 0.00871911
FDVSG, a = 10 1 0.02299748 0.02299740

a=2 a = 10
Table II.
t FDVSG error MFE/FDVSG error FDVSG error MFE/FDVSG error
The accuracy of
FDVSG and MFE/ 0.1 5.8152E-06 5.2901E-06 2.0061E-04 1.6277E-04
FDVSG methods at 0.3 3.5864E-05 2.7875E-05 1.8041E-02 5.7222E-03
various times (t) 0.5 1.3531E-04 1.0180E-04 4.0452E-01 1.0177E-01
 the interval t0 = 0.01 till tmax = 0.5. We also took identical time-step values used by Savovi
c Numerical
and Caldwell (2009), D t = 1.0E06 for a = 2 and D t = 2.0E–06 for a = 10. This choice of the solution of
time-step value guarantee stability of differential scheme used by the pure FDVSG method.
In Table I, we compare temperature the distributions for various time points in the case of
Stefan problem
a = 10.

2689
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

(a)

Figure 3.
The accuracy of
MFE/FDVSG and
FDVSG methods for
(a) a = 2 and for
(b) a = 10
(b)
 HFF Because Table I is insufficient for the exact comparison of the two numerical methods, we
27,12 considered the root mean square error defined by:
vffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi
u N 
u1 X 2
Error ¼ t umethod
i  uexact
i (17)
N i¼1
2690
where N is the total number of observed points. Equation (17) was taken as the error
function, representing accuracy evaluation of the method. As the error value decreases, a
method gives better heat distribution over a given time interval.
Table II and Figure 3 represent accuracy of the two methods for three time instants. It
can be noted that for earlier times, the FDVSG method provides solution on par with the
MFE/FDVSG, but with increasing time, the MFE/FDVSG gives better stability and better
match the exact solution, especially for a = 10.
The accuracy of the temperature distribution for the MFE/FDVSG method for a = 2 is
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

significantly better than the accuracy of the results achieved for a = 10, as shown in
Figure 3. Because a values for almost all technically applicable materials are less than 5, it
can be assumed that both methods are sufficiently accurate for most practical applications.
Therefore, in all further tests, the value of parameter a will be fixed to the value of 2.
The second accuracy test was the time step value sensitivity. Due to the stability
requirements of the FDVSG method, we used value Dt = 1.0E–06 for all FDVSG runs. In

Table III.
Error of MFE/ Dt Error MFE/FDVSG Error FDVSG
FDVSG and FDVSG 1.0E-03 1.6189E-04 –
method for various 1.0E-04 1.2982E-04 –
time step values at 1.0E-05 1.3180E-04 –
t = 0.5 s 1.0E-06 1.3228E-04 1.3531E-04

t MFE/FDVSG FDVSG Error MFE/FDVSG Error FDVSG Exact value

s(t)
0.02 0.0399652 0.0399999 3.480E-05 5.160E-08 0.04
0.05 0.0999859 0.0999991 1.410E-05 8.141E-07 0.1
0.1 0.1999999 0.1999936 1.000E-07 6.377E-06 0.2
0.2 0.3999930 0.3999509 7.000E-06 4.907E-05 0.4
0.3 0.5999228 0.5998399 7.720E-05 1.600E-04 0.6
0.4 0.7999311 0.7996321 6.890E-05 3.679E-04 0.8
Table IV. 0.5 0.9998902 0.9993009 1.098E-04 6.990E-04 1
Position of moving
boundary and its ds/dt
velocity obtained 0.02 2.0017569 1.9999920 1.756E-03 7.967E-06 2.0
0.05 2.0002564 1.99995148 2.564E-04 4.852E-05
using MFE/FDVSG
0.1 1.9999762 1.9998113 2.378E-05 1.887E-04
and FDVSG methods 0.2 1.9996765 1.9992818 3.234E-04 7.182E-04
compared to the 0.3 1.9993043 1.9984515 6.956E-04 1.548E-03
exact solution for 0.4 1.9988485 1.9973471 1.151E-03 2.653E-03
a=2 0.5 1.9983220 1.9959886 1.677E-03 4.011E-03
 contrast, as can be seen from Table III, the MFE/FDVSG method remains stable, even with a Numerical
significantly higher time step value. solution of
The third test was about determining the location and speed of the moving boundary
through time, based on equation (3). Although the same equation is used for both numerical
Stefan problem
schemes, we expect better results than those reported by Savovi
c and Caldwell (2009), due to
a more accurate temperature distribution achieved using finite element method. The results
2691
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

(a)

Figure 4.
(a) The absolute error
of the boundary
position and (b) the
absolute error of the
boundary speed
(b)
 HFF that confirm this assumption are presented in Table IV and Figure 4. The MFE/FDVSG is
27,12 clearly more accurate than the pure FDVSG, which is especially visible considering higher
time points, with accumulated errors.
Finally, we performed the fourth test, which considers variable time step value and
variable finite element mesh density, only for the MFE/FDVSG method. The tests were
carried out with four time step values taken as a geometric sequence (D t) and the nominal
2692 number of finite elements 10, 20, 50 and 100 (N), using equation (17) as the error descriptor.
The obtained results are presented in Table V and Figure 5. It should be noted that in
Figure 5, we single out and show only italic items from Table V, where the stability for
certain N is reached (no significant improvement with further time step decrease). The

Dt 1.0E-03
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

N 10 20 50 100
Error 1.6189E-04 1.4839E-04 1.7269E-04 1.9599E-04
Dt 1.0E-04
N 10 20 50 100
Error 1.2982E-04 2.3488E-05 1.4841E-05 1.5383E-05
Dt 1.0E-05
Table V. N 10 20 50 100
The accuracy of the Error 1.3180E-04 2.0955E-05 2.7023E-06 1.2072E-06
MFE/FDVSG method Dt 1.0E-06
for various time step N 10 20 50 100
Error 1.3228E-04 2.1220E-05 2.7437E-06 9.4844E-07
values and nominal
number of elements Note: The italic cells designate points where there’s no significant improvement with further time step
(extra node present) decrease for certain N

Figure 5.
The accuracy of the
MFE/FDVSG method
for time/space step
combinations
denoted as Dt
N (italic Note: For the sake of comparison, we also include the case when
cells in Table V)
no extra node is introduced
 Figure 5 also considers the case when no extra node is introduced. It is clear that the addition Numerical
of the extra node improves the accuracy for over two orders of the magnitude. solution of
Stefan problem
6. Conclusions
We report the implementation of the variable space grid method based on mixed finite
element/finite difference approach. The MFE/FDVSG method itself represents a variation of
the finite difference FDVSG scheme proposed by Savovi
c and Caldwell (2009). The 2693
difference between FDVSG and MFE/FDVSG scheme consists in using the finite element
method with remeshing, instead of pure finite difference approach, for the purpose of
solving the heat equation. The method is applied to 1D Stefan problem with the time-
dependent Dirichlet boundary conditions describing melting process. The novel MFE/
FDVSG exhibits very good agreement with exact solution and a clear benefit in terms of
accuracy and stability over pure finite difference FDVSG method described in Savovi
c and
Caldwell (2009). The benefit is especially visible in the case of a = 10. The proposed method
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

is also significantly less sensible to the time-step value choice. All presented results
encourage us to apply proposed approach to the problems and cases where analytical
solutions are not available. Moreover, slightly modified scheme has a potential to be applied
for 2D and even 3D Stefan problems, with the inclusion of the automatic mesh generation
procedure.

References
Caldwell, J. and Chan, C.-C. (2000), “Spherical solidification by the enthalpy method and the heat
balance integral method”, Applied Mathematical Modelling, Vol. 24 No. 1, pp. 45-53.
Caldwell, J. and Kwan, Y.Y. (2004), “Numerical methods for one-dimensional Stefan problems”,
Communications in Numerical Methods in Engineering, Vol. 20 No. 7, pp. 535-545.
Churchill, S.W. and Gupta, J.P. (1977), “Approximations for conduction with freezing or melting”,
International Journal of Heat and Mass Transfer, Vol. 20 No. 11, pp. 1251-1253.
Comini, G., Del Guidice, S., Lewis, R.W. and Zienkiewicz, O.C. (1974), “Finite element solution of non-
linear heat conduction problems with special reference to phase change”, International Journal
for Numerical Methods in Engineering, Vol. 8 No. 3, pp. 613-624.
Esen, A. and Kutluay, S. (2004), “A numerical solution of the Stefan problem with a Neumann-type
boundary condition by enthalpy method”, Applied Mathematics and Computation, Vol. 148
No. 2, pp. 321-329.
Franke, R. (1982), “Scattered data interpolation: tests of some methods”, Mathematics of Computation,
Vol. 38 No. 157, pp. 181-200.
Furzeland, R.M. (1980), “A comparative study of numerical methods for moving boundary problems”,
IMA Journal of Applied Mathematics, Vol. 26 No. 4, pp. 411-429.
Gloth, O., Hänel, D., Tran, L. and Vilsmeier, R. (2003), “A front tracking method on unstructured grids”,
Computers & Fluids, Vol. 32 No. 4, pp. 547-570.
Goodman, T.R. (1958), “The heat-balance integral and its application to problems involving a
change of phase”, Transactions of the American Society of Mechanical Engineers, Vol. 80
No. 2, pp. 335-342.
Hill, J. (1987), One-Dimensional Stefan Problems: An Introduction, Longman Scientific and Technical,
Harlow, Essex.
Holmes, A.D. and Yang, H. (2008), “A front-fixing finite element method for the valuation of american
options”, SIAM Journal on Scientific Computing, Vol. 30 No. 4, pp. 2158-2180.
 HFF Irawan, B., Wardana, I., Wahyudi, S. and Argo, B.D. (2013), “Numerical methods for predicting melting
through the gap applied to the one-dimensional Stefan problem using Taylor series”, Applied
27,12 Mathematical Sciences, Vol. 7 No. 91, pp. 4509-4519.
Javierre, E., Vuik, C., Vermolen, F.J. and Van der Zwaag, S. (2006), “A comparison of numerical models
for one-dimensional stefan problems”, Journal of Computational and Applied Mathematics,
Vol. 192 No. 2, pp. 445-459.
2694 Kosec, G., Založnik, M., Šarler, B. and Combeau, H. (2011), “A meshless approach towards solution of
macrosegregation phenomena”, Computers, Materials, & Continua, Vol. 22 No. 2, pp. 169-195.
Kutluay, S., Bahadir, A.R. and Özdes , A. (1997), “The numerical solution of one-phase classical stefan
problem”, Journal of Computational and Applied Mathematics, Vol. 81 No. 1, pp. 135-144.
Lewis, R.W., Nithiarasu, P. and Seetharamu, K.N. (2004), Fundamentals of the Finite Element Method
for Heat and Fluid Flow, John Wiley & Sons.
Lewis, R.W., Morgan, K., Thomas, H.R. and Seetharamu, K. (1996), The Finite Element Method in Heat
Transfer Analysis, John Wiley & Sons.
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

McKinley, S. and Levine, M. (1998), “Cubic spline interpolation”, College of the Redwoods, Vol. 45,
pp. 1049-1060.
Marshall, G. (1986), “A front trackingmmethod for one-dimensional moving boundary problems”,
SIAM Journal on Scientific and Statistical Computing, Vol. 7 No. 1, pp. 252-263.
Mennig, J. and Özis ik, M.N. (1985), “Coupled integral equation approach for solving melting or
solidification”, International Journal of Heat and Mass Transfer, Vol. 28 No. 8, pp. 1481-1485.
Minkowycz, W.J., Sparrow, E.M. and Murthy, J. (2006), Handbook of Numerical Heat Transfer, number
(Sirsi) i9780471348788, Wiley Online Library.
Mitchell, S.L. and Vynnycky, M. (2012), “An accurate finite-difference method for ablation-type
stefan problems”, Journal of Computational and Applied Mathematics, Vol. 236 No. 17,
pp. 4181-4192.
Nguyen, V.T., Peraire, J., Khoo, B.C. and Persson, P.O. (2010), “A discontinuous galerkin front
tracking method for two-phase flows with surface tension”, Computers & Fluids, Vol. 39
No. 1, pp. 1-14.
Nithiarasu, P., Lewis, R.W. and Seetharamu, K.N. (2016), Fundamentals of the Finite Element Method
for Heat and Mass Transfer, John Wiley & Sons.
Papac, J., Helgadottir, A., Ratsch, C. and Gibou, F. (2013), “A level set approach for diffusion and stefan-
type problems with robin boundary conditions on quadtree/octree adaptive cartesian grids”,
Journal of Computational Physics, Vol. 233, pp. 241-261.
Reutskiy, S.Y. (2011), “A meshless method for one-dimensional Stefan problems”, Applied Mathematics
and Computation, Vol. 217 No. 23, pp. 9689-9701.
Rizwan-uddin, (1999), “One-dimensional phase change with periodic boundary conditions”, Numerical
Heat Transfer: Part A: Applications, Vol. 35 No. 4, pp. 361-372.
Savovi
c, S. and Caldwell, J. (2003), “Finite difference solution of one-dimensional Stefan problem with
periodic boundary conditions”, International Journal of Heat and Mass Transfer, Vol. 46 No. 15,
pp. 2911-2916.
Savovi
c, S. and Caldwell, J. (2009), “Numerical solution of Stefan problem with time-dependent boundary
conditions by variable space grid method”, Thermal Science, Vol. 13 No. 4, pp. 165-174.
Singh, J., Gupta, K.P. and Rai, N.K. (2011), “Variational iteration method to solve moving boundary
problem with temperature dependent physical properties”, Thermal Science, Vol. 15 No. S2,
pp. 229-239.
Vitorino, N., Abrantes, J.C.C. and Frade, J.R. (2010), “Numerical solutions for mixed controlled
solidification of phase change materials”, International Journal of Heat and Mass Transfer,
Vol. 53 No. 23, pp. 5335-5342.
 Voller, V.R. and Cross, M. (1985), “Applications of control volume enthalpy methods in the solution of Numerical
Stefan problems”, Computational Techniques in Heat Transfer, Pineridge Press Ltd., Mambles,
Swansea, pp. 245-275. solution of
Vrankar, L., Kansa, E.J., Ling, L., Runovc, F. and Turk, G. (2010), “Moving-boundary problems solved Stefan problem
by adaptive radial basis functions”, Computers & Fluids, Vol. 39 No. 9, pp. 1480-1490.
Vynnycky, M. and Mitchell, S.L. (2015), “On the numerical solution of a Stefan problem with finite
extinction time”, Journal of Computational and Applied Mathematics, Vol. 276, pp. 98-109.
2695
Corresponding author
Milos Ivanovic can be contacted at: mivanovic@kg.ac.rs
Downloaded by 178.149.6.158 At 09:51 28 November 2017 (PT)

For instructions on how to order reprints of this article, please visit our website:
www.emeraldgrouppublishing.com/licensing/reprints.htm
Or contact us for further details: permissions@emeraldinsight.com