Monolithic simulation of convection-coupled
phase-change - verification and reproducibility

arXiv:1801.03429v1 [physics.flu-dyn] 10 Jan 2018

Alexander G. Zimmerman and Julia Kowalski
AICES Graduate School, RWTH Aachen University,
Schinkelstr. 2, 52062 Aachen, Germany
zimmerman@aices.rwth-aachen.de

Abstract. Phase interfaces in melting and solidification processes are
strongly affected by the presence of convection in the liquid. One way
of modeling their transient evolution is to couple an incompressible flow
model to an energy balance in enthalpy formulation. Two strong nonlinearities arise, which account for the viscosity variation between phases
and the latent heat of fusion at the phase interface.
The resulting coupled system of PDE’s can be solved by a single-domain
semi-phase-field, variable viscosity, finite element method with monolithic system coupling and global Newton linearization [9]. A robust computational model for realistic phase-change regimes furthermore requires
a flexible implementation based on sophisticated mesh adaptivity. In this
article, we present first steps towards implementing such a computational
model into a simulation tool which we call Phaseflow [1].
Phaseflow utilizes the finite element software FEniCS [3], which includes
a dual-weighted residual method for goal-oriented adaptive mesh refinement. Phaseflow is an open-source, dimension-independent implementation that, upon an appropriate parameter choice, reduces to classical benchmark situations including the lid-driven cavity and the Stefan
problem. We present and discuss numerical results for these, an octadecane PCM convection-coupled melting benchmark, and a preliminary
3D convection-coupled melting example, demonstrating the flexible implementation. Though being preliminary, the latter is, to our knowledge,
the first published 3D result for this method. In our work, we especially
emphasize reproducibility and provide an easy-to-use portable software
container using Docker [7].

Keywords: incompressible flow, finite element method, Newton method,
phase-change, reproducibility

2

A. G. Zimmerman, J. Kowalski

Nomenclature
t
tf
x
x∗
p
u
T
fB
r
Tr
Pr
Ra
Ste
µ
φ
∆t
()n

Time
Final time

Spatial coordinates x = x y z
phase-change interface’s position
Pressure field
Velocity field
Temperature field
Buoyancy force
Regularization smoothing factor
Regularization central temperature
Prandtl number
Rayleigh number
Stefan number
Dynamic viscosity
Semi-phase-field
Time step size
Values from discrete time n

1

Introduction

γ Coefficient for penalty stabilization
h Finite element cell diameter
F The vector-valued strong form
F Functional for the variational form
Ω The spatial domain
∂Ω Boundary of the spatial domain
∂ΩT Dirichlet boundary for T
V The scalar solution function space
V The vector solution function space
W Mixed finite element function
 space
w System solution w = p u T
ψ Finite element basis functions
M Adaptive goal functional
M Adaptive solver tolerance
δw Residual of linearized system
ω Newton method relaxation factor
()k Values from Newton iteration k

The melting and solidification of so-called phase-change materials (PCM’s) are
relevant to many applications ranging from the design of latent heat based energy
storage devices [12], to ice-ocean coupling and its effects on Earth’s climate
[10], to the evolution of ocean worlds on the icy moons of our solar system
[14] and the design of robotic melting probes for their exploration [16]. The
predictive modeling of phase-change systems is, however, challenging due to 1.
strong nonlinearities at the phase-change interface (PCI), 2. the coupling of
several physical processes (i.e multi-physics), and 3. a large range of relevant
scales both in space and time (i.e. multi-scale). Any mathematical model of
a complex phase-change process hence manifests as a multi-scale and multiparameter, nonlinear PDE system. We aim to develop a robust and flexible
model to simulate these systems.
In this work, we will focus on phase-change in the presence of liquid convection. Convection can have a tremendous effect on the evolution of phaseinterfaces, as shown in Figure 1. A comprehensive introduction to melting and
freezing without convection is given in [2]. A mathematical model that accounts
for convection-coupled phase-change is presented in [5]. Therein mushy layer theory is introduced, which provides a model for understanding the PCI at macroscale. The physical system is mathematically modeled by considering balance
laws for mass, momentum, and energy, with the momentum and energy balances coupled via buoyancy (which forces natural convection).
This constitutes a difficult multi-physics problem. In a PCM domain, certain physics dominate in the solid and liquid subdomains. To solve the coupled

3

Hot plate

Monolithic simulation of convection-coupled phase-change

𝒈
Ice

Fig. 1. Sketch of the physical model superimposed on a melting experiment [17].
Heating from the left advances the phase-change interface (PCI) rightward. Buoyancy
forces a convection cell, causing the upper PCI to advance more rapidly. In the absence
of convection, the PCI would propagate as a planar front, indicating the significant
impact of convection on the PCI’s evolution.

