\pdfoutput=1\relax\pdfpagewidth=8.26in\pdfpageheight=11.69in\pdfcompresslevel=9
\documentclass[reqno]{amsart}
\usepackage{graphicx} % for including postscript files
\usepackage{subfigure}
\AtBeginDocument{{\noindent\small
Fifth Mississippi State Conference on
Differential Equations and Computational Simulations,
{\em Electronic Journal of Differential Equations},
Conference 10, 2003, pp. 55--69.\newline
ISSN: 1072-6691. URL: http://ejde.math.swt.edu or http://ejde.math.unt.edu
\newline ftp ejde.math.swt.edu (login: ftp)}
\thanks{\copyright 2003 Southwest Texas State University.}
\vspace{1cm}}
\begin{document}
\setcounter{page}{55}
\title[\hfilneg EJDE/Conf/10\hfil Tin melting computations]
{Tin melting: Effect of grid size and scheme on the numerical solution}
\author[V. Alexiades, N. Hannoun, \& T. Z. Mai \hfil EJDE/Conf/10\hfilneg]
{Vasilios Alexiades, Noureddine Hannoun, \& Tsun Zee Mai}
\address{Vasilios Alexiades \hfill\break
Mathematics Department, University of Tennessee,
Knoxville, TN 37996-1300 USA, \hfill\break
and\\
Oak Ridge National Laboratory, Oak Ridge, TN 37831, USA}
\email{alexiades@utk.edu}
\address{Noureddine Hannoun \hfill\break
Mathematics Department, University of Tennessee, Knoxville, TN
37996-1300, USA}
\email{hannoun@math.utk.edu}
\address{Tsun Zee Mai \hfill\break
Mathematics Department, Box 870350,
The University of Alabama, Tuscaloosa, AL 35487, USA}
\email{tmai@gp.as.ua.edu}
\thanks{Published February 28, 2003.}
\subjclass{80A22, 74S10, 76D05, 65N99}
\keywords{Phase change, melting, numerical simulation, enthalpy method,
\hfill\break\indent benchmark, fluid flow.}
\begin{abstract}
Benchmark solutions in Computational Fluid Dynamics are necessary for
testing and for the verification of newly developed algorithms and codes.
For flows involving heat transfer coupled to solid-liquid phase change
and convection in the melt, such benchmark solutions do not exist. A
set of benchmark problems has recently been proposed for melting of
metals and waxes and several researchers responded by providing
solutions for the given problems. It was shown in two recent
publications that there were large discrepancies in the results obtained
by those contributors. In the present, work we focus on one of the four
test problems, tin melting at Rayleigh number $Ra=2.5\times 10^{5}$.
Solutions obtained for several grids and two discretization schemes
are presented and compared. Our results are used to explain the
origin of the discrepancies in earlier results. Suggestions for
future work are also provided.
\end{abstract}
\maketitle
%\date{May 30, 2001 and, in revised form, November 14, 2002.}
\section{Introduction}
Problems of \textit{ melting} frequently arise in natural and
industrial processes \cite{kn:AS93}. Applications include river, lake,
and polar melting, magma dynamics, thermal energy storage, casting and
molding, crystal growth, welding, coating, and vapor spray deposition, laser
processing, thawing of frozen food, preserved organs, tissue cultures,
and many more.
Mathematical models describing such \textit{ phase change problems} are
not amenable to standard analytical tools when flow in the melt
is involved. The problems are (de facto nonlinear) \textit{ moving
boundary problems}. For this reason, the
research community has been actively working, mostly during the last
two decades, on the development of numerical techniques suitable for
phase-change problems. This was motivated by the increasing
availability of highly powerful computers along with efficient
software.
Numerical methods dealing with phase change problems may be classified
into three main categories. The first includes \textit{ front tracking
methods}~\cite{kn:EB83,kn:C84,kn:Y93,kn:MZ94,kn:HAC95,kn:SURS96} that
consider the solid-liquid interface as a boundary between two
domains (solid and liquid) where distinct sets of conservation
equations are solved. An additional boundary condition is specified
at the interface to account for heat transfer between the two
domains. These methods require a moving adaptive grid that distorts
as time proceeds to mimic both shape and position of the interface.
Therefore, they are not suitable for problems involving highly
distorted interfaces or mergers and breakups. The second category
of numerical method includes the so called \textit{ fixed grid
techniques}~\cite{kn:VCM87,kn:VP87,kn:BVR88,kn:VS90,kn:VS91,kn:KLYM92%
,kn:SR94,kn:SURS96} for which a single set of equations is used for
the entire computational domain (solid + liquid). Transfer of latent
heat and extinction of velocities in the solid are accounted for
through the inclusion of appropriate source terms in the transport
equations. In addition to using a simpler model, fixed grid techniques
(commonly known as \textit{ enthalpy methods}) offer the highly
acknowledged benefit of the possibility of using a fixed Cartesian
grid for the entire calculations. A third category of numerical
methods known as \textit{ Eulerian-Lagrangian methods}
\cite{kn:JT96,kn:UMS99} attempts to combine features of the other
two methods.
The problem of ``\emph{Gallium melting in a rectangular enclosure
heated from the side}'' has been one of the most popular problems
for assessing the validity of newly developed techniques for phase
change coupled to convection in the liquid. This problem has received
considerable interest from both experimentalists
\cite{kn:GV86,kn:CK94,kn:CPK94} and CFD (computational fluid dynamics)
practitioners
\cite{kn:BVR88,kn:D89,kn:VJ93,kn:HC93,kn:MZ94,kn:SR94,kn:RM96,kn:SG00}.
However, all comparisons were mostly qualitative, exhibiting rather large
discrepancies among numerical methods and experimental results
\cite{kn:BVR88,kn:MZ94}. The controversial work of
Dantzig~\cite{kn:D89}, and later the illuminating work of Stella and
Giangi~\cite{kn:SG00}, have cast a doubt on the validity of past
numerical simulations. Both authors have found flow patterns with
several rolls in contrast to other publications where only a single
roll was observed. However, this seems to be contradicted by experimental
observations.
As a result of the lack of a reference solution (\textit{ benchmark solution}) for
phase change problems, Gobin and Lequere~\cite{kn:LG98} have promoted
a comparison exercise with the purpose of obtaining reference
solutions that could be used by code developers to
\textit{ verify}~\cite{kn:R98} their own code. In that exercise, four
test problems were suggested, two for paraffin waxes (high Prandtl number
$Pr$), and two for metals (low $Pr$). The first set of results has
been presented in two publications~\cite{kn:BBCC99,kn:BL00}, the emphasis
being mostly on a qualitative comparison. The authors concluded that
there were still too large discrepancies between the results from
various participants and argued it would be necessary to perform
additional simulations. They suggested that participants perform
grid refinement studies and validate their own codes on simpler
configurations. The most striking results were the fact that the less
accurate numerical solutions were actually those which best agreed
with experiments.
In this paper, we focus on one of the four test problems of the
comparison exercise~\cite{kn:BBCC99,kn:BL00}, Case 2 for tin melting at
Rayleigh number $Ra=2.5\times 10^{5}$. We present results obtained with the
same code for two discretization schemes. We want to provide an answer
to whether or not the solution presented by most contributors to the
Benchmark exercise were converged or not. We conclude the study by a
discussion about the possible origin of the discrepancies between
numerical results and experiments as well as between results of the
contributors to the comparison exercise.
\section{Physical problem}
The configuration under study is well documented
in~\cite{kn:BBCC99,kn:BL00}. A sketch of the problem is shown in
Figure~\ref{fig:benchconfig}.
\begin{figure}[!t]
\includegraphics[width=3.9in]{benchconfig.pdf}
\caption{Configuration under study : melting of tin in
a square cavity}
\label{fig:benchconfig}
\end{figure}
A square cavity (width $W=H$ height) filled with tin (pure metal) is
initially at freezing temperature $T=T_{f}$ (tin is solid). The top
and bottom boundaries are adiabatic (thermally insulated), while the
right boundary is maintained at constant temperature $T_{f}$. At time
$t=0$, the left boundary is suddenly brought to a hot temperature
$T_{h}$ larger than the melting temperature. Heat transfer by
conduction results in the melting of tin in the cavity. Thermal
gradients generate density gradients which in turn drive convection in
the liquid tin. Hence, the solid-liquid interface is no longer
planar (as for the pure conduction regime). We are interested in the
motion as well as the shape of the solid-liquid interface as melting
proceeds. The study will also focus on the flow pattern in the liquid
as well as the heat transfer process.
\section{Mathematical model -- Governing equations}
The numerical simulations are carried out based on the transport
equations~(\ref{eq:conserv}a--d) which were developed
in~\cite{kn:BVR88}.
\noindent
{\bf Mass}
\begin{subequations} \label{eq:conserv}
\begin{flalign}
& \frac{\partial \rho}{\partial t}+\nabla\cdot(\rho\vec{V}) = 0
\label{eq:modcontinuity} \tag{\ref{eq:conserv}a} \\
\intertext{\bf Momentum} & \frac{\partial \rho u}{\partial
t}+\nabla\cdot(\rho\vec{V}u) = \nabla\cdot(\mu \nabla u)-\frac{\partial P}{\partial
x}-Au \label{eq:modumoment}
\tag{\ref{eq:conserv}b} \\
& \frac{\partial \rho v}{\partial t}+\nabla\cdot(\rho\vec{V}v) = \nabla\cdot(\mu
\nabla v) -\frac{\partial P}{\partial y}-Av+\rho_{\mathrm{ref}} g
\beta(h-h_{\mathrm{ref}})/c_{p}
\label{eq:modvmoment} \tag{\ref{eq:conserv}c} \\
\intertext{\bf Energy} & \frac{\partial \rho h}{\partial t}+\nabla\cdot(\rho\vec{V}h)
= \nabla\cdot(\alpha \nabla h)-\frac{\partial \rho \Delta H}{\partial
t}-\nabla\cdot(\rho\vec{V} \Delta H)
\label{eq:modenergy} \tag{\ref{eq:conserv}d}
\end{flalign}
\end{subequations}
In these equations, $\rho$ stands for density, $T$ for temperature,
$P$ for pressure, $\vec{V}$ for velocity vector, $u$ and $v$ for $x$
and $y$ components of velocity, $h$ for sensible enthalpy, $\Delta H$
for latent enthalpy content, $\mu$ for dynamic viscosity, $\lambda$ for
thermal conductivity, $c_{p}$ for specific heat, $\alpha$ for
$\lambda/c_{p}$, $\beta$ for coefficient of thermal expansion, $g$ for
gravitational acceleration, $t$ for time and $x$ and $y$ for Cartesian
coordinates. The subscript ``$\mathrm {ref}$'' is used for reference
quantities.
The flow is assumed to two dimensional, unsteady, and obeys the
Navier-Stokes equations for incompressible Newtonian fluids written in
Cartesian coordinates (eq.~\ref{eq:conserv}a-c). The thermophysical
properties ($c_{p}$, $\lambda$, $\alpha$, $\mu$, $\beta$) are independent
of temperature and are the same for both solid and liquid. The model also
assumes that density $\rho$ may vary with temperature. However the
present simulations assume a constant density along with the
Boussinesq approximation (density variations due to temperature gradients
are accounted for in buoyancy terms).
The present formulation is a one-domain method wherein the same set of
equations is used for both solid and liquid. The material in the
cavity is regarded as a porous medium with porosity varying with
liquid fraction through Carman-Kozenay's law. The
constant $A$ in the source term of the momentum Equations~%
(\ref{eq:modumoment}, \ref{eq:modvmoment}) has the form:
\begin{equation}
A = -\frac{C(1-f_{L})^{2}}{f_{L}^{3}+q}
\label{eq:modcarmankozeny}
\end{equation}
where $f_{L}$ is the liquid fraction and $C$ and $q$ are two
constants chosen to ensure driving the velocities to zero in the solid,
while maintaining a convergent algorithm. The source term $S_{h}$:
\begin{equation}
S_{h}= - \: \frac{\partial \rho \Delta H}{\partial t}-\nabla\cdot(\rho\vec{V}
\Delta H)
\label{eq:modenergysource}
\end{equation}
in the right hand side of the energy
equation~(\ref{eq:modenergy}) accounts for latent enthalpy transfer
during phase change.
\section{Numerical method}
A finite volume method is used to discretize
Equations~(\ref{eq:conserv}a--d) on a Cartesian uniform grid with
staggered arrangement for the velocities~\cite{kn:HW65}. Each equation
is integrated on a control volume centered at a node of the main
variable for that equation. Second order accuracy is retained for
quadratures, source terms and diffusion terms. Convective fluxes are
approximated with the formalism of Patankar~\cite{kn:P80} which allows
one to choose among five different schemes (upwind, hybrid, centered,
exponential, and power law). Time discretization is fully implicit
(Euler Backward). Nonlinearity and coupling between the various equations is
handled by the SIMPLER algorithm of Patankar~\cite{kn:P88}.
The energy equation requires a special treatment due to the presence
of latent enthalpy content in the source terms arising from
phase change. To obtain the new enthalpy value $h^{k+1}$
for the current outer iteration ($k+1$), one needs $\Delta H^{k+1}$
which itself depends on the unknown solution $h^{k+1}$. The process
for handling this problem is described in~\cite{kn:BVR88}. At each
outer iteration ($k+1$) the latent enthalpy content is updated from
the values at the previous outer iteration ($k$) through the
formula:
\begin{equation}
\Delta H^{k+1}=\Delta H^{k}+\omega_{\Delta
H}\frac{a_{P}^{h}}{a_{P}^{oh}}
\biggl\{h^{k}-c_{p}\biggl[\frac{T_{L}-T_{S}}{L}\Delta
H^{k}+T_{S}\biggr]\biggr\}
\label{eq:DHupdate}
\end{equation}
In equation~(\ref{eq:DHupdate}), $a_{P}^{h}$ and $a_{P}^{oh}$ stand
for the central coefficient of the discretized energy equation and the
unsteady term coefficient respectively~\cite{kn:BVR88}, while
$\omega_{\Delta H}$ is a
relaxation factor used to avoid divergence of outer iterations. Latent
enthalpy content $\Delta H$ is related to temperature $T$ through a
linear relationship over a small interval $\epsilon=T_{L}-T_{S}$ where
$T_{L}$ and $T_{S}$ are the liquidus and solidus temperatures
respectively. for $T>T_{L}$ the latent enthalpy content is equal to
the latent heat, $L$ while for $T