# Two Simple Numerical Methods for the Free Boundary in One-Phase Stefan Problem

Authors: Seung Hyun Kim
Journal: Journal of Applied Mathematics
 (2014)

---

## Abstract

We present two simple numerical methods to find the free boundary in one-phase Stefan problem. The work is motivated by the necessity for better understanding of the solution surface (temperatures) near the free boundary. We formulate a log-transform function with the unfixed and fixed free boundary that has Lipschitz character near free boundary. We solve the quadratic equation in order to locate the free boundary in a time-recursive way. We also present several numerical results which illustrate a comparison to other methods.


---

## Body

## 1. Introduction

The free boundary in a Stefan problem can give prediction of the boundary between two opposing materials. For a famous example, the boundary of an ice in water can be found according to time. In this case, water above zero degrees Celsius and ice at zero degrees Celsius are the two materials. Stefan [1] first gave conditions necessary to find the free boundary. Stefan problems can be used for the computation of diffusion [2], preservation of human blood [3], solidification of metals [4], and estimation of the thermal state in the mantle [5]. Many researchers have studied to figure out the effective method in order to approximate the free boundary.Methods can be largely separated into three groups: analytic methods, mixed methods of numerical and analytic methods, and another method which is numerical methods. Analytic methods are developed by researchers, such as exact solution by using a fractional differential equation of Liu and Xu [6] and an analytical solution with variable latent heat of Voller et al. [7]. Mixed methods are developed by researchers, such as the approximate analytical and a numerical solution method by using a linear perturbation method of Yigit [8] and the analytical and numerical solution method with moving boundary of Lorenzo-Trueba and Voller [9]. Numerical methods are developed by researchers, such as the numerical method with automatic differentiation of Asaithambi [10] and the finite element approximation method of Nedoma [11], Słota [12], and Hinze and Ziegenbalg [5].Among many methods, a boundary immobilization method (BIM), suggested by Landau [13], fixes the free boundary by using a transform method. The benefit of this function is its simple domain. Whatever the shape of the free boundary is, this method changes the original domain to a rectangular domain. After the BIM method, Kutluay et al. [14] suggested another method called variable space grid (VSG) method.The main contribution of this paper is the development of two simple numerical methods to find the free boundary in a time-recursive way. Our result is motivated by the necessity for better understanding of the solution surface near free boundary. We adopt the boundary immobilization [13] to change the unknown free boundary to a known and fixed boundary. On the other hand, we adopt the constant number of space intervals between two boundaries [14].We exploit a log-transform function with the unfixed and fixed free boundary that has Lipschitz character which reduces the accumulating error of the solution surface near the free boundary.We apply extrapolation near the free boundary. Thus, we can determine the free boundary by solving a quadratic equation in a time-recursive way. Our method also provides fast and accurate results for calculating the free boundary and determining solutions.The structure of the paper is as follows. Section2 presents the model formulation. The log-transform function with the unfixed and fixed free boundary to calculate the free boundary is presented in Section 3. Numerical results and comparative studies are presented in Section 4. Section 5 summarizes the paper.


## 2. Problem Formulation