problem, out of many approaches in the literature, we employ the single-domain
semi-phase-field enthalpy method. With this approach, the convection itself is
handled by an incompressible Navier-Stokes flow model, while the energy balance is modeled as the convection and diffusion of an enthalpy field. The key
idea is to solve the same equations on the entire domain. This is commonly referred to as the fixed grid approach, for which various techniques are reviewed
in [19]. Using an enthalpy method, phase-change latent heat effects are isolated
to a source term in the energy balance [18]. Single-domain methods require a
velocity correction scheme, for which we select the variable viscosity method,
thereby treating the solid as a highly viscous fluid.
Efficiently applying this approach requires local mesh refinement to resolve
the moving PCI. This can be accomplished either with mesh refinement or fronttracking [12]. Adopting the former approach, we leverage the existing work on
goal-oriented adaptive mesh refinement (AMR) methods, particularly the dualweighted residual method [4]. These days, goal-oriented AMR is widely practiced,
and multiple open-source software libraries provide this capability.
In the following, Section 2 describes the mathematical model. Section 3
presents numerical methods for the model’s discretization and linearization. Section 4 presents our implementation, Phaseflow [1], based on the finite element
software FEniCS [3]. Finally, Section 5 presents verification via comparison to
benchmark problems, a convergence study for the 1D Stefan problem, and a preliminary result for convection-coupled melting in a 3D domain. We close with
conclusions and an outlook.

4

2
2.1

A. G. Zimmerman, J. Kowalski

Mathematical Phase-Change Model
The Governing Equations

For the coupled phase-change system, out of many approaches from the literature, we adopt an enthalpy formulated [18], single-domain semi-phase-field [6],
variable viscosity model. The mass and momentum balances (1) and (2) are given
by the incompressible Navier-Stokes equations, with velocity field u = u(x, t)
and pressure field p = p(x, T ). Invoking the Boussinesq approximation, we extend the momentum equation with a temperature-dependent buoyancy forcing
term fB (T ) which couples it to the energy equation. This approach is well established in the context of natural convection [20]. Furthermore, we consider the
phase-dependent viscosity µ(φ) in the momentum equation, where φ = φ(T ) is
the temperature-dependent phase. For constant heat capacity, the energy balance
in enthalpy form reduces to (3), which is an extended form of the convectiondiffusion equation for the temperature field T = T (x, t). The diffusion term in (3)
1 ∂
is scaled by the Prandtl number Pr. The nonlinear source term Ste
∂t φ accounts
for the phase-change latent heat, where Ste is the Stefan number. Altogether the
system of governing equations is
∇·u=0

(1)

∂
u + (u · ∇) u + ∇p − ∇ · (2µ(φ)D(u)) + fB (T ) = 0
(2)
∂t
∂
1 ∂
1
T−
φ + ∇ · (T u) −
∆T = 0
(3)
∂t
Ste ∂t
Pr


T
where the symmetric part of the rate-of-strain tensor is D(u) = 21 ∇u + (∇u) .
The equations are unitless per the normalization in [9], shifted always such that
T = 0 corresponds to the temperature of fusion.
We refer to φ as a semi -phase-field [6], because we do not treat φ as an
additional unknown. Rather, φ maps the temperature field to values between
zero and one. For this we use


1
Tr − T
φ(T ) =
1 + tanh
,
(4)
2
r
where Tr is the central temperature and r is a smoothing parameter. The Newton
method requires the differentiability of (4). Subject to this requirement, there
are many other useful regularizations from which to choose. Figure 2 plots φ and
φ0 . For the Stefan problem in Section 5.3, we simply set Tr equal to the physical
temperature of fusion (i.e. Tr = 0 for the normalized system). Following [9], for
the coupled convection problem, we set Tr such that the strongest variation is
localized near the newly appearing phase, e.g. liquid for melting. Physically, the
region where 0 < φ < 1 is analogous to a mushy layer. In this sense, we can
view φ as the solid volume-fraction. Given φ, we define the phase-dependent
viscosity as functions of the liquid and solid values (assuming constants µL and

Monolithic simulation of convection-coupled phase-change

5

1
0.8
0.6
0.4
0.2
0
-0.2
-0.4
-0.6
-0.8
-1
-0.1

-0.05

0

0.05

0.1

Fig. 2. The regularized semi-phase-field (4) with Tr = 0 and two values of r. For
smaller r, φ(T ) steepens and its derivative φ0 (T ) approaches the Dirac delta function.
Here the derivative is scaled by 2r for convenience.

µS respectively) with
µ(φ) = µL + (µS − µL )φ,

(5)

With the variable viscosity method, the solid is treated as a highly viscous
fluid, forcing its flow velocity to zero. In general, not only viscosity varies in
the PCM domain. For example, the thermal material properties of water-ice
vary significantly between the solid and liquid phases. In the case of octadecane
PCM’s, such as the benchmark in Section 5.4, we may treat these material
properties as constants [20].
2.2

The Initial Boundary Value Problem

To simulate the time evolution of the PCM system, we will solve (8) as an initial
boundary value problem subject to the initial values
T
T
p(x, t), u(x, t), T (x, t) = 0, u0 (x), T0 (x)

∀(x, t) ∈ Ω × 0

(6)

and boundary conditions
u(x, t) = uD (x, t) ∀(x, t) ∈ ∂Ω × (0, tf ],
T (x, t) = TD (x, t) ∀x ∈ ∂ΩT × (0, tf ],

(7)

(n̂ · ∇)T (x, t) = 0 ∀x ∈ (∂Ω \ ∂ΩT ) × (0, tf ]
For the applications in Section 5, we apply for all boundaries the velocity no-slip
condition. This implies homogeneous Dirichlet boundary conditions, uD = 0, for

6

A. G. Zimmerman, J. Kowalski

