\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