In this section, we present a mathematical formula for determining solution of one-phase Stefan problem.The temperature of liquid is denoted byu(x,t), where x is the liquid/solid space for x∈[0,∞) and t is the time to expiration for t∈[0,T].As seen in the previous article by Stefan [1], the temperature of liquid is considered the solution to a free boundary problem with a parabolic PDE. We suppose that the free boundary S(t) is continuously increasing with S(0)=1. The region where it is free boundary to hold, generally called the liquid region, is defined as ΩL=[0,S(t)]×[0,T], and the region where it is free boundary, generally called the solid region, is defined as ΩS=(S(t),∞)×[0,T]. Then, u(·,·) and S(·) uniquely solve
(1)∂u∂t=∂2u∂x2inΩL,
being subject to the boundary conditions
(2)-ux(0,t)=f(t),(3)u(S(t),t)=0.
The location of the liquid/solid interface is given by the heat balance equation known as Stefan condition
(4)-ux(S(t),t)=dS(t)dt,(5)u(x,t)=0,inΩS.
Initially there is liquid region which implies the condition
(6)S(0)=1,u(x,0)=0.
Referring to [14], this problem has the following exact solution for the liquid temperature distribution u(x,t) and the interface location S(t), respectively:
(7)u(x,t)=exp⁡(t-x+1)-1,inΩL,(8)S(t)=t+1,
where T=1 and plug f(t)=exp⁡(t+1) into (2). The exact solution (7) can be used to compare the numerical solutions and can also be used to initialize the numerical schemes.A VSG method, proposed in Kutluay et al. [14], which uses a fixed node number between a fixed boundary (x=0) and a moving boundary (x=s(t)), was kept constant and equal to N so that the moving boundary always lays on the Nth grid. They derive the equation and the boundary conditions with respect to u(x,t) as follows. First, the process finds the difference between the S(tm) and S(tm+1) and divides it by N, the number of nodes. In the case below, they assign N+1 as the number of x and M+1 as the number of t. Therefore, Δt=T/M so that tm=mΔt where T=1 is the maximum of t. Similarly, Δx=S(t)/N so that xn=nΔx. We also define unm=u(xn,tm), Sm=S(tm), S′=dS/dt, and xnm=xn(tm). Consider
(9)dS=Sm+1-Sm,dxn=xnm+1-xnmn=0,1,2,…,N,m=0,1,2,…,M-1.
The ratio dS:dxn equals the ratio Sm:xnm. Thus,
(10)dxndt=xnmSmdSdtut|n=ux|tx′(t)|n+ut|x∴ut=uxxnSdSdt+uxx.
Since -ux(0,t)=f(t),
(11)u-1m=u1m+2Δxmf(t).
When f(t)=exp⁡(t+1),
(12)u-1m=u1m+2Δxmexp⁡(tm+1).
By discretizing (10), they get
(13)unm+1=unm+xnmS′Δt2SmΔxm(un+1m-un-1m)+Δt(Δxm)2(un+1m-2unm+un-1m).
Plugging (12) into (13), they obtain the following:
(14)u0m+1=(1-2Δt(Δxm)2)u0m+2Δt(Δxm)2u1m+(2ΔtΔxm-ximS′ΔtSm)exp⁡(tm+1).
In order to get S(t), they use
(15)Sm+1=Sm-Δt2Δxm(3uNm-4uN-1m+uN-2m)wherem=0,1,2,….A boundary immobilization (BIM) method, proposed in Landau [13], uses a change in variables to transform the free boundary problem into a nonlinear problem on a fixed domain. The following transformation of state variable serves for such a purpose. The main point of BIM method is transformation. The transform function ξ=x/S(t) changes any shape of domain to a rectangular domain. ξ=0 when x=0, and ξ=1 when x=S(t). This fixes the free boundary to ξ=1. In the case below, he assigns N+1 as the number of ξ(0≤ξ≤1) and M+1 as the number of t. Therefore, Δt=T/M so that tm=mΔt where T is the maximum of t. Similarly, Δξ=1/N so that ξn=nΔξ. He also define unm=u(ξn,tm), Sm=S(tm), and S′=dS/dt(16)ut|x=uξξt+ut|ξ=-xS2S′(t)+ut|ξ.
Thus,
(17)ut=ξSS′uξ+1S2uξξ,0<ξ<1,t>0.
When ξ=0,
(18)uξ=-Sexp⁡(t+1).
By using central scheme, he gets
(19)u-1m=u1m+2ΔξSmexp⁡(tm+1).
Discretizing, he obtains
(20)unm+1=unm+ξnS′Δt2SmΔξ(un+1m-un-1m)+ΔtSm2Δξ2(un+1m-2unm+un-1m).
Plugging (19) into (20) when n=0, he gets
(21)u0m+1=(1-2ΔtSm2Δξ2)u0m+2ΔtSm2Δξ2u1m+(2ΔtSmΔξ-S′ξ0Δt)exp⁡(tm+1).
He uses the following equation in order to trace Sm:
(22)Sm+1=Sm-Δt2Δξ(3uNm-4uN-1m+uN-2m)wherem=0,1,2,…,M-1.


## 3. Log-Transform Function