all cases except the lid-driven cavity in Section 5.1, where the ”moving” lid makes
uD non-homogeneous. The notation ∂ΩT refers to the boundary subdomain with
non-homogeneous Dirichlet boundary conditions on the temperature. Physically
this models a thermal reservoir beyond the wall, keeping the wall at a constant
temperature. All other walls are prescribed homogeneous Neumann boundary
conditions on the temperature, which physically models thermally insulated (i.e.
adiabatic) walls.
The generality of this initial boundary value problem allows us to solve a
variety of interesting benchmarks by setting the appropriate parameters. In Section 5 we will demonstrate the 2D lid-driven cavity, 2D heat-driven cavity, 1D
Stefan problem, 2D convection-coupled melting of an octadecane PCM, and 3D
convection-coupled melting. Our implementation which we present in Section
4 is as versatile as the model. Testing these benchmarks is accomplished with
short Pythons scripts specifying the parameters.

3

Numerical Methods

We base our work on the numerical approach from [9]. Therefore, we discretize
in time via finite differences, discretize in space via the finite element method
(FEM) with a penalty formulation for stabilization, couple the system monolithically, and solve the nonlinear system globally via Newton’s method.

3.1

The Nonlinear Variational Problem


We denote the system’s solution as w = p u T and write the strong form (1),
(2), and (3) as the vector-valued functional
F(w) = 0 ∀(x, t) ∈ Ω × (0, tf ],

(8)

with the three components respectively being the mass, momentum, and energy
equations. To stabilize the finite element method, following the penalty formulation in [9], we add a pressure stabilizing term γp to the mass component. The
coefficient γ is generally small, and in this case is γ = 10−7 . We consider 1D, 2D,
and 3D spatial domains, respectively allowing for Ω ⊂ R1 , Ω ⊂ R2 , or Ω ⊂ R3 .
Per the standard Galerkin finite element method, we write the variational (i.e.
weak) problem whose solution approximates (8). To obtain the monolithic system, we use mixed finite elements
 [8]. Therefore, we multiply (8) from the left by
test functions ψ = ψp ψ u ψT and integrate by parts over the domain. Finally,
we employ the fully implicit Euler method by substituting the time discretizations ∂t u = ∆t−1 (un+1 − un ), ∂t T = ∆t−1 (Tn+1 − Tn ), ∂t φ = ∆t−1 (φn+1 − φn ).

Monolithic simulation of convection-coupled phase-change

7

Altogether, for homogeneous Neumann boundary conditions, this yields the
time-discrete nonlinear variational form
Z

ψp ψ u ψT F(p, u, T )dx =
F (ψ, w) =
Ω

b(u, ψp ) − (ψp , γp)

1
+ ψu ,
(un+1 − un ) + fB (T ) + c(u; u, ψ u ) + b(ψ u , p) + a(µ(φ); u, ψ u )
∆t

 

1
1
1
+
ψT , Tn+1 − Tn −
(φ(Tn+1 ) − φ(Tn )) + ∇ψT , ∇T − T u
∆t
Ste
Pr
(9)
R
R
where we use the short-hand (u, v) = Ω uvdx, (u, v) = Ω u · vdx for integrating inner products. Additionally, we use a common notation [11] for the linear,
bilinear, and trilinear forms of the variational Navier-Stokes equations
Z
a : V × V → R, a(µ; u, v) = 2
µD(u) : D(v)dx
Z Ω
b : V × V → R, b(u, p) = −
p∇ · udx
(10)
Ω
Z
c : V × V × V → R, c(u; z, v) =
vT (∇z) udx


Ω

Given (9), we write the variational problem as
Find w ∈ W such that
F (ψ, w) = 0

∀ψ ∈ Ŵ

(11)

where W = V ×V×V and Ŵ = V̂ ×V̂×V̂ . This distinction comes from how nonhomogeneous boundary conditions are often handled in finite element method
implementations. The solution is split into homogeneous and non-homogeneous
parts, the former is found, and then the
 latter is reconstructed. Therefore, for
example, T belongs to the space V = v ∈ H01 (Ω) : v = TD on ∂Ω , while ψT

belongs to the space V̂ = v ∈ H01 (Ω) : v = 0 on ∂Ω , where H01 (Ω) is the classical Hilbert space fulfilling requirements for continuity and compact support.
The vector-valued function space V depends on the spatial domain’s dimensionality, e.g. V = V × V in 2D. For the incompressible Navier-Stokes solution,
we use the Taylor-Hood element [11], i.e. we restrict the pressure solution to
piece-wise continuous linear Lagrange polynomials, and we restrict the velocity
solution to piece-wise continuous quadratic Lagrange polynomials. We use the
same polynomial space for temperature and pressure.
3.2

Linearization via Newton’s Method

We apply Newton’s method by solving a sequence (indexed by superscript k) of
linear problems

8

A. G. Zimmerman, J. Kowalski

Find δw ∈ W such that
Dw F (ψ, wk ; δw) = F (ψ, wk ) ∀ψ Ŵ

(12)

for the residual δw which updates the solution wk , converging wk to the approximate solution of the nonlinear problem. The Gâteaux derivative of F (ψ, wk ),
d
F (ψ, wk + δw)|=0 , is given by
defined as Dw F (ψ, wk ; δw) ≡ d
Dw F (ψ, wk ; δw) = b(δu, ψp ) − (ψp , γδp)

1
0
k
+ ψu ,
δu + δT fB (T ) + c(uk ; δu, ψ u ) + c(δu; uk , ψ u )
∆t


+b(ψ u , δp) + a δT µ0 (T k ); uk , ψ u + a µ(T k ); δu, ψ u

 

