TVE-F 19007
Examensarbete 15 hp Juni 2019
Solution of the Stefan problem
with general time-dependent boundary conditions using a random walk
method
Daniel Stoor
Teknisk- naturvetenskaplig fakultet UTH-enheten
Besöksadress:
Ångströmlaboratoriet Lägerhyddsvägen 1 Hus 4, Plan 0
Postadress:
Box 536 751 21 Uppsala
Telefon:
018 – 471 30 03
Telefax:
018 – 471 30 00
Hemsida:
http://www.teknat.uu.se/student
Abstract
Solution of the Stefan problem with general
time-dependent boundary conditions using a random walk method
Daniel Stoor
This work deals with the one-dimensional Stefan problem with a general time- dependent boundary condition at the fixed boundary. The solution will be obtained using a discrete random walk method and the results will be compared qualitatively with analytical- and finite difference method solutions. A critical part has been to model the moving boundary with the random walk method. The results show that the random walk method is competitive in relation to the finite difference method and has its advantages in generality and low effort to implement. The finite difference method has, on the other hand, higher accuracy for the same computational time with the here chosen step lengths. For the random walk method to increase the accuracy, longer execution times are required, but since the method is generally easily adapted for parallel computing, it is possible to speed up. Regarding applications for the Stefan problem, there are a large range of examples such as climate models, the diffusion of lithium-ions in lithium-ion batteries and modelling steam chambers for oil extraction using steam assisted gravity drainage.
Ämnesgranskare: Maria Strömme Handledare: Magnus Ögren
Populärvetenskaplig sammanfattning
Ett Stefanproblem är ett värmeledningsproblem med en rörlig rand, vilket betyder att domänen där temperaturfördelningen utvärderas förändras med tiden. Som exempel på tillämpning av problemet kan tänkas ett smältande isblock där man önskar bestämma smältvattnets och/eller isens temperatur samt hur mycket is som smält efter en viss tid. Problemet är uppkallat efter Josef Stefan som formulerade problemet i en artikel 1891 och har idag flera tillämpningsområden, såsom klimatmodeller, modellering av hur litiumjoner diffunderar i litiumjonbatterier och vid användning av metoder för oljeutvin- ning där ånga injiceras i marken. I Stefans originalartikel presenterades ex- akta lösningar för problemet förutsatt att temperaturen på domänens yta inte förändrades, men i många tillämpningar är det intressant att ta fram lösningar då temperaturen på domänens yta istället är tidsberoende. En naturlig tillämpning för detta kan tänkas vara smältande is där ytans tem- peratur förändras efter dygnstemperaturen. För att hitta lösningar i dessa fall är det nödvändigt att använda numeriska metoder.
Den numeriska metod som här undersöks är en slumpvandringsmetod och bygger på kopplingen mellan sannolikhetsteori och värmeledningsekvationen.
Värmeledning som fysikalisk process uppkommer genom termisk rörelse hos små partiklar i ett medium. Genom att modellera detta matematiskt med slumppartiklar som vandrar i olika riktingar med vissa sannolikheter kan man under vissa villkor få fram ekvationer som uppfyller värmeledningsek- vationen. Målsättningen med detta projekt har varit att anpassa en så- dan slumpvandringsmetod till Stefanproblemet för godtyckliga tidsberoende temperaturer på domänens yta. En viktig del har handlat om att modellera hur gränsen mellan fast- och flytande fas förändras med hjälp av slumpvan- dringsmetoden.
Resultaten har sedan jämförts med resultat från finita differensmetoden, som är en annan numerisk metod, och exakta lösningar. Här visar sig slumpvan- dringsmetoden ge goda resultat med rätt val av parametrar och i jämförelse med finita differensmetoden är den fördelaktig när det gäller flexibilitet och enkelhet att implementera. I detta arbete har endast ett endimensionellt fall studerats, men fördelarna för slumpvandringsmetoden gäller i ännu högre grad i högre dimensioner.
Contents
1 Introduction 5
1.1 Objectives . . . . 5
1.2 Background . . . . 6
1.3 Classic heat conduction problem . . . . 6
2 The Stefan problem in one dimension 7 2.1 The Stefan condition . . . . 8
3 Random walk modelling 10 3.1 Brownian motion as the limit of random walks . . . 12
3.2 Random walk and the heat equation . . . 14
3.3 A random walk model with boundary conditions . . . 16
3.4 Modelling the moving boundary . . . 17
4 Results 19 4.1 Heat equation with fixed boundary . . . 19
4.2 Stefan problem with constant boundary condition . . . 20
4.3 Stefan problem with boundary condition f(t) = et− 1 . . . . 21
4.4 Stefan problem with oscillating boundary condition . . . 22
4.5 Stefan problem with boundary condition according to daytime temperature variations . . . 23
5 Discussion 24 5.1 Applications . . . 24
6 Conclusion 25
References 26
Appendix 28
1 Introduction
In many applications in science and technology there are complicated par- tial differential equations (PDE:s) that need to be solved. Since analytical solutions may not always exist, different types of numerical methods are necessary. One of these methods is the random walk method (from here ab- breviated RWM), which has seen an increase in physical applications since probabilistic methods have been extended to fit a larger number of physical problems [1]. Among the strengths of the RWM is that it is very suitable for programming since it is easy to implement and does not create any large linear system to be solved. It often also supports parallel computing well.
Another benefit is that the RWM is a local method, i.e. a solution can be carried out at specific points without having to obtain a solution for the whole domain [1] [2]. At last the RWM is also suitable for higher dimensions since complicated geometries are not as hard to model and memory usage is lower than with e.g. the finite difference method. In this work though, we will be restricted to the one-dimensional case.
We will here use a discrete RWM to examine solutions to the Stefan problem, which is a PDE consisting of the heat equation defined in a phase changing medium. There are different types of formulations of this problem, but one of the characteristics is that it has a free or moving boundary governed by the Stefan condition, which describes the position of the interface between the phases. Beyond the moving boundary, the general formulation of the problem usually also includes a fixed boundary with a boundary condition different from the moving one. For physical reasons the boundary condition at the moving boundary is set to the melting point TM, but at the fixed boundary the conditions may be set to an arbitrary function f(t). For the case where we have a constant temperature at the fixed boundary f(t) = T0
and for one other particular form of f(t), there are analytical solutions to the Stefan problem. Otherwise, in most cases we need numerical approximations to evaluate a solution. As an example of such a case, it would be meaningful to model the melting of a medium where the surface temperature is an oscillating function according to the variations in day temperature.
1.1 Objectives
The aim of this thesis is to solve the Stefan problem for arbitrary functions f (t) using the RWM. Moreover, we want to do a qualitative comparison analysis between the RWM, the exact solutions (when they exist) and finite difference method solutions made by Jonsson [3]. To make the choice of method comprehensible, we also want to present some theory about the connection between the RWM and the heat equation.
1.2 Background
The RWM has it roots in probabilistic problems such as "Gambler’s ruin"
described by e.g. Pascal [4]. An important contribution for the modern applications of the RWM was the formulation of Brownian motion that sat- isfies the diffusion equation in an article by Einstein in 1905 [5]. As we will see in section 3.1, random walk is closely associated with Brownian motion.
From 1920s RWM:s have been used for analyzing PDE:s including diffusion processes and equilibrium states and since probability theory and computer science have developed, RWM:s have been used in a broader range of ap- plications [1]. Among these are population genetics, description of financial markets and polymer chains. A recent study shows that diffusion patterns and RWM:s also can be applied when analyzing network structures, such as the internet [4].
The Stefan Problem has its name from Josef Stefan (1835-1893) who was first to describe problems including a moving boundary in detail. This was investigated in his report on ice formation in polar seas [6], where he also presents the analytical solution (2.8) to the problem where the fixed bound- ary has constant temperature. Although, for the general case where f(t) is an arbitrary function, Stefan could not give an explicit solution except for some ideas and approximations. An explicit solution for this case has not been obtained yet, though there are some attempts to formulations, e.g. [7], which may require effort to implement. Except for the ice formation prob- lem examined by Stefan, moving boundary problems have many applications which will be mentioned later in the discussion section 5.1. Since the mod- elling of moving boundary problems in these applications plays an important role and has not yet been evaluated with a RWM, this will be the subject of this thesis.
1.3 Classic heat conduction problem
As a starting problem, the classic heat equation with fixed boundary and homogenous Dirichlet conditions is considered.
∂T
∂t = α∂2T
∂x2 , 0 < x < L , t > 0 , (1.1)
T (0, t) = T (L, t) = 0 , t > 0 , (1.2)
T (x, 0) = g(x) . (1.3)
Using separation of variables, the general solution to this problem is [8]
T (x, t) =
∞
X
n=1
cne−α(nπL)2tsinnπx L
. (1.4)
We write the initial condition as
g(x) =
∞
X
n=1
cnsin
nπx L
. (1.5)
By recognizing this as the Fourier series expansion of g(x) on 0 < x < L, cn
can be determined.
cn= 2 L
L
Z
0
g(x) sin
nπx L
dx . (1.6)
Here we set g(x) = 1 which gives the solution on the form
T (x, t) = 4 π
∞
X
k=0
exp (−απL22(2k + 1)2t)
2k + 1 sin (2k + 1)πx L
. (1.7)
2 The Stefan problem in one dimension
In our model for the one-dimensional Stefan problem we consider a block of ice with infinite extent and one surface to air. At t = 0 there is no water phase and the ice phase is at Tice = 0°. But for t > 0 the ice has started melting and thus we have a water phase on top of the ice. The temperature of the surface, i.e. the interface between air and water for t > 0 is changing over time according to f(t) and to simulate a melting process, we assume f (t)≥ 0 ∀ t. This yields the following equations:
∂T
∂t = αL∂2T
∂x2 , 0 < x < s(t) , t > 0 , (2.1)
T (0, t) = f (t) , t > 0 , (2.2)
T (x, 0) = 0 , (2.3)
lρds
dt = kL ∂T
∂x x=s(t)
, t > 0 , (2.4)
s(0) = 0 , (2.5)
T (s(t), t) = TM = 0 , t > 0 . (2.6)
Here the thermal diffusivity in the liquid part, αL in (2.1), is defined as
αL= kL
ρLcL, (2.7)
where kL is the heat conductivity, ρL the density and cL the specific heat capacity in the liquid part. Note that these properties differ between the solid and liquid part, αL6= αS. But since T = 0° in the liquid part and the temperature distribution only is evaluated in the solid part, αS is not taken into consideration here. In equation (2.4) l is the specific latent heat and ρ is the density. Here it is assumed ρ = ρL = ρS which is not the real case, but makes the model easier to handle [3].
The analytical solution to the problem when f(t) = T0 is (see e.g. [3])
T (x, t) = T0
1−erf(x/(2
√αLt)) erf(λ)
s(t) = 2λ√
αLt β√
πλeλ2erf λ = T0,
(2.8)
where β = cL/l and erf (x) is the error function defined as √2π Rx
0 e−y2dy. A special case when an analytical solution also can be found is when f(t) = et− 1. We will use this for comparison with our RWM in the result section.
Provided β = 1, the solution is [9]
T (x, t) = et−x− 1
s(t) = t . (2.9)
2.1 The Stefan condition
The free boundary, i.e. the interface between the two phases, is time- dependent and denoted as x = s(t). At time t0 the domain A = {x > 0}
is divided into two subdomains consisting of A1 = A∩ {x < s(t0)} and A2 = A∩ {x > s(t0)}, where A1 is the water phase and A2 the ice phase.
Here we will consider a one phase problem which means that the temper- ature in one of the phases (here the ice phase) is constant at the critical temperature T0 = 0°. A more general study would be a two phase problem where the temperature distribution in both phases are taken into account [10].
We will here briefly derive the Stefan condition found in (2.4), which will later also be useful in the numerical estimation of the interface s(t). More details on the derivation of the Stefan condition can be found in e.g. [10] and [3]. In the case of melting ice, the water phase at time t1 > t0 will increase resulting in s(t1) > s(t0). If we think of a block of ice with cross sectional
S
s(t0) s(t1)
ex
−→
water ice
Figure 2.1: Volume of melted ice between time t0 and t1.
area S, the volume of the melted ice on the time interval t0 < t < t1 is S × (s(t1) − s(t0)), see fig. 2.1. The heat [J] required for the melting is determined according to
Q = V · ρ · l = S(s(t1)− s(t0))· ρl , (2.10) where l is the latent heat. As we assume the heat is only spread by diffusion, the heat transport will obey Fourier’s law:
q =−ki
dT
dx , (2.11)
where q is the local heat flux density [W/m2]. By energy conservation and the expressions for the heat fluxes from the liquid and solid part, Q can be formulated as:
Q =
t1
Z
t0
Z
S
−kL
∂T (s(τ ), τ )
∂x · ex− kS
∂T (s(τ ), τ )
∂x · (−ex)
dA dτ
= S
t1
Z
t0
−kL
∂T (s(τ ), τ )
∂x + kS∂T (s(τ ), τ )
∂x
dτ . (2.12)
Putting equation (2.10) and (2.12) together, dividing by t1− t0 and letting t1→ t0 will yield the equation (2.4) for the Stefan condition:
ρl lim
t1→t0
S(s(t1)− s(t0)) t1− t0
= lim
t1→t0
1 t1− t0
S
t1
Z
t0
−kL
∂T (s(τ ), τ )
∂x + kS∂T (s(τ ), τ )
∂x
dτ (2.13)
⇒ ρlds
dt =−kL
∂T (s(t), t)
∂x + kS∂T (s(t), t)
∂x . (2.14)
Here t0 is replaced by t since t0 can be chosen arbitrarily. In our case where we assume T = 0° for x > s(t), diffusion only occur in the liquid phase and (2.14) reduces to [3]
ρlds
dt =−kL
∂T (s(t), t)
∂x . (2.15)
3 Random walk modelling
The basic idea of the one-dimensional RWM is to let a ’walker’ start at the origin and by a certain probability take either one step to the right or one step to the left. For two or three dimensions there will be four or six possible paths. Further on we will see that this model with proper choices of param- eters will be suitable for the diffusion process, which governs heat transport.
But to get there, a presentation on the connection between the RWM and the diffusion equation will be made.
At first, consider a random walker moving on a one-dimensional integer lat- tice. If Sn denotes the position after n steps starting at x0 and Xj =±1 is a random step in position between two time steps tj−1 and tj, we define
Sn= x0+
n
X
j=1
Xj, (3.1)
where both results of Xj are here chosen to be equally probable, P {Xj = 1} = P {Xj =−1} = 1/2, see fig. 3.1. To see how walkers disperse with time we will now investigate how far from the origin (x0= 0) a walker have gone after n time steps. It is natural to assume that the average position or mean
Figure 3.1: Random walk on the integer lattice.
distance with many walkers will be zero since the probability for a move to the left and the right are equal. But to measure the average distance from the origin after n steps we use the variance σ2or mean square distance hSn2i.
The definition is
σ2 =hS2ni =
n
X
i=1
pi(xi− µ)2, (3.2) with pi denoting the probability to end up at position xi. Since the mean distance is µ = 0 we get
σ2 =hSn2i =
n
X
i=1
pix2i . (3.3)
Here we will use a simple proof by Feynman [11] to show that hSn2i = n.
For n = 1, hS12i is easily obtained:
hS12i = 1
2· (−1)2+1
2 · 12 = 1 . (3.4)
For n > 1 we get Sn2 by starting with Sn−1. The distance after n − 1 steps is Sn−1 and thus after n steps the distance is Sn−1− 1 or Sn−1+ 1. Taking the square leads us to
Sn2 =
(Sn−1− 1)2 = Sn−12 − 2Sn−1+ 1 or
(Sn−1+ 1)2 = Sn−12 + 2Sn−1+ 1 .
(3.5)
Since both cases are equally probable, the mean square distance calculated as in (3.5), will be
hSn2i = 1
2(hSn−12 i−h2Sn−1i+1)+1
2(hSn−12 i+h2Sn−1i+1) = hSn−12 i+1 . (3.6) By induction it follows that
hS12i = 1
hS22i = hS12i + 1 = 2 ...
hSn2i = hSn−12 i + 1 = n .
(3.7)
3.1 Brownian motion as the limit of random walks
The objective for our RWM is to simulate Brownian motion. While random walk is time discrete the Brownian motion is a continuous random movement which describes the motion of small particles caused by thermal motion of atoms. This is the physical process behind diffusion and heat transport and thus the goal for our model. Standard Brownian motion is mathematically described as a Wiener process, which is a collection of random variables {Wt} with time index t. One of the important characteristics of a Wiener process Wt is that for each t, the random variable Wtis N(0, t), i.e. it has a normal distribution with mean zero and variance t. From the Wiener process we can define the general Brownian motion, X(t), by
X(t) = σWt+ µ (3.8)
where X(t) is N(µ, σ2t) [12]. To make the RWM approximate Brownian motion, we will start by giving it the property of the Wiener process, i.e.
N (0, t). As we will see this will set some requirements on the scaling between the step sizes ∆x and ∆t. In our first model, we had a random walk Sn on the integer lattice meaning ∆x = ∆t = 1. If we instead look at a specific time interval 0 < t < t0 we might discretize the interval with ∆t = δ = t0/N for a large integer N, see figure 3.2. Then we can evaluate Wtat t = jδ and compare the result with the random walk scaled with ∆x:
Wt= Wjδ ≈ ∆xSj. (3.9)
Figure 3.2: Random walks with t0 = 1 and N =100, 1000 and 10000 steps.
The question here is how to choose ∆x. According to the definition of the Wiener process we want the variance for the RWM to equal t0 for Wt0. The variance of the right side of (3.9) is
h(∆xSn)2i = (∆x)2hSn2i = (∆x)2N , (3.10)
where we have used the result in (3.7). Equaling the variance for the left and right side of (3.9) gives the right scaling for ∆x
t0 = (∆x)2N
⇒ ∆x =q
t0
N =√ δ =√
∆t . (3.11)
To normalize we set t0 = 1 and achieve t = jδ = j/N. Then we can reformulate (3.9) as
Wjδ = Wj/N ≈ Sj
√N . (3.12)
To prove our random walk on the right hand side in (3.12) is normal dis- tributed as in the definition of the Wiener process we can use the central limit theorem [13]
Theorem. (Central Limit Theorem, CLT) If X1, X2. . . are independent random variables with mean µ and variance σ2 the central limit theorem states that for large n the distribution of
Zn= (X1+ . . . + Xn)− nµ
√σ2n (3.13)
is approximately that of a normal distribution with mean 0 and variance 1.
More precisely for every a < b,
n→∞lim P{a < Zn< b} =
b
Z
a
√1
2πe−x2/2dx. (3.14) To use CLT we again rewrite (3.12)
Wjδ = Wj/N ≈ Sj
√N = Sj
√j r j
N = Sj
√j
√t . (3.15)
Since each random step Xj has the variance 1 we can say from CLT that as j→ ∞, Sj/√
jwill approach a normal distribution with mean 0 and variance 1. The CLT also tells us that if a distribution Y has the variance 1, σY will have the variance σ2. This yields that the right hand side of (3.15) is normal distributed with mean 0 variance t as j → ∞ and thus satisfies the definition for a Wiener process. Since the random walk Sj where j → ∞ is achieved by choosing δ → 0, we can summarize these steps by saying that a Wiener process is a limit of random walks when we choose a grid where ∆x =√
∆t and by letting ∆t → 0 [13]. For the general Brownian motion (3.8) we have to adjust the scaling with the variance term σ2, which according to (3.11) yields
∆x =√
σ2∆t . (3.16)
3.2 Random walk and the heat equation
In this section we will study the physical constants related to the heat equa- tion and see how to translate them into our RWM. In our one-dimensional model, we want to let one walker represent the temperature difference of 1
°C on the volume element ∆x × 1 × 1 m3. To make a simple illustration for the heat equation, consider the number of walkers in volume element i with width ∆x as Ni(t) at time t. If we let the probability for a walker to go either to the left or the right to be equal during a time step ∆t we say P = 1/2. Then we expect to have P · Ni walkers going to box i + 1 and the same amount going to box i − 1. At the same time walkers from boxes i − 1 and i + 1 will walk into box i giving the balancing equation for Ni [14]
Ni(t + ∆t) = Ni(t)− (walkers to i − 1) − (walkers to i + 1) +(walkers from i − 1) + (walkers from i + 1)
⇒ Ni(t + ∆t) = Ni(t)− P · Ni− P · Ni
+P · Ni−1+ P · Ni+1.
(3.17)
By making the substitution P = R∆t for some constant R we can rewrite this as
Ni(t + ∆t) = Ni(t) + R∆t(Ni+1(t)− 2Ni(t) + Ni−1(t)) . (3.18) Rearranging the terms gives
Ni(t + ∆t)− Ni(t)
∆t = R(∆x)2(Ni+1(t)− 2Ni(t) + Ni−1(t))
(∆x)2 , (3.19)
which is the same as a discretized partial differential equation for N(x, t):
∂N
∂t = R(∆x)2∂2N
∂x2 . (3.20)
We can see that (3.20) has the same appearance as the heat equation:
∂T
∂t = α∂2T
∂x2 . (3.21)
By recognizing α = R(∆x)2 = (P/∆t)(∆x)2 we can update our scaling for
∆xfrom (3.16) to suit the thermal diffusivity α for our specific heat equation
∆x =
rα∆t
P =√
2α∆t (3.22)
where we have set σ2= 2α.
For higher dimensions, the same arguments can be made to derive the general equation in D dimensions. By first looking at the two-dimensional case we have
Ni,j(t+∆t) = Ni,j(t)+R∆t [Ni+1,j(t) + Ni−1,j(t) + Ni,j+1(t) + Ni,j−1(t)− 4Ni,j(t)]
(3.23) which can be written as
Ni,j(t + ∆t)− Ni,j(t)
∆t
= h2 Ni+1,j(t)− 2Ni,j+ Ni−1,j(t)
h2 +Ni,j+1(t)− 2Ni(t) + Ni,j−1(t) h2
(3.24) where h2 = ∆x2 = ∆y2. From here we see that
∂N
∂t = Rh2 ∂2N
∂x2 +∂2N
∂y2
. (3.25)
In the two-dimensional case there are 4 possible directions for each walker, as denoted in (3.23), and we say each direction is equally probable, P = 1/4.
This gives a the new scaling h = ∆x = ∆y = √
4α∆t. Thus for the D- dimensional heat equation,
∂T
∂t = α∇2T , (3.26)
the scaling for the grid must be chosen as h = ∆x1 = ∆x2 = . . . = ∆xD =
√2Dα∆t.
To summarize this section we will look at an example. Consider the heat conduction problem for an infinite rod, with a heat impulse at t = 0
∂T
∂t = α∂∂x2T2 , x∈ R , t > 0 , T (x, 0) = δ(x) ,
T (x, t)→ 0 when t → ∞ .
(3.27)
From [8] we can see that the solution to this problem is T (x, t) = 1
√4παte−x2/(4αt). (3.28) This problem is easy to model with the RWM since it is defined for all real values and thus has no boundaries. If we approach this problem as a single
(a) N=5 time steps (b) N=100 time steps
Figure 3.3: Analytical vs. random walk solution for problem (3.27), t = 1 s.
Brownian motion starting at x0 = 0and instead use the scaling (3.22), we can in the same way as in the section above approximate this with a random walk
Wjδ ≈ Sj
√j r α
Pt . (3.29)
This random walk will have mean 0 and variance (α/P )t as j → ∞ which gives the probability density function
√P
√2παte−P x2/(2αt). (3.30)
By choosing P = 1/2 as in the one-dimensional random walk this gives the same result as for the analytical derivation (3.28). In figure 3.3 we can see a comparsion between the analytical solution and the discrete probability density function for a random walk.
3.3 A random walk model with boundary conditions
In the previous problem we adapted a RWM to a problem defined on the whole x-axis, if we instead want to treat our starting problem (1.1)-(1.3), it is necessary to take care of the boundary. This is done by discretizing the space and time according to (3.22). Here we choose the following, where we for simplicity set α = 1:
∆x = xi+1− xi = 10−2, i = 0, 1, . . . , N− 1 , x0 = 0 , xN = 1 ,
∆t = tj+1− tj = 5· 10−5, j = 0, 1, . . . , M − 1 , t0= 0 , tM = 1 . (3.31) The initial condition T (x, 0) = 1 will here be represented by one walker
starting at T (xi, t0) for all i. By the next timestep t1 all walkers will have moved one step either to the right or to the left. In our starting problem, where we have homogenous Dirichlet conditions, the walkers that reach the boundaries will be absorbed and disappear. In the case of inhomogenous Dirichlet conditions, T (0, t) 6= 0 as in the Stefan problem (2.2), we will also have walkers starting from the boundary. The number of walkers starting at T (x0, tj)will here be set according to f(t), where f(t) = 1 will be represented by one walker starting at T (x0, tj)for all j. We will iterate over the timesteps until all walkers have reached the boundaries or the maximum time tM is attained. A problem so far is that the result of our model with one walker, representing a temperature difference of 1°C per volume unit, might differ a lot depending on how each random walk turns out. In reality, the moving particles causing thermal diffusion representing that raise of temperature are a large number. Therefore, to get an accurate result we multiply the number of walkers starting at all points defined by initial- or boundary conditions with a large number n and at the end we divide the value at all points with n.
3.4 Modelling the moving boundary
To solve the Stefan problem with the RWM, the critical part is how to handle the moving boundary s(t). To make up a model for the movement of the boundary s(t) we start in the same way as in section 2.1. In (2.10) we established that the heat required to move the boundary a small step ∆s is
Q = S· ∆s · ρl , (3.32)
and thus
∆s = Q
Sρl . (3.33)
We want to compare this with the heat represented by one walker as it raises
S
∆x
Figure 3.4: One walker raises the temperature of the above volume with 1°C.
the temperature 1°C of the volume ∆x × S m3, see fig. 3.4. This could be expressed as
Qwalker= c· ρ · V · ∆T = c · ρ · ∆x · S · 1°C . (3.34) By combining (3.33) and (3.34), we see that
∆s = cρ∆xS ρlS = c
l∆x . (3.35)
So for every walker absorbed by the boundary the boundary will move the increment ∆s. To adjust for the multiplication with n at the starting points, as mentioned above, we also need to correct the step length ∆s by dividing with n. With this said, we can see that the moving boundary will have the position i in the x-grid when
i≤ sk
∆x < i + 1 , (3.36)
where we denote sk as the sum of the s-steps sk= sk−1+∆s
n . (3.37)
It might be interesting to see how the ratio between ∆s and ∆x turns out as we insert physical parameters to c and l. For water at 0°C we have the values c = 4.22 kJ/(kg·K) and l = 334 kJ/kg, which gives us
c
l ≈ 0.0126 , (3.38)
and we see that ∆s < ∆x. Note that in the opposite case, if ∆s ∆x, the boundary will move several x-steps as it is reached by one walker and this will lead to poor results when modelling the movement of the boundary.
Thus, if we choose c/l > 1 we have to compensate by increasing the number of iterations n and thereby decreasing the step size ∆s/n. So a rule of thumb to yield a good approximation of the boundary is to choose n such that
∆s
n ∆x . (3.39)
4 Results
4.1 Heat equation with fixed boundary
In (1.7) we presented an analytical solution for the heat conduction problem (1.1)-(1.3) for the initial temperature g(x) = 1°C. Fig. 4.1 shows the tem- perature distributions T (x, t) for the analytical result and the RWM solution with α = 1. In fig. 4.2 we see a comparison between the analytical result and the RWM in the cross-section at x = 0.5.
(a) Analytical solution from (1.7) with 100 Fourier terms.
(b) RWM solution with n = 104 and
∆x = 0.01.
Figure 4.1: Solutions for the heat conduction problem with initial condition g(x) = 1 for t∈ [0, 0.4].
Figure 4.2: Comparison between RWM with n = 104, ∆x = 0.01 and ana- lytical solution for t∈ [0, 0.4].
4.2 Stefan problem with constant boundary condition
In (2.8) we presented the analytical solution for the Stefan problem (2.1) - (2.6) when f(t) = T0. Fig. 4.3 shows the temperature distributions T (x, t) for the analytical result and the RWM solution with T0 = 1°C, α = 1 and β = 1. The green respective the red line denotes the solid-liquid interface.
In fig. 4.4a we compare different number of iterations n for the RWM in a cross-section at t = 0.5 s. The last figure, 4.4b compares different sizes of the step length ∆x in a plot of the moving boundary s(t).
(a) Analytical solution from (2.8). (b) RWM solution with n = 104 and
∆x = 0.01.
Figure 4.3: Solutions for Stefan problem with boundary condition f (t) = 1 for t∈ [0, 1].
(a) RWM solutions with ∆x = 0.01 and different number of iterations n, FDM and analytical solutions, x ∈ [0, 0.4].
(b) Solution for the moving boundary s(t) with analytical values and RWM for n = 104and different step lengths
∆t, t∈ [0, 0.3].
Figure 4.4
4.3 Stefan problem with boundary condition f (t) = et− 1 In (2.9) we presented the analytical solution for the Stefan problem (2.1) - (2.6) when f(t) = et−1. Fig. 4.5 shows the temperature distributions T (x, t) for the analytical result and the RWM solution with α = 1 and β = 1. The green respective the red line denotes the solid-liquid interface. In fig. 4.6 we see a comparison between the analytical result, the RWM and FDM in the cross-section at t = 0.5 s.
(a) Analytical solution from (2.9). (b) RWM solution with n = 104, ∆x = 0.01.
Figure 4.5: Solutions for Stefan problem with boundary condition f (t) = et− 1 for t ∈ [0, 1].
Figure 4.6: Comparison between RWM with n = 104, ∆x = 0.01, FDM and analytical solutions for x∈ [0, 0.4].
4.4 Stefan problem with oscillating boundary condition In the introduction we discussed that it would be interesting to evaluate an oscillating function in the fixed boundary, and since we want f(t) ≥ 0 ∀ t we choose f(t) = 1 + sin (3πt). Then the RWM solution for the Stefan problem (2.1) - (2.6) yields the temperature distribution T (x, t) as seen in fig. 4.7.
Here we set α = 1 and β = 2. The red line in fig. 4.7 denotes the solid-liquid interface. Due to limitations in the implementation of the finite difference method, we are restricted to consider f(t) = sin t on the boundary when comparing the two methods. The result for the latter case is shown in the cross-section at t = 1 s in fig. 4.8.
Figure 4.7: RWM solution for f (t) = 1 + sin (3πt) with n = 104, ∆x = 0.01 on the interval t∈ [0, 1].
Figure 4.8: Comparison between RWM with n = 104, ∆x = 0.01 and FDM for f (t) = sin t at t = 1 on the interval x∈ [0, 0.4].
4.5 Stefan problem with boundary condition according to daytime temperature variations
To apply our RWM model in a simulation of melting ice we here assume the phyiscal constants for water, α ≈ 0.1429 mm2/s and β ≈ 0.0126 K−1. We want to model the melting of ice according to daytime temperature variations and therefore we set the values at the fixed boundary to observed air temperatures from Örebro airport 1-3 March 2019 [15]. Assuming the air temperature at the ice surface is a simplification and does not deal with the temperature gradient between air and ice or heat transport by convection or radiation. Nevertheless, fig. 4.9 and 4.10 gives a qualitative view of the evolution of the melted water.
Figure 4.9: RWM model where f (t) is set to observed day temperatures at Örebro airport 1-3 March 2019. x in mm and t∈ [0, 62 h].
Figure 4.10: Another view of fig. 4.9.
5 Discussion
To evaluate the results from the previous section, we can see qualitatively from fig. 4.4a that the RWM solution converges to the analytical as n → ∞.
As we proved in section 3.1, we also have a convergence criteria ∆t → 0, which is confirmed by fig. 4.4b. Due to limitations in the project time, we have left the work with a more thoroughly error analysis for the future.
For a similar error analysis on the RWM, we direct to a work by Lockby [16].
There are a few simplifications in our model for the Stefan problem that can be improved in a more detailed study. Among the physical simplifications we have mentioned our assumption that we use the same density for water and ice, ρL = ρS which is not the real case. To adjust this property we would have to make a different approach to the Stefan condition. We may also want to consider a temperature distribution in the ice phase Tice(x, t) 6= 0, which would lead us to a two phase problem with a system of PDE:s. Some formulations of two phase problems have analytical solutions, see e.g. [17], to compare with for future studies.
Regarding the RWM, this work has proved that in comparison with the finite difference method from Jonsson [3], the RWM has required considerably less programming effort, measured in lines of code (see appendix). The programming code for the RWM is also flexible and compared to the finite difference method, it requires small effort to change the boundary- and initial conditions for the heat equation arbitrarily.
5.1 Applications
As mentioned in the introduction, there are several applications for the Ste- fan problem. Only by looking at the original purpose of Stefan’s article in 1891, which was to model the arctic ices, we see that this is highly relevant today due to the demand of accurate climate models. According to Hunke et al. [18], Stefan’s one dimensional thermodynamical model is still in use for global climate models, although the complete thermodynamical sea-ice models are more complex. Hunke et al. also states that "Fresh numerical approaches and algorithm improvements play a critical role in the develop- ment process, as climate models continue to push the limits of computational power" [18]. Hence thermodynamical ice-sea models might be a subject for future work with the RWM approach.
Other areas where we might want to model a solid-liquid interface is in 3D-printing, freezing of food, solidifying of building components, lithium-ion batteries and in oil extraction using steam assisted gravity drainage. For the latter, Chen et al. [19] have solved an inverse Stefan problem in order to
model the velocity of the expanding steam chamber used to decrease the vis- cosity of bitumen. Regarding lithium-ion batteries, the diffusion of lithium ions in the battery is separated into two phases, one where lithium ions are evenly distributed and one where they are not present. To be able to better compute the properties of the battery, such as life-time and capacity, it is helpful to estimate the movement of the interface between these two phases, which can be interpreted as a Stefan problem [20].
6 Conclusion
In accordance with our opening objectives we have succesfully numerically evaluated a solution for the one-dimensional Stefan problem with a general time-dependent boundary condition at the fixed boundary using a RWM. A qualitative error analysis gives that the error can be reduced by increasing the iterations, n → ∞, and by decreasing the timestep, ∆t → 0. In com- parison with the finite difference method, the RWM has proven to be easy to implement and more general in sense of effort needed to substitute the boundary- and initial conditions.
References
[1] B.V Budaev and D.B. Bogy. “Probabilistic approach to the Lamé equa- tions of linear elastostatics”. In: International Journal of Solids and Structures 40 (2003), pp. 6285–6306.
[2] M.K. Chati et al. “Random walk method for the two- and three- dimensional Laplace, Poisson and Helmholtz’s equations”. In: Inter- national Journal for Numerical Methods in Engineering 51 (2001), pp. 1133–1156.
[3] T. Jonsson. “On the one dimensional Stefan problem with some nu- merical analysis”. Umeå universitet, 2013. url: http://www.diva- portal.se/smash/record.jsf?pid=diva2:647481.
[4] N. Masuda et al. “Random walks and diffusion on networks”. In: Physics Reports 716-717 (2017), pp. 1–58.
[5] A. Einstein. “Über die von der molekularkinetischen Theorie der Wärme geforderte Bewegung von in ruhenden Flüssigkeiten suspendierten Teilchen”.
In: Ann. Physik 322 (1905), pp. 549–560.
[6] J. Stefan. “Über die Theorie der Eisbildung, insbesondere über die Eisbildung im Polarmeere”. In: Ann. Physik Chemie 42 (1891), pp. 269–
286.
[7] L.N. Tao. “The Stefan problem with arbitrary initial and boundary conditions”. In: Quarterly of applied mathematics October (1978), pp. 223–
233.
[8] G. Sparr and A. Sparr. Kontinuerliga system. Studentlitteratur, 2000.
isbn: 9144013558.
[9] S.L. Mitchell and M. Vynnycky. “Finite-difference methods with in- creased accuracy and correct initialization for one-dimensional Ste- fan problems”. In: Applied Mathematics and Computation 215 (2009), pp. 1609–1621.
[10] D. Andreucci. Lecture notes on the Stefan problem. Università di Roma La Sapienza. 2004. url: http://www.dmmm.uniroma1.it/pubblicazioni/
doc/phd_quaderni/02-1-and.pdf. (accessed: 26.05.2019).
[11] R. Feynman. The Feynman lectures on physics. Chapter 6: Probabil- ity. url: http://www.feynmanlectures.caltech.edu/I_06.html.
(accessed: 08.05.2019).
[12] K. Sigman. IEOR 4700: Notes on Brownian Motion. 2006. url: http:
/ / www . columbia . edu / ~ks20 / FE - Notes / 4700 - 07 - Notes - BM . pdf. (accessed: 28.05.2019).
[13] G. F. Lawler. Random Walk and the Heat Equation. University of Chicago, Department of Mathematics, 2010. isbn: 9780821884713.
[14] D. K. Dysthe. Chapter 4: Introduction to diffusion and random walks.
url: https://www.uio.no/studier/emner/matnat/fys/FYS2160/
h17/simuleringsopgaver/virrevandrer_diffusjon.pdf. (accessed:
19.05.2019).
[15] SMHI. Mina observationer - WOW. Örebro flygplats 2019-03-01 - 2019- 03-03. 2019. url: https : / / www . smhi . se / vadret / vadret - i - sverige/mina-observationer-wow#id=95130. (accessed: 31.05.2019).
[16] A. Lockby. “En slumpvandringsmetod för värmeledningsekvationen med rörlig rand”. Örebro universitet, 2016. url: http://www.diva-portal.
se/smash/record.jsf?pid=diva2:935914.
[17] R.M. Furzeland. “A Comparative Study of Numerical Methods for Moving Boundary Problems”. In: J. Inst. Maths Applies 26 (1980), pp. 411–429.
[18] E.C. Hunke et al. “Sea-ice models for climate study: retrospective and new directions”. In: Journal of Glaciology 56.200 (2010), pp. 1162–
1172.
[19] X. Chen et al. “The Application of Stefan Problem in Calculating the Lateral Movement of Steam Chamber in SAGD”. In: Mathematical Problems in Engineering (2015). doi: http://dx.doi.org/10.1155/
2015/372581.
[20] T. Hoffman et al. “Numerical simulation of phase separation in cathode materials of lithium ion batteries”. In: International Journal of Solids and Structures 100-101 (2016), pp. 456–469.