In this section, we present a log-transform function with the unfixed free boundary and fixed free boundary that can determine the free boundary by solving a quadratic equation in a time-recursive way. Under the assumption of the Stefan problem, the time for free boundary can be shown to be the first hitting time of a boundary, the free boundary, in the plane consisting of pairs of the space and the time to expiration. Namely, the temperature curve of Stefan problem touches the line representing the intrinsic value tangentially. With a careful examination of the solution surface near the free boundary, we find a Lipschitz surface which reduces the accumulating error of the solution surface near free boundary. To find the free boundary, we present a log-transform function with the unfixed and fixed free boundary that has Lipschitz character near the free boundary as follows:(23)w(x,t)∶=ln⁡⁡(u(x,t)+1).
The transformed function w(x,t) provides that the solution surface in ΩS is a horizontal plane, and it is an inclined plain in ΩL. Namely, this function forms a gradual slope with the hyperplane corresponding to the solid region, thereby making the interpolation error more easily reducible when finding the free boundary (see Figure 1). w(x,t) also has a Lipschitz character with nonsingularity in (0,S(t)). More precisely, we have |S′(t)|/(f(t)+1)<|wx(S(t),t)|=|wx(S(t),t)/(w(S(t),t)+1)|=|-S′(t)/(w(S(t),t)+1)|≤|S′(t)|<∞. Then, w(x,t) is Lipschitz continuous and a natural candidate function for computation in [0,S(t)]. We also calculate the partial derivative with respect to x in (10), (17).Log-transformation fromu(x,t) to w(x,t).
(a)
(b)Figure1(a) shows that u(x,t) is transformed to w(x,t) and Figure 1(b) shows reduced interpolation error near the free boundary.We find the log-transform function with the unfixed and fixed free boundary to decide the free boundary by the Taylor series.


### 3.1. Log-Transform Function with Unfixed Free Boundary (LTUF)

From (1), (3), and (4), we obtain the following relations near free boundary:
(24)u(S(t),t)=0ux(S(t),t)=-S′(t)uxx(S(t),t)=(1+xNS)S′(t)2.
From u(x,t)=ew(x,t)-1, we obtain the following relations near free boundary (w(S(t),t)=0):
(25)ut=wt,ux=wx,uxx=wx2+wxx.
Plugging (25) into (10), we obtain
(26)wt=xiS(t)S′(t)wx+wx2+wxx.
Using u(x,t)=ew(x,t)-1, (24), and (25), we obtain the following:
(27)w(S(t),t)=0wx(S(t),t)=-S′(t)wxx(S(t),t)=xNSS′(t)2.
We lower the error by using the optimized safety parameter ρ for a Taylor series with moving node, a method suggested by Kim et al. [15], in order to control the error:
(28)w(x,t)=w(S(t),t)+wx(S(t),t)1!(x-S(t))+wxx(S(t),t)2!(x-S(t))2+O((x-S(t))3),(29)w(S~(t),t)=w(S(t),t)+wx(S(t),t)1!(S~(t)-S(t))+wxx(S(t),t)2!(S~(t)-S(t))2.
Using S~(t)=S(t)-ρΔx, (27), and (29), we get
(30)xN2S(-ρΔx)2S′(t)2+ρΔxS′(t)-w~=0,vvw~=w(S~(t),t).
Combining (30) with ux(S(t),t)<0, S′(t)>0, we have
(31)S′(t)=-b+b2-4ac2a,
where
(32)a=xN2S(-ρΔx)2,b=ρΔx,c=-w~.
For discretization (Δx,Δt), we introduce a two-dimensional mesh in the first quadrant of the x-t plane. From (31) we have
(33)Sm+1=Sm+-b~+b~2-4a~c~2a~Δt,
where
(34)a~=xNm2Sm(-ρΔxm)2,b~=ρΔxm,c~=-w(NΔx-ρΔx,mΔt),vvvjvm=0,1,2,…,M-1.
When the initial values are given by w(NΔx-ρΔx,t0) (transformed value of the temperature) and S0 (free boundary) at t=0, we can determine S1 (free boundary) at Δt using (33) and find w from (26). More importantly, for updating the free boundary our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until MΔt and obtain the free boundary in a time-recursive way.


### 3.2.  Log-Transform Function with Fixed Free Boundary (LTFF)