1
1 0 k
1
k
k
+
ψT δT, 1 −
φ (T ) + ∇ψT , ∇δT − T δu − δT u ,
∆t
Ste
Pr


(13)

For a given wk (13) and hence (12) are linear with respect to the unknown
Newton residual δw. In typical fashion, we define the test functions, solution and
Newton residual as linear combinations of the same basis. Upon selecting a mesh
and concrete basis, this allows (12) to be re-written as a discrete linear system
of the form Ax = b, which can be efficiently solved on a computer by standard
methods and software. In this work, we directly solve each linear system with
LU decomposition (using the interface of FEniCS, discussed in Section 4).
We use the latest discrete time solution as the initial guess for the Newton
solver, i.e. we initialize w0 = wn . After each iteration solving for δw, we update
the solution with the relaxed residual wk+1 := wk + ωδw, where 0 < ω ≤ 1 is a
relaxation factor. Relaxing Newton’s method is often useful for highly nonlinear
problems. The 1D and 2D results in Section 5 use the full Newton (ω = 1)
method. So far only the preliminary 3D result requires relaxation, where we will
set ω = 0.8.
3.3

Adaptive Mesh Refinement

The single-domain approach requires local mesh refinement, or else the computational cost would quickly become impractical, especially in 3D. Furthermore,
we cannot a priori prescribe where to locally refine the mesh, since our goal is to
predict the position of the PCI. This means we must employ adaptive mesh refinement, now commonly referred to as AMR. The theory of AMR fundamentally
requires an error estimator, which exist for many interesting problems [15]. In
the context of phase-change simulations, hierarchical error estimators have been
derived for the Stefan problem [6]; but no such rigorous work has been completed
for the problem with coupled convection. Promising results are reported in [9],
wherein they used a mesh adaptivity procedure by metric control which was
particular to the software library and its Delaunay mesh generation procedure.
Unfortunately, it is unclear exactly which metrics were used for adaptivity in
those results.

Monolithic simulation of convection-coupled phase-change

9

We instead employ goal-oriented AMR [4]. To discuss AMR, let us briefly
write the spatially discrete problem which is dependent on a mesh Wh ⊂ W(Ω):
Find wh ∈ Wh ⊂ W(Ω) such that
F (ψ h , wh ) = 0 ∀ψ h ∈ Ŵh ⊂ Ŵ

(14)

Goal-oriented AMR requires some goal functional M (w) to be integrated over
the domain Ω. The goal-oriented adaptive solution of (14) can then be written:
Find Wh ⊂ W(Ω) and wh ∈ Wh such that
|M (w) − M (wh )| < M

(15)

where M is some prescribed tolerance. Since M (w) is unknown, we still require
an error estimator. To this end, we use the dual-weighted residual method, as
implemented in FEniCS [3]. Computing cell-wise error estimates requires solving
a linearized adjoint (with respect to the goal) problem. The primal and adjoint
problems are solved on a hierarchy of meshes. Computing the linearized adjoint
solution on each mesh is relatively cheap compared to solving the nonlinear
primal problem. See [4] for a full explanation of the method. For the adaptive
solutions in this work, we set the goal functional
Z
M (w) =
φ(T )dx
(16)
Ω

which represents the volume of solid material remaining in the domain.

4

Implementation

The method presented in this paper, along with its application to a series of test
problems, were implemented by the author into an open-source Python module
named Phaseflow [1], using the finite element library FEniCS [3].
4.1

FEniCS

The abstract interface of FEniCS [3] makes it an ideal library to quickly prototype models and methods which use FEM. Additionally, standard algorithms
for Newton linearization and goal-oriented AMR are already implemented. FEniCS is an umbrella project. A major back-end component is the C++ library
DOLFIN, which implements most of the classes and methods seen by the user.
Python interfaces are generated mostly automatically, and indeed FEniCS is
primarily used as a Python module. From a conceptual perspective, we are most
interested in two components of FEniCS: the Unified Form Lanugage (UFL) and
the FEniCS Form Compiler (FFC). UFL allows us to write down the abstract
discrete variational (i.e. weak) form in a way that is understood by FEniCS.
This means that to implement variational forms, we write them as source code
almost character-for-character. The FEniCS Form Compiler (FFC) then automatically implements the FEM matrix assembly routine in optimized C++ code
using just-in-time (JIT) compiling.

10

4.2

A. G. Zimmerman, J. Kowalski

Phaseflow

Phaseflow [1], created by the author, is a Python module maintained openly on
GitHub. Phaseflow implements the methods of this paper using the open-source
finite element library FEniCS [3]. The interface allows users to run these methods for the variational form (9) with any set of similarity parameters Ste, Pr,
buoyancy model fB (T ), regularization φ(T ), liquid and solid viscosities µL and
µS , time step size ∆t, and stabilization coefficient γ. The initial values (6) are
set with a function in the solution space. Combined with the built-in FEniCS
method, this allows the user to interpolate general mathematical expressions
written in the C syntax, or to use an existing solution. Dirichlet boundary conditions (7) are also set with general C mathematical expressions. Additionally
the interface allows users to control algorithm parameters, such as tolerances for
the solvers, the Newton relaxation parameter ω, and other options.
The interface accepts any FEniCS mesh object, which can either be generated from scratch or can be converted from a standard format via the FEniCS
library. Leveraging the abstract interfaces of FEniCS, Phaseflow is dimensionindependent, and has been applied to 1D, 2D, and 3D spatial domains. In this
work, unit square and rectangular prism domains have been used. Phaseflow
writes all solutions to disk in the efficient XDMF+HDF5 format using the builtin FEniCS interface. Furthermore, this is combined with the H5Py library to
write checkpoint files for easily restarting simulations at later times. Phaseflow
has other ease-of-use features, such as a stopping criterion for steady state solutions. This was useful for the heat-driven cavity test in Section 5.2, where the
Newton method required a small time step size.
4.3

Reproducibility with Docker

Phaseflow leverages open technology, and it is meant to contribute further to
open scientific and engineering research. Open research in the computational
sciences is often hindered by the difficulty of reproducing results, primarily because of the ever-increasing complexity of computing environments [7]. Recently,
the software Docker has emerged as a technology which greatly facilitates reproducibility in the computational sciences.
Built upon an official FEniCS Docker image, we provide a Phaseflow Docker
with all dependencies. Enabled by the Docker image and Phaseflow’s test suite,
changes to the master branch on GitHub are continuously tested using the
Travis-CI continuous integration service [7]. Docker has been central to the
development of Phaseflow. Initially, a primary motivation for its use was the
existence of Docker images where most of the dependencies of Phaseflow are
already installed and maintained by the FEniCS developers. This allows us to
spend more time implementing models, and less time building tool chains. From
there, maintaining a Docker image which runs Phaseflow takes little effort. With
the Phaseflow Docker image in hand, we were then able to set up a continuous
integration process with remarkably little effort.

Monolithic simulation of convection-coupled phase-change

11

The vast majority of Phaseflow’s source code is covered with a suite of unit
tests and integration tests written with PyTest. Most results in this paper are
included in Phaseflow’s test suite, while the others are covered by the example
scripts and notebooks in its repository [1]. Such results are easy to reproduce,
and we highly encourage the reader to do this. To further aid reproducibility, a
specific version of Phaseflow has been archived (see the release versions at [1])
to coincide with this publication.

5

Verification

To verify our implementation, we consider a series of examples, including the liddriven cavity benchmark from [13] extended with a solid subdomain, the heatdriven cavity benchmark from [20], an approximation of the analytical 1D Stefan
problem from [2] with and without AMR, a preliminary 2D convection-coupled
melting demonstration with AMR, and a preliminary 3D convection-coupled
melting demonstration without AMR. Furthermore, for the Stefan problem, we
verify the convergence orders of the temporal and spatial discretizations, and of
the Newton method. Phaseflow’s flexible interface described in Section 4 allows
us to implement each of these applications with minimal effort, and to test the
same lines of source code, increasing our confidence with the implementation.
5.1

Lid-Driven Cavity with Solid Subdomain

To test phase-dependent viscosity, we consider an extension of the lid-driven
cavity benchmark based on data published in [13]. The standard benchmark uses
a unit square geometry. We extend the geometry below the bottom wall, and
set the temperature in the new region below the freezing temperature, such that
the liquid subdomain per (4) covers the original unit square geometry. Therefore
the initial temperature values are
(
Tc for y ≤ 0,
(17)
T0 (x) =
Th otherwise
where we chose Th = 1 and Tc = −1. For the semi-phase-field (4) we set central
temperature Tr = −0.01 and smoothing parameter r = 0.01. We set variable
viscosity (5) with µL = 1 and µS = 108 . To capture this strong variation, we
refined all cells of the initial mesh which touched coordinate y = 0 for four
refinement cycles. Note that in this case, the local refinement is not adaptive.
This problem at steady state fits our general model (9) with null buoyancy fB = 0
and with an arbitrarily large Prandtl number to nullify thermal conduction.
As boundary conditions uD , we set positive horizontal velocity on the lid, zero
horizontal velocity away from the lid, and zero vertical velocity everywhere. We
set initial velocity u0 similarly, which serves as a suitable initialization for solving
until steady state with Newton’s method. Because Phaseflow implements only
the unsteady problem, we solve a single large time step to obtain the approximate

12

A. G. Zimmerman, J. Kowalski

1

0.8

0.6

0.4

0.2

0

-0.2

-0.4
-0.2

0

0.2

0.4

0.6

0.8

1

Fig. 3. Left: Result from lid-driven cavity test (Re = 100). The mesh (which is locally
refined near the PCI) is shown in translucent gray. Velocity streamlines are shown in
black. As expected, the moving lid causes circulation in the cavity, and we see the
expected recirculation zone in the bottom-right corner of the liquid subdomain.
Right: Horizontal velocity sampled from the vertical centerline, compared to benchmark
data published in [13].

steady state solution. It is advantageous to model the steady problem this way,
when the goal is to verify components of our unsteady implementation. Figure 3
shows and discusses the result. The solution agrees very well with the published
benchmark data.
5.2

Heat-Driven Cavity

To verify the coupled energy equation without phase-change, we compare to
the heat-driven cavity benchmark data published in [20]. This problem fits our
general model with constant viscosity µS = µL = 1 and arbitrarily large Stefan
number. We handle buoyancy with an idealized linear 1 Boussinesq model
 
Ra
0
T
,
(18)
fB (T ) =
−1
PrRe2
The momentum equation (2) is scaled such that the Reynolds number Re is
always unity. For this benchmark, the Rayleigh number is Ra = 106 and the
Prandtl number is Pr = 0.71. This Rayleigh number is considered to be high,
and demonstrates substantial natural convection. We set homogeneous Dirichlet
(i.e. no slip) boundary conditions on the velocity uD = 0, and non-homogeneous
Dirichlet boundary conditions for the temperature with hot and cool temperatures, Th = 0.5 and Tc = −0.5, respectively on the left and right vertical walls,
1

The method handles nonlinear bouyancy with (13). Phaseflow’s test suite includes
a benchmark with the nonlinear density anomaly of water.

Monolithic simulation of convection-coupled phase-change

13

1
Wang 2010
Phaseflow

0.8

0.6

0.4

0.2

0
-0.1

-0.05

0

0.05

0.1

Fig. 4. Left: Result from heat-driven cavity test (Ra = 106 ). The mesh is shown in
translucent gray. Velocity streamlines are shown in black. As expected, we see the flow
circulating, being forced upward near the hot wall and downward near the cold wall.
Right: The horizontal component of velocity is plotted along the vertical centerline.
The solution agrees very well with the data published in [20].

i.e.
(
Th
TD (x) =
Tc

for x = 0,
for x = 1

(19)

Again we solve the unsteady problem until it reaches steady state. For this we
initialize the velocity field to zero (i.e. u0 = 0) and the temperature field to vary
linearly between the hot and cold walls (i.e. T0 (x) = Th +x(Tc −Th )). Unlike with
the lid-driven cavity, here we cannot obtain the steady state solution in a single
time step, because the initial guess is not sufficient for the Newton method to
converge. Therefore we solve a sequence of time steps using ∆t = 0.001. Figure
4 shows the successful result with further discussion.
5.3

The Stefan Problem

To verify the energy equation with phase-change, we compare to the analytical 1D Stefan problem as written in [2], with parameters comparable to the
octadecane PCM melting benchmark in Section 5.4, including the Stefan number Ste = 0.045. This problem fits our general model with nullified buoyancy
fB = 0, unity Prandtl number, and zero initial velocity u0 = 0 (which remains
zero). Solving the Stefan problem as a special case of the coupled problem results
in an unnecessarily large system. In this case, the vast majority of degrees of
freedom are trivial. From an implementation perspective, this exercise is quite
valuable, because we test the exact lines of code which are used for the coupled
problem.

14

A. G. Zimmerman, J. Kowalski

We again set temperature boundary conditions (19), this time with Th = 1
and Tc = −0.01. The initial temperature field is set such that a thin layer of melt
exists near the hot wall, with the rest of the domain at the cold wall temperature,
i.e.
(
Th for x ≤ x∗0 ,
T0 (x) =
(20)
Tc otherwise
We parameterize this initial PCI position as
x∗0 =

L 1−q
2
N0

(21)

where q is the number of initial hot wall refinement cycles. This ensures that
for the given initial mesh, the thinnest possible layer of melt exists. For the
semi-phase-field regularization (4) we set parameters Tr = 0, r = 0.01.
We simulate until time tf = 0.1, using the time step size ∆t = 0.001. In
testing, this size was needed to bound the point-wise error between the numerical and analytical solutions from above by T (x) − Texact (x) < 0.01 ∀x ∈
{0, 0.025, 0.05, 0.075, 0.1, 0.5, 1} at time t = tf (shown in Figure 5). This discrete
set of verification points is a good proxy for the global solution. Figure 5 shows
the successful result and further discussion.
We also take this opportunity to demonstrate AMR. We set the adaptive
goal (16) with tolerance M = 10−6 . During the time-dependent simulation, new
cells are only added, and never removed, i.e. the mesh is never coarsened. This
is an unfortunate limitation of the current AMR algorithm. We see that cells
have been clustered near the PCI, and these clusters remain everywhere the PCI
has been during the time-dependent simulation. Rightward of the PCI, the cells
grow much larger, e.g. the right half of the domain is covered by only two cells.
Even without coarsening, already in 1D this is a large improvement. The gains
in 2D and 3D will be even greater.
5.4

2D Convection-Coupled Melting of Octadecane PCM

To demonstrate the entire coupled system, we present a preliminary result for the
convection-coupled melting of the octadecane PCM benchmark from [20] and [9].
This problem uses all aspects of our general model (9). We set variable viscosity
(5) (with µL = 1 and µS = 108 ), the buoyancy model (18) (with Ra = 3.27×105 ,
Pr = 56.2, and Re = 1), no-slip velocity boundary conditions uD = 0, and again
the temperature boundary conditions (19) with Th = 1 and Tc = −0.01. Again
we initialize the temperature field with (20) such that the simulation begins with
a thin layer of melt, and initialize a stationary velocity field u0 = 0. For the semiphase-field (4), we set regularization parameters Tr = 0.01 and r = 0.025. For the
time-dependent simulation, we set ∆t = 1. For AMR, we again set the adaptive
3

For the uniform case, this is the minimum number of cells for which the Newton
method would converge. Similarly for the AMR case, starting with a coarser initial
mesh disrupts the Newton method.

Monolithic simulation of convection-coupled phase-change

15

0.1

1
Analytical
Phaseflow, Uniform
Phaseflow, Adaptive

0.8

0.08

0.6
0.06
0.4
0.04
0.2
Analytical
Phaseflow, Uniform
Phaseflow, Adaptive

0.02

0

0

-0.2
0

0.5

1

0

0.05

0.1

Fig. 5. Left: Comparison at t = 0.1 between the analytical Stefan problem solution,
Phaseflow using a uniform mesh of 311 cells3 , and Phaseflow using AMR (with 133
cells in the final adapted mesh). The AMR solution is shown with a marker at every
mesh vertex. This emphasizes the much smaller cell sizes in the wake of the PCI.
Right: Comparison of the PCI positions at three times between the analytical, uniform,
and AMR solutions. We observe a bias in the uniform mesh solution which shrinks over
time. This is primarily because of the mesh-dependent initial PCI position (21).

goal (16) from our Stefan problem example, but with tolerance M = 10−5 until
t = 36 and M = 0.5 × 10−5 after that time.
The result in Figure 6 is promising; but compared to the benchmark, we
see that the PCI has not advanced far enough by simulated time t = 80. We
still need to investigate the effects of our choices for ∆t, r, x∗0 , M (w), M , and
the initial mesh refinement. Of these, we know that a much smaller smoothing
parameter r was used in [9].

5.5

3D Convection-Coupled Melting

The dimension-independent implementation, facilitated by the abstractions from
FEniCS, allow us to quickly demonstrate a 3D example. We consider the a
problem similar to the 2D convection-coupled melting in Section 5.4, but with
the domain extruded in the z direction, adiabatic no-slip boundary conditions
on the walls parallel to the z plane, similarity parameters Ste = 1, Ra = 106 ,
Pr = 0.71, and numerical parameters µS = 104 , r = 0.05. For the previous
examples, we employed the ”full” Newton method with ω = 1. For this problem,
we relaxed the Newton method with a factor of ω = 0.8. Figure 7 shows the
preliminary result.

16

A. G. Zimmerman, J. Kowalski

Fig. 6. Preliminary result for the 2D convection-coupled octadecane PCM melting
benchmark. Solutions are shown at simulated times t = 36 and t = 80 (on the left
and right). The mesh is shown in translucent gray to highlight the adaptive mesh
refinement (AMR). The PCI is shown as a thick black temperature isoline. Velocity
streamlines are shown as thin black lines. This result looks promising when compared
to octadecane PCM melting results from [20] and [9]. We see the flow circulating, and
the PCI advanced more quickly at the top than at the bottom. We also see that the
lack of grid coarsening is becoming expensive. By the final time, many refined cells
from earlier times are likely not needed.

Monolithic simulation of convection-coupled phase-change

17

Fig. 7. Preliminary result for 3D convection-coupled melting, with front and isometric
views (on the left and right). The mesh is shown in translucent gray to highlight the
local mesh refinement. The PCI is a darker translucent gray iso-surface. The streamlines
are colored by a grayscale from white to black as velocity magnitude increases. We use
a coarse mesh and omit AMR, because the refined 3D problem quickly exceeded the
capabilities of the desktop-scale computers used for this work. Despite the coarse mesh,
we still observe the primary features we expect in the solution. The flow is circulating,
and we see that the top of the PCI advancing more quickly than the bottom.

5.6

Convergence

To verify the accuracies of the finite difference time discretization, finite element
space discretization, and Newton linearization methods, we consider the 1D Stefan problem from Section 5.3. The mixed finite element formulation was shown
to be second order accurate for the incompressible Navier-Stokes equations in [8].
We are not aware of such a result for the energy-coupled problem. Here, as a
first step, we focus only on the energy equation with phase-change.
Based on the choices of discretizations, with fully implicit Euler for time and
finite elements for space, we expect first order convergence in time and second
order in space. Figure 8 compares this with the actual convergence orders of
Phaseflow’s solution from Section 5.3. With respect to the time step sizes δt, the
observed convergence order is only slightly higher than expected. With respect
to the grid spacing h, for a sufficiently smooth solution with well-behaved data,
we should expect second order accuracy. Though there is some deviation from
second order near h = 0.001, the results in Figure 8 show good agreement.
Next, we verify the accuracy of the Newton method. Figure 8 shows its convergence behavior for multiple time steps, highlighting the difference in performance between the first and last time steps. Ideally, Newton’s method will converge quadratically for well-posed problems with a suitable initial guess. At earlier times, a few Newton iterations pass before reaching approximately quadratic
convergence. Later times converge slightly better than quadratically.

18

A. G. Zimmerman, J. Kowalski

100

100

10-1

10-1

10-2

10-3
Actual
First Order

10-2
10-3

10-2

Actual
Second Order

10-1

10 0

10-4
10-4

10-2

100

10 -5
10-5

10 -10
10-10

10 -15
Actual
Quadratic

10 -20

0

2

4

Actual
Quadratic

6

8

10

10-15

1

2

3

4

5

Fig. 8. Convergence study results for Phaseflow applied to the Stefan problem from
Section 5.3. Top: Temporal (on the left) and spatial (on the right) convergence results
based on the solution at the final time tf .
Bottom: Newton method convergence shown for multiple time steps. Each sub-plot has
a series of lines which represent ideal (quadratic) iterations.

Monolithic simulation of convection-coupled phase-change

6

19

Conclusions and Outlook

In this work, we presented the computational tool Phaseflow for simulating the
convection-coupled melting and solidification of PCM’s. Our work is based on the
numerical method proposed in [9], and so we used an enthalpy formulated, singledomain semi-phase-field, variable viscosity, finite element method, with monolithic system coupling and global Newton linearization. The primary difference
between our numerical approach and that in [9] is the mesh adaptivity algorithm,
where we employ the dual-weighted residual method for goal-oriented AMR. We
implemented the method into our open source Python module Phaseflow, based
on the finite element software library FEniCS, and verified this against a series
of classical benchmarks. We obtained a promising result for the octadecane PCM
convection-coupled melting benchmark. Leveraging our dimension-independent
implementation, we contributed detailed convergence results for the 1D phasechange problem, and have applied the method to a preliminary 3D convectioncoupled melting example. Furthermore, we openly shared a Docker container
which allows anyone to reproduce our results in the same software environment.
It is our hope that this facilitates the further development of this and related
methods. The method appears promising, and Phaseflow is ready for application
to interesting problems.
FEniCS was a good choice as the base of our implementation, allowing us to
focus on the models and numerical methods rather than on implementing standard algorithms. Furthermore, the existing Docker software container allowed us
to quickly leverage the library. Some difficulties do remain, primarily: 1. FEniCS
lacks mesh coarsening capability, prohibiting the efficient application of AMR to
the moving PCI problem. 2. The adaptive solver (using dual-weighted residual
goal-oriented AMR) has not been implemented for distributed memory systems,
therefore currently limiting the implementation’s applicability to a single compute node, prohibiting realistic 3D applications. On the other hand, there is yet
potential to simplify Phaseflow with existing FEniCS features. Most interestingly, there is an automatic differentiation capability which could serve as an
alternative to computing the Gâteaux derivative (13) directly. During the development of Phaseflow, we have successfully applied this feature to some simplified
cases; and we would like to further explore this capability.
Next steps of our work include the development of robustness features, including... 1. a priori bounds on the time step size ∆t, based on the properties
of the linear system (12) and initial values (6), 2. an algorithm for obtaining
the initial mesh on which to obtain the first solution to begin AMR. The latter is important, because we discovered a limitation with relying on AMR to
resolve strong nonlinearities. Typically, goal-oriented AMR is used to provide
either the most accurate solution (with respect to the goal) for a given cost, or
the minimum-cost discretization for a given accuracy (with respect to the goal).
In the case of the presented method, AMR is required for the Newton method to
converge and to hence obtain the first solution in the hierarchy. Handling these
issues is necessary before practically applying this method to a wider range of
realistic problems.

20

A. G. Zimmerman, J. Kowalski

References
1. Phaseflow. https://github.com/geo-fluid-dynamics/phaseflow-fenics.
2. V Alexiades. Mathematical modeling of melting and freezing processes. CRC Press,
1992.
3. M Alnæs, J Blechta, J Hake, et al. The fenics project version 1.5. Archive of
Numerical Software, 3(100):9–23, 2015.
4. W Bangerth and R Rannacher. Adaptive Finite Element Methods for Differential
Equations. Springer Science & Business Media, 2003.
5. GK Batchelor, HK Moffatt, MG Worster, and TR Osborn. Perspectives in Fluid
Dynamics. Cambridge University Press Cambridge, 2000.
6. Y Belhamadia, A Fortin, and É Chamberland. Anisotropic mesh adaptation for
the solution of the stefan problem. J. Comput. Phys., 194(1):233–255, 2004.
7. C Boettiger. An introduction to docker for reproducible research. ACM SIGOPS
Operating Systems Review, 49(1):71–79, 2015.
8. F Brezzi and M Fortin. Mixed and hybrid finite element methods, volume 15.
Springer Science & Business Media, 1991.
9. I Danaila, R Moglan, F Hecht, and S Le Masson. A newton method with adaptive finite elements for solving phase-change problems with natural convection. J.
Comput. Phys., 274:826–840, 2014.
10. MS Dinniman, XS Asay-Davis, BK Galton-Fenzi, et al. Modeling ice shelf/ocean
interaction in antarctica: A review. Oceanography, 29(4):144–153, 2016.
11. J Donea and A Huerta. Finite element methods for flow problems. John Wiley &
Sons, 2003.
12. Y Dutil, DR Rousse, NB Salah, et al. A review on phase-change materials: mathematical modeling and simulations. Renewable and Sustainable Energy Reviews,
15(1):112–130, 2011.
13. UKNG Ghia, KN Ghia, and CT Shin. High-re solutions for incompressible flow
using the navier-stokes equations and a multigrid method. J. Comput. Phys.,
48(3):387–411, 1982.
14. HW Hsu, F Postberg, Y Sekine, et al. Ongoing hydrothermal activities within
enceladus. Nature, 519(7542):207–210, 2015.
15. DW Kelly, De SR Gago, OC Zienkiewicz, et al. A posteriori error analysis and
adaptive processes in the finite element method: Part i—error analysis. Int. J.
Numer. Methods Eng., 19(11):1593–1619, 1983.
16. J Kowalski, P Linder, S Zierke, et al. Navigation technology for exploration of
glacier ice with maneuverable melting probes. Cold Reg. Sci. Technol., 123:53–70,
2016.
17. K Schüller and J Kowalski. Integrated modeling and validation for phase change
with natural convection. Pending, 2018.
18. VR Voller, M Cross, and NC Markatos. An enthalpy method for convection/diffusion phase change. Int. J. Numer. Methods Eng., 24(1):271–284, 1987.
19. VR Voller, CR Swaminathan, and BG Thomas. Fixed grid techniques for phase
change problems: a review. Int. J. Numer. Methods Eng., 30(4):875–898, 1990.
20. S Wang, A Faghri, and TL Bergman. A comprehensive numerical model for melting with natural convection. International Journal of Heat and Mass Transfer,
53(9):1986–2000, 2010.