In the case below, we assignN+1 as the number of y(0≤y≤1) and M+1 as the number of t. Therefore, Δt=T/M so that tm=mΔt, where T is the maximum of t. Similarly, Δy=1/N so that yn=nΔy. We also define u^nm=u^(yn,tm), Sm=S(tm), and S′=dS/dt. The method starts with
(35)y=xS(t),u^(y,t)=u(x,t).
We obtain from (17) and (35)
(36)u^(1,t)=0,u^y(1,t)=-S(t)S′(t),u^yy(1,t)=2S(t)2S′(t)2.
From u^(y,t)=exp⁡(w(y,t))-1, we obtain the following relations near free boundary (w(1,t)=0):
(37)u^t=wt,u^y=wy,u^yy=wy2+wyy.
Plugging (35) and (37) into (17), we obtain
(38)wt=ξSS′wξ+1S2wξ2+1S2wξξ.
Using u^(x,t)=exp⁡(w(x,t))-1, (36), and (37), we obtain the following:
(39)w(1,t)=0,wy(1,t)=-S(t)S′(t),wyy(1,t)=S(t)2S′(t)2.
We lower the error by using the optimized safety parameter ρ for a Taylor series with fixed free boundary, a method suggested by Kim et al. [16], in order to control the error. We use Taylor series so that we can manage the error by employing the safety parameter ρ:
(40)w(y,t)=w(1,t)+wy(1,t)1!(y-1)+wyy(1,t)2!(y-1)2+O((y-1)3),w(y~,t)=w(1,t)+wy(1,t)1!(y~-1)+wyy(1,t)2!(y~-1)2.
Using y~=1-ρΔy, (39), and (40), we get
(41)12S(t)2(-ρΔy)2S′(t)2+S(t)ρΔyS′(t)-w~=0.
Combining (41) with ux(S(t),t)<0, S′(t)>0, we have
(42)S′(t)=-b+b2-4ac2a,
where
(43)a=12S(t)2(-ρΔy)2,b=S(t)ρΔy,c=-w~.
For discretization (Δy,Δt), we introduce a two-dimensional mesh in the first quadrant of the y-t plane. From (42) we have
(44)Sm+1=Sm+-b~+b~2-4a~c~2a~Δt,
where
(45)a~=12Sm2(-ρΔy)2,b~=SmρΔy,c~=-w(1-ρΔy,mΔt),jjvm=0,1,2,…,M-1.
When the initial values are given by w(1-ρΔy,t0) (transformed value of the temperature) and S0 (free boundary) at t=0, we can determine S1 (free boundary) at Δt using (44) and find w from (38). More importantly, for updating the free boundary, our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until MΔt and obtain the free boundary in a time-recursive way.


## 3.1. Log-Transform Function with Unfixed Free Boundary (LTUF)

From (1), (3), and (4), we obtain the following relations near free boundary:
(24)u(S(t),t)=0ux(S(t),t)=-S′(t)uxx(S(t),t)=(1+xNS)S′(t)2.
From u(x,t)=ew(x,t)-1, we obtain the following relations near free boundary (w(S(t),t)=0):
(25)ut=wt,ux=wx,uxx=wx2+wxx.
Plugging (25) into (10), we obtain
(26)wt=xiS(t)S′(t)wx+wx2+wxx.
Using u(x,t)=ew(x,t)-1, (24), and (25), we obtain the following:
(27)w(S(t),t)=0wx(S(t),t)=-S′(t)wxx(S(t),t)=xNSS′(t)2.
We lower the error by using the optimized safety parameter ρ for a Taylor series with moving node, a method suggested by Kim et al. [15], in order to control the error:
(28)w(x,t)=w(S(t),t)+wx(S(t),t)1!(x-S(t))+wxx(S(t),t)2!(x-S(t))2+O((x-S(t))3),(29)w(S~(t),t)=w(S(t),t)+wx(S(t),t)1!(S~(t)-S(t))+wxx(S(t),t)2!(S~(t)-S(t))2.
Using S~(t)=S(t)-ρΔx, (27), and (29), we get
(30)xN2S(-ρΔx)2S′(t)2+ρΔxS′(t)-w~=0,vvw~=w(S~(t),t).
Combining (30) with ux(S(t),t)<0, S′(t)>0, we have
(31)S′(t)=-b+b2-4ac2a,
where
(32)a=xN2S(-ρΔx)2,b=ρΔx,c=-w~.
For discretization (Δx,Δt), we introduce a two-dimensional mesh in the first quadrant of the x-t plane. From (31) we have
(33)Sm+1=Sm+-b~+b~2-4a~c~2a~Δt,
where
(34)a~=xNm2Sm(-ρΔxm)2,b~=ρΔxm,c~=-w(NΔx-ρΔx,mΔt),vvvjvm=0,1,2,…,M-1.
When the initial values are given by w(NΔx-ρΔx,t0) (transformed value of the temperature) and S0 (free boundary) at t=0, we can determine S1 (free boundary) at Δt using (33) and find w from (26). More importantly, for updating the free boundary our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until MΔt and obtain the free boundary in a time-recursive way.


## 3.2.  Log-Transform Function with Fixed Free Boundary (LTFF)

In the case below, we assignN+1 as the number of y(0≤y≤1) and M+1 as the number of t. Therefore, Δt=T/M so that tm=mΔt, where T is the maximum of t. Similarly, Δy=1/N so that yn=nΔy. We also define u^nm=u^(yn,tm), Sm=S(tm), and S′=dS/dt. The method starts with
(35)y=xS(t),u^(y,t)=u(x,t).
We obtain from (17) and (35)
(36)u^(1,t)=0,u^y(1,t)=-S(t)S′(t),u^yy(1,t)=2S(t)2S′(t)2.
From u^(y,t)=exp⁡(w(y,t))-1, we obtain the following relations near free boundary (w(1,t)=0):
(37)u^t=wt,u^y=wy,u^yy=wy2+wyy.
Plugging (35) and (37) into (17), we obtain
(38)wt=ξSS′wξ+1S2wξ2+1S2wξξ.
Using u^(x,t)=exp⁡(w(x,t))-1, (36), and (37), we obtain the following:
(39)w(1,t)=0,wy(1,t)=-S(t)S′(t),wyy(1,t)=S(t)2S′(t)2.
We lower the error by using the optimized safety parameter ρ for a Taylor series with fixed free boundary, a method suggested by Kim et al. [16], in order to control the error. We use Taylor series so that we can manage the error by employing the safety parameter ρ:
(40)w(y,t)=w(1,t)+wy(1,t)1!(y-1)+wyy(1,t)2!(y-1)2+O((y-1)3),w(y~,t)=w(1,t)+wy(1,t)1!(y~-1)+wyy(1,t)2!(y~-1)2.
Using y~=1-ρΔy, (39), and (40), we get
(41)12S(t)2(-ρΔy)2S′(t)2+S(t)ρΔyS′(t)-w~=0.
Combining (41) with ux(S(t),t)<0, S′(t)>0, we have
(42)S′(t)=-b+b2-4ac2a,
where
(43)a=12S(t)2(-ρΔy)2,b=S(t)ρΔy,c=-w~.
For discretization (Δy,Δt), we introduce a two-dimensional mesh in the first quadrant of the y-t plane. From (42) we have
(44)Sm+1=Sm+-b~+b~2-4a~c~2a~Δt,
where
(45)a~=12Sm2(-ρΔy)2,b~=SmρΔy,c~=-w(1-ρΔy,mΔt),jjvm=0,1,2,…,M-1.
When the initial values are given by w(1-ρΔy,t0) (transformed value of the temperature) and S0 (free boundary) at t=0, we can determine S1 (free boundary) at Δt using (44) and find w from (38). More importantly, for updating the free boundary, our method does not include any iteration until the numerical solution is obtained. So, we repeat the previously mentioned process until MΔt and obtain the free boundary in a time-recursive way.


## 4. Numerical Result

In this section, we provide two numerical examples to illustrate our method. We also provide runtimes and computation errors used to compare with the results obtained by other numerical methods such as the VSG method developed by Kutluay et al. [14] and the BIM method developed by Landau [13]. In this section, we compare the four methods based on runtime, L2 norm, and root mean square error (RMSE). All calculations are carried out using a C++ implementation with the a 2.66 Ghz Intel 4 Core CPU with 4 GB RAM. A finite difference method with an explicit finite difference approximation is proposed for our methods. The benchmark results are obtained using the exact solution.The parameter values used to calculate the free boundary and the values of temperature areT=1.0, f(t)=exp⁡(t+1), S(0)=1, and a discrete mesh of 10000×40 nodes.In Figure2, we find a numerical optimization ρ=0.00706. Table 1 also shows the results of the free boundary and the values of temperature with the safety parameters. One can see from Table 1 that the free boundary monotonically increases as the value of ρ increases, but the values of temperature monotonically decrease when ρ increases. They are so gradual that they are not very susceptible to change in the value of ρ.Table 1
Values in the solution.
Safety parameterρ
RMSE
S
(
t
=
0.5
)
u
(
0.5
,
T
)
u
(
1.0
,
T
)
u
(
1.5
,
T
)
0.0060
0.000407836
1.4998795
3.4807649
1.7177973
0.6484301
0.0062
0.000407824
1.4998786
3.4807643
1.7177965
0.6484289
0.0064
0.000407814
1.4998777
3.4807638
1.7177957
0.6484277
0.0066
0.000407806
1.4998767
3.4807632
1.7177949
0.6484264
0.0068
0.000407802
1.4998758
3.4807627
1.7177941
0.6484252
0.0070
0.000407799
1.4998748
3.4807621
1.7177933
0.6484240
0.0072
0.000407800
1.4998739
3.4807616
1.7177925
0.6484227
0.0074
0.000407803
1.4998729
3.4807610
1.7177917
0.6484215
0.0076
0.000407809
1.4998720
3.4807604
1.7177909
0.6484203
0.0078
0.000407817
1.4998711
3.4807599
1.7177901
0.6484190
0.0080
0.000407828
1.4998702
3.4807594
1.7177893
0.6484178
Maximum difference
0.000000037
0.0000093
0.0000055
0.0000080
0.0000123Figure 2
RMSE of the safety parameterρ in the solution.In Figure3, Table 2, Figure 4, and Table 3, we take the parameter values used in Figure 2 except for the discrete mesh, plot runtimes, and computational errors compared with the various methods.Table 2
Comparison of RMSE and runtime in the solution.
Mesh size(N×M)
(Runtime, RMSE)
VSG
BIM
LTUF
LTFF
1250
×
5
(0.045, 0.02346430)
(0.086, 0.02214381)
(0.030, 0.02340634)
(0.057, 0.02206451)
2500
×
10
(0.082, 0.00589595)
(0.091, 0.00552885)
(0.083, 0.00588832)
(0.098, 0.00551475)
5000
×
20
(0.203, 0.00152029)
(0.296, 0.00142859)
(0.227, 0.00151975)
(0.289, 0.00142755)
10000
×
40
(0.820, 0.00042633)
(1.005, 0.00040835)
(0.867, 0.00042561)
(1.037, 0.00040779)
20000
×
80
(3.640, 0.00010951)
(4.122, 0.00010726)
(3.700, 0.00010865)
(3.969, 0.00009299)Table 3
Comparison of L2 norm and runtime in the solution.
Mesh size(N×M)
(Runtime, L2 norm)
VSG
BIM
LTUF
LTFF
1250
×
5
(0.045, 0.33100481)
(0.086, 0.31159113)
(0.030, 0.33101571)
(0.057, 0.31125825)
2500
×
10
(0.082, 0.08358954)
(0.091, 0.07818977)
(0.083, 0.08348142)
(0.098, 0.07799352)
5000
×
20
(0.203, 0.02155389)
(0.296, 0.02020339)
(0.227, 0.02154625)
(0.289, 0.02018861)
10000
×
40
(0.820, 0.00604438)
(1.005, 0.00577497)
(0.867, 0.00603415)
(1.037, 0.00576715)
20000
×
80
(3.640, 0.00155269)
(4.122, 0.00151694)
(3.700, 0.00154049)
(3.969, 0.00131512)Figure 3
Comparison of RMSE and runtime in the solution.Figure 4
Comparison of L2 norm and runtime in the solution.Figure5 and Table 4 show the RMSE and L2 norm error of the LTUF method and compare it with sufficiently inclined slope.Table 4
Comparison of RMSE and L2 norm with LTUF method in the solution.
Mesh size(N×M)
LTUF
(Runtime, RMSE)
(Runtime, L2 norm)
1250
×
5
(0.030, 0.02340634)
(0.030, 0.33101571)
2500
×
10
(0.083, 0.00588832)
(0.083, 0.08348142)
5000
×
20
(0.227, 0.00151975)
(0.227, 0.02154625)
10000
×
40
(0.867, 0.00042561)
(0.867, 0.00603415)
20000
×
80
(3.700, 0.00010865)
(3.700, 0.00154049)Figure 5
Comparison of RMSE and L2 norm with LTUF method in the solution.Figure6 and Table 5 show the RMSE and L2 norm error of the LTFF method and compare it with sufficiently inclined slope. Especially, Figures 5 and 6 and Tables 4 and 5 show the numerical convergence of our methods. Note that the discrete meshes of 1250×5, 2500×10, 5000×20, 10000×40, and 20000×80 nodes are plotted in Figures 3, 4, 5, 6, 7, 8, and 9, and Tables 2, 3, 4, 5, 6, and 7.Table 5
Comparison of RMSE and L2 norm with LTFF method in the solution.
Mesh size(N×M)
LTFF
(Runtime, RMSE)
(Runtime, L2 norm)
1250
×
5
(0.057, 0.02206451)
(0.057, 0.31125825)
2500
×
10
(0.098, 0.00551475)
(0.098, 0.07799352)
5000
×
20
(0.289, 0.00142755)
(0.289, 0.02018861)
10000
×
40
(1.037, 0.00040779)
(1.037, 0.00576715)
20000
×
80
(3.969, 0.00009299)
(3.969, 0.00131512)Table 6
Comparison of VSG, BIM, LTUF, LTFF method, and exact solution whenf(t)=exp⁡(t+1).
Positionx
Numerical solution
Exact solution
VSG
BIM
LTUF
LTFF
0.00
6.389056
6.389056
6.389056
6.389056
6.389056
0.10
5.685461
5.685448
5.685467
5.685440
5.685894
0.20
5.049275
5.049261
5.049281
5.049254
5.049647
0.30
4.473626
4.473613
4.473632
4.473605
4.473947
0.40
3.952754
3.952741
3.952761
3.952734
3.953032
0.50
3.481448
3.481434
3.481454
3.481428
3.481689
0.60
3.054989
3.054975
3.054997
3.054978
3.055199
0.70
2.669112
2.669098
2.669120
2.669094
2.669296
0.80
2.319954
2.319939
2.319962
2.319937
2.320117
0.90
2.004020
2.004006
2.004030
2.004006
2.004116
1.00
1.718151
1.718137
1.718162
1.718139
1.718281
1.10
1.459484
1.459470
1.459498
1.459474
1.459603
1.20
1.225432
1.225418
1.225446
1.225425
1.225540
1.30
1.013652
1.013638
1.013667
1.013648
1.013752
1.40
0.822024
0.822010
0.822042
0.822024
0.822119
1.50
0.648631
0.648618
0.648651
0.648635
0.648721
1.60
0.491738
0.491725
0.491760
0.491746
0.491824
1.70
0.349775
0.349762
0.349799
0.349787
0.349859
1.80
0.221320
0.221308
0.221347
0.221338
0.221403
1.90
0.105089
0.105077
0.105119
0.105112
0.105171
2.00
0.000000
0.000000
0.000000
0.000000
0.000000
RMSE
0.000109
0.000107
0.000108
0.000092
L2 norm
0.001552
0.001516
0.001540
0.001315
Runtime (secs)
3.640
4.122
3.700
3.969Table 7
Comparison of VSG, BIM, LTUF, LTFF method, and exact solution whenf(t)=1.
positionx
Numerical solution
Exactsolution
VSG
BIM
LTUF
LTFF
0.00
1.128379
1.128379
1.128379
1.128379
1.128379
0.10
1.030558
1.030846
1.031035
1.031171
1.031198
0.20
0.939323
0.939562
0.939461
0.938557
0.939644
0.30
0.853239
0.853118
0.852186
0.853340
0.853672
0.40
0.772252
0.772710
0.773105
0.772877
0.773215
0.50
0.696855
0.697913
0.698005
0.697264
0.698177
0.60
0.627581
0.628210
0.627142
0.628511
0.628436
0.70
0.563013
0.562986
0.563804
0.562138
0.563851
0.80
0.503018
0.503362
0.504110
0.504430
0.504255
0.90
0.447454
0.449016
0.448385
0.448748
0.448465
1.00
0.397071
0.398605
0.394381
0.399291
0.399282
1.10
0.351200
0.352020
0.353327
0.353470
0.353492
1.20
0.309050
0.310587
0.310942
0.311511
0.311870
1.30
0.270431
0.272901
0.274082
0.274519
0.274184
1.40
0.235147
0.238341
0.239877
0.239340
0.240196
1.50
0.203922
0.206799
0.208701
0.210078
0.209664
1.60
0.175475
0.179454
0.181921
0.182140
0.182347
1.70
0.149508
0.154482
0.157249
0.158290
0.158005
1.80
0.125812
0.131752
0.135028
0.136610
0.136403
1.90
0.104323
0.111398
0.116116
0.117326
0.117311
2.00
0.085031
0.093533
0.098785
0.100843
0.100509
RMSE
0.0055537
0.0024662
0.0007687
0.0006550
L2 norm
0.0789335
0.0350518
0.0109265
0.0093102
Runtime (secs)
0.795
0.860
0.783
0.821Figure 6
Comparison of RMSE and L2 norm with LTFF method in the solution.Figure 7
Comparison of RMSE and L2 Norm in the solution.Figure 8
Increasing size of differential error of VSG and LTUF according to final timeT in the solution.Figure 9
Increasing size of differential error of BIM and LTFF according to final timeT in the solution.For the sake of clarity, the pairs of (RMSE, L2 norm) of all four methods are graphed in Figure7. Figure 7 illustrates the pairs in the solution. The figure shows that, in terms of errors (RMSE, L2 norm), LTUF approximates more accurately than VSG does, and LTFF approximates more accurately than BIM does. In Figures 8 and 9, we investigate the functional behaviors of the difference between RMSE (VSG) and RMSE (LTUF), L2 norm (VSG) and L2 norm (LTUF), RMSE (BIM) and RMSE (LTFF), and L2 norm (BIM) and L2 norm (LTFF) with respect to final time change.Table6 reports the values of the temperatures for the specific parameter set associated with the discrete mesh of 20000×80.In Figure10 and Table 7, we investigate the free boundary behavior with respect to time change when T=1.0,f(t)=1, and S(0)=0. From (1)~(5) this problem has the following exact solution u(x,t):
(46)u(x,t)=2tπexp⁡(-x24t)-xerfc⁡(x2t).Numerical solution and contours whenf(t)=1.
(a)
(b)


## 5. Final Remarks

We suggested LTUF and LTFF methods. There are approximated solutions that show reasonable agreement.VSG method suggested by Kutluay et al. [14] shows that accumulating error occurred during calculation by the three-term backward difference near the free boundary, while LTUF method adopts a log-transform function to reduce the accumulating error in solution surface that causes gradual slope. Our method employing a log-transform function with the fixed free boundary solves a nonlinear problem on a fixed domain derived from a free boundary problem. Since the computation process depends on Lipschitz surface, we need to focus on the motion of the solution surface which would be simple to see the minute behavior of solution surface.Kutluay et al. [14] use three-term backward difference scheme during calculating the free boundary with BIM method suggested by Landau [13]. However, we emphasize that there is no three-term backward difference scheme in LTFF method. Instead, we solve the quadratic equation and directly find the free boundary. In an environment with the rapidly changing free boundary, our straight forward method is a very powerful tool to understand Stefan problem.Numerical study also shows that overall speed and accuracy comparisons have demonstrated some superiority of our methods over other methods. Further investigations on the extending of the various models could be interesting and useful. Also, more comparative and extensive studies on our method should be helpful for completing the mission of the present work. This remains as our future research work.



---

*Source: 764532-2014-06-03.xml*