A numerical scheme for the Stefan problem on adaptive Cartesian grids with supralinear convergence rate

doi:10.1016/j.jcp.2009.04.044

Journal of Computational Physics

Volume 228, Issue 16, 1 September 2009, Pages 5803-5818

https://doi.org/10.1016/j.jcp.2009.04.044 Get rights and content

Abstract

We present a level set approach to the numerical simulation of the Stefan problem on non-graded adaptive Cartesian grids, i.e. grids for which the size ratio between adjacent cells is not constrained. We use the quadtree data structure to discretize the computational domain and a simple recursive algorithm to automatically generate the adaptive grids. We use the level set method on quadtree of Min and Gibou [C. Min, F. Gibou, A second order accurate level set method on non-graded adaptive Cartesian grids, J. Comput. Phys. 225 (2007) 300–321] to keep track of the moving front between the two phases, and the finite difference scheme of Chen et al. [H. Chen, C. Min, F. Gibou, A supra-convergent finite difference scheme for the poisson and heat equations on irregular domains and non-graded adaptive Cartesian grids, J. Sci. Comput. 31 (2007) 19–60] to solve the heat equations in each of the phases, with Dirichlet boundary conditions imposed on the interface. This scheme produces solutions that converge supralinearly $(\sim 1.5)$ in both the $L^{1}$ and the $L^{\infty}$ norms, which we demonstrate numerically for both the temperature field and the interface location. Numerical results also indicate that our method can simulate physical effects such as surface tension and crystalline anisotropy. We also present numerical data to quantify the saving in computational resources.

Introduction

The Stefan problem is a moving interface model where the main physical process is diffusion, and is at the center of the study of crystal growth. Moreover, it has important applications in the growing field of semi-conductors, in combustion, in bio-nanotechnology, in tissue engineering and many others. The difficulty in solving the Stefan problem stems from the fact that the interface position must be computed as part of the solution. In addition, boundary conditions, some of which are depending on the interface shape, must be satisfied at the evolving front.

Successful numerical methods for this type of problems need to be able to efficiently solve the heat equation on irregular domains and keep track of a moving interface that may undergo complex topological changes such as the merging and the pinching of two fronts. The interface that separates the two phases can be either explicitly tracked or implicitly captured. The main advantage of explicit approaches, e.g. front tracking [18], is their accuracy. The main disadvantage is that special care is needed for topological changes. In turn, the explicit treatment of connectivity makes the method challenging to extend to three spatial dimensions. Implicit representations such as the level set method [29], [30], [35] and the phase-field method [20] represent the front as an iso-contour of a continuous function. Topological changes are consequently handled in a natural fashion, and thus these methods are readily implemented in both two and three spatial dimensions.

In this paper we use the level set method, first introduced by Osher and Sethian [30], [29], [35], which is a sharp interface model that can be used to exactly locate the interface in order to apply discretizations that depend on the exact interface location. This method can handle discontinuous material properties easily. The earliest level set method for solidification-type problems was presented by Sethian and Strain [34], where the authors recast the equations of motion into a boundary integral formulation and used the level set method to update the location of the interface. Chen et al. [11] introduced a level set approach where the diffusion equation was solved directly with finite difference schemes. A similar approach was proposed by Udaykumar et al. [39]. In both cases, the solution is second-order accurate and the resulting linear system is non-symmetric. Kim et al. [19] demonstrated that this method produces results in agreement with the solvability theory [23]. Gibou et al. proposed a second-order accurate symmetric discretization [15], [14] and third-order nonsymmetric discretization [13]. The methods presented in all these works are on uniform grids.

Physical phenomena modelled by the Stefan problem present differences in length scales and numerical approximations on uniform grids are in such cases extremely inefficient in terms of C.P.U. and memory requirement since only a small fraction of the domain needs high grid resolution, while other parts of the domain can produce accurate solutions on much coarser grids. In fact, since the Stefan problem is described by a parabolic partial differential equation, the solution is smooth far from the interface where coarse grids are sufficient to adequately capture the variations of the solution. On the other hand, the solution presents a kink near the interface, due to Dirichlet boundary conditions. It is therefore desirable to design a scheme that refines automatically near the interface in order to capture the sharp features in the solution’s gradients, while leaving coarser cells far away.

Adaptive mesh refinement techniques are advantageous in those cases: AMR (adaptive mesh refinement) on Cartesian grids has been pioneered by the work of Berger and Oliger [9] and Berger and Colella [8] for compressible flows. Within this framework, the domain is first discretized as a coarse grid and then rectangular blocks of uniform grids with finer resolution are added where the solution requires higher accuracy. Several applications for elliptic, parabolic and hyperbolic systems have been proposed [4], [5], [6], [17], [22]. This line of work was also applied to two-phase flows in Sussman et al. [37].

Implementations based on recursive data structures have become more popular in recent years. These approaches do not impose patches of uniform grids but rather allows the grid cells to continuously change in size. Quadtree (in 2D) and octree (in 3D) data structures [32], [33] are convenient in representing such grids, and have been proven to be optimal [1]. Recently, Min et al. [27] solved the variable coefficient Poisson equation on a rectangular domain with non-graded adaptive Cartesian grids, i.e. grids for which the size ratio between adjacent cells is not constrained, and obtained second-order accurate solutions with second-order accurate gradients in the $L^{1}$ and $L^{\infty}$ norms. This work was then extended in Chen et al. [10] to irregular domains with non-graded Cartesian grids to achieve second-order accuracy for both the solution and its gradients. The advantage of non-graded grids is that mesh generation is efficient and straightforward. It also save computational resources.

In this paper, we build on the work of [10], [26] to present a level set approach to the numerical simulation of the Stefan problem on adaptive Cartesian grids. We note that the extremely simple refinement scheme we are using produces grids that are mostly graded, but not necessarily. We emphasize that the numerical method we present is applicable to non-graded Cartesian grids and that other refinement criteria, potentially more efficient or more adapted to particular problems, can be used. We demonstrate the supra-linear convergence of our method in the $L^{1}$ and the $L^{\infty}$ norms and demonstrate that it can simulate physical effects such as thermal conductivity, crystalline anisotropy, surface tension, molecular kinetics and undercoolings. Finally, we illustrate the efficiency and accuracy by comparing with computations on uniform grids.

Section snippets

Equations for the Stefan problem

Consider a domain $Ω = Ω^{-} \cup Ω^{+}$ , where the two subdomains $Ω^{-}$ and $Ω^{+}$ are separated by an interface $Γ$ . The boundary of $Ω$ is denoted by $\partial Ω$ (see Fig. 1). The Stefan problem describes the evolution of a scalar T, equals to $T_{s}$ in $Ω^{-}$ and $T_{l}$ in $Ω^{+}$ , such that: $\frac{\partial T_{s}}{\partial t} = \nabla \cdot (D_{s} \nabla T_{s}) in Ω^{-},$ $\frac{\partial T_{l}}{\partial t} = \nabla \cdot (D_{l} \nabla T_{l}) in Ω^{+},$ where the subscripts s and l denote the solid and liquid phases, respectively. The diffusion constants $D_{s}$ and $D_{l}$ can be discontinuous across $Γ$ . The temperature on the solid–liquid interface satisfies the following

Automatic grid generation

We represent the grid in two spatial dimensions using the quadtree data structure [32], [33]: Referring to Fig. 2, the root of the tree corresponds to the entire domain and we add four children cells as we split the domain in four equal parts. This process is repeated recursively in order to refine parts of the domain where interesting features develop.

By definition, a quadtree is graded if the difference between two adjacent cell levels is at most one, i.e. a cell can be at most twice as big

Treating T-junction nodes

The main difficulty in designing numerical schemes on non-graded adaptive Cartesian grids comes from T-junctions nodes, i.e. nodes for which adjacent nodes are not necessarily aligned in a Cartesian direction. In Min et al. [27], [25], [26], we introduced a new approach to define the values at these ghost nodes that automatically produces second-order accurate schemes. In particular, referring to Fig. 3, we define a ghost node $T_{4}$ as: $T_{4} = \frac{s_{5} T_{6} + s_{6} T_{5}}{s_{5} + s_{6}} - \frac{s_{5} s_{6}}{s_{2} + s_{3}} (\frac{T_{2} - T_{0}}{s_{2}} + \frac{T_{3} - T_{0}}{s_{3}}) .$ The last term in

Accuracy and efficiency

We discuss the efficiency and accuracy of our numerical method on the known Frank-sphere exact solution [12]: In two spatial dimensions, the region described by a disk with zero temperature, is growing into the supercooled liquid. The radius of the growing disk is given by $R (t) = s_{0} \sqrt{t}$ , parameterized by $s_{0}$ . The temperature field is given by $T (r, t) = T (s) = \{\begin{matrix} 0 & if s ⩽ s_{0} \\ T_{\infty} (1 - \frac{F (s)}{F (s_{0})}) & if s > s_{0}, \end{matrix}$ where r is the distance to the center of the disk, $s = r / \sqrt{t}, F (s) = E_{1} (s^{2} / 4)$ , with $E_{1} (z) = \int_{z}^{\infty} \frac{e^{- t}}{t} dt$ . Here, $F (s)$ is the similarity

Conclusions

We have presented a level set approach to the numerical simulation of the Stefan problem on non-graded adaptive Cartesian grids, i.e. grids for which the size ratio between adjacent cells is not constrained. Compared with uniform grids, adaptive grids relocate computational resources efficiently by refining only where needed, thus reducing both the amount of computational time and the memory usage significantly. Compared with traditional graded adaptive grids, the non-graded property of our

Acknowledgments

The research of H. Chen and F. Gibou was supported in part by a Sloan Research Fellowship in Mathematics, by the National Science Foundation Under Grant Agreement DMS 0713858 and by the Department of Energy Under Grant Agreement DE-FG02-08ER15991. The research of C. Min was supported by the Korea Research Foundation Grant funded by the Korean Government (KRF-2008-331-C00045).

References (40)

A. Almgren et al.
A conservative adaptive projection method for the variable density incompressible Navier–Stokes equations
J. Comput. Phys.
(1998)
A. Almgren et al.
A fast adaptive vortex method in three dimensions
J. Comput. Phys.
(1994)
T. Aslam
A partial differential equation approach to multidimensional extrapolation
J. Comput. Phys.
(2004)
M. Berger et al.
Local adaptive mesh refinement for shock hydrodynamics
J. Comput. Phys.
(1989)
M. Berger et al.
Adaptive mesh refinement for hyperbolic partial differential equations
J. Comput. Phys.
(1984)
S. Chen et al.
A simple level set method for solving Stefan problems
J. Comput. Phys.
(1997)
F. Gibou et al.
A fourth order accurate discretization for the Laplace and heat equations on arbitrary domains, with applications to the Stefan problem
J. Comput. Phys.
(2005)
F. Gibou et al.
A second-order-accurate symmetric discretization of the Poisson equation on irregular domains
J. Comput. Phys.
(2002)
H. Johansen et al.
A Cartesian grid embedded boundary method for Poisson’s equation on irregular domains
J. Comput. Phys.
(1998)
D. Juric et al.
A front tracking method for dendritic solidification
J. Comput. Phys.
(1996)

P. McCorquodale et al.

A node-centered local refinement algorithm for Poisson’s equation in complex geometries

J. Comput. Phys.

(2004)

C. Min et al.

A second order accurate projection method for the incompressible Navier–Stokes equation on non-graded adaptive grids

J. Comput. Phys.

(2006)

C. Min et al.

A second order accurate level set method on non-graded adaptive Cartesian grids

J. Comput. Phys.

(2007)

C. Min et al.

A supra-convergent finite difference scheme for the variable coefficient Poisson equation on non-graded grids

J. Comput. Phys.

(2006)

C.-H. Min

Local level set method in high dimension and codimension

J. Comput. Phys.

(2004)

S. Osher et al.

Fronts propagating with curvature-dependent speed: algorithms based on Hamilton–Jacobi formulations

J. Comput. Phys.

(1988)

J. Sethian et al.

Crystal growth and dendritic solidification

J. Comput. Phys.

(1992)

C.-W. Shu et al.

Efficient implementation of essentially non-oscillatory shock capturing schemes

J. Comput. Phys.

(1988)

M. Sussman et al.

An adaptive level set approach for incompressible two-phase flow

J. Comput. Phys.

(1999)

M. Sussman et al.

A level set approach for computing solutions to incompressible two-phase flow

J. Comput. Phys.

(1994)

Cited by (58)

A metric-based adaptive mesh refinement criterion under constrain for solving elliptic problems on quad/octree grids
2024, Journal of Computational Physics
In this work we propose and investigate the performance of a metric-based refinement criteria for adaptive meshing used for improving the numerical solution of an elliptic problem. We show that in general, when solving elliptic equations such as the Poisson-Helmholtz equation, the minimization of the interpolation error often used as local refinement criteria does not always guarantee the minimization of the total numerical error. Numerical and theoretical arguments are given to unveil the critical role of the mesh compression – the size aspect ratio between the finest cell size and the mean cell size of an adapted mesh – to determine whether the estimated error is purely local meaning that the interpolation error is a good enough error model for the total error or if other, non-local, sources of error need to be accounted for. We show through particular examples that a slightly sub-optimal mesh in terms of interpolation error may significantly reduce the total error of a numerical solution, depending on the value of the compression ratio and not on the number of grid points. Based on this observation, we propose a new method to exclude the grids where non-local errors can control the accuracy of the solution. This is achieved by an automatic estimation of the optimal compression ratio, which is imposed as an additional constraint in the minimal element size during the mesh adaptation process. The method is tested on quadtree and octree grids, showing very satisfactory performances in reducing the total numerical error despite the additional constrain imposed.
Machine learning algorithms for three-dimensional mean-curvature computation in the level-set method
2023, Journal of Computational Physics
We propose a data-driven mean-curvature solver for the level-set method. This work is the natural extension to $R^{3}$ of our two-dimensional error-correcting strategy presented in Larios-Cárdenas and Gibou (October 2022) [1] and the hybrid inference system of Larios-Cárdenas and Gibou (August 2022) [2]. However, in contrast to [1], [2], which built resolution-dependent dictionaries of curvature neural networks, here we develop a single pair of neural models in $R^{3}$ , regardless of the mesh size. The core of our system comprises two feedforward neural networks. These models ingest preprocessed level-set, gradient, and curvature information to fix numerical mean-curvature approximations for selected nodes along the interface. To reduce the problem's complexity, we have used the local Gaussian curvature $κ_{G}$ to classify stencils and fit these networks separately to non-saddle (where $κ_{G} ⪆ 0$ ) and saddle (where $κ_{G} ⪅ 0$ ) input patterns. The Gaussian curvature is essential for enhancing generalization and simplifying neural network design. For example, non-saddle stencils are easier to handle because they exhibit a mean-curvature error distribution characterized by monotonicity and symmetry. While the latter has allowed us to train only on half the mean-curvature spectrum, the former has helped us blend the machine-learning-corrected output and the baseline estimation seamlessly around near-flat regions. On the other hand, the saddle-pattern error structure is less clear; thus, we have exploited no latent characteristics beyond what is known. Such a stencil distinction makes learning-tuple generation in $R^{3}$ quite involving. In this regard, we have trained our models on not only spherical but also sinusoidal and hyperbolic paraboloidal patches at various configurations. Our approach to building their data sets is systematic but gleans samples randomly while ensuring well-balancedness. Furthermore, we have resorted to standardization and dimensionality reduction as a preprocessing step and integrated layer-wise regularization to minimize outlying effects. In addition, our strategy leverages mean-curvature rotation and reflection invariance to improve stability and precision at inference time. The synergy among all these features has led to a performant mean-curvature solver that works for any grid resolution. Experiments with several interfaces confirm that our proposed system can yield more accurate mean-curvature estimations in the $L^{1}$ and $L^{\infty}$ norms than modern particle-based interface reconstruction and level-set schemes around under-resolved regions. Although there is ample room for improvement, this research already shows the potential of machine learning to remedy geometrical problems in well-established numerical methods. Our neural networks are available online at https://github.com/UCSB-CASL/Curvature_ECNet_3D.
Adjoint-based optimization of two-dimensional Stefan problems
2023, Journal of Computational Physics
A range of optimization cases of two-dimensional Stefan problems, solved using a tracking-type cost-functional, is presented. A level-set method is used to capture the interface between the liquid and solid phases and an immersed boundary (cut cell) method coupled with an implicit time-advancement scheme is employed to solve the heat equation. A conservative implicit-explicit scheme is then used for solving the level-set transport equation. The resulting numerical framework is validated with respect to existing analytical solutions of the forward Stefan problem. An adjoint-based algorithm is then employed to efficiently compute the gradient used in the optimization algorithm (L-BFGS). The algorithm follows a continuous adjoint framework, where adjoint equations are formally derived using shape calculus and transport theorems. A wide range of control objectives are presented, and the results show that using parameterized boundary actuation leads to effective control strategies in order to suppress interfacial instabilities or to maintain a desired crystal shape.
Quasi-equilibrium approximation of the solidification process for micro phase change materials taking into account curvature and surface tension
2023, International Journal of Thermal Sciences
Phase change materials (PCMs) have superior energy storage properties, which can be used to solve the problem of energy absorption and release at various times, intensities and locations. Energy storage using PCMs is a promising way to make the rational use of energy, which could help improving our environment. Owing to their significant applications, we investigate the solidification process of microparticles, which is also known as the Stefan problem. Due to the tiny size of microparticles, quasi-equilibrium approximation of such Stefan’s problem is considered, where the effect of the curvature and the surface tension are incorporated. Here, we determine a simple but yet effective iterative model for the moving boundaries of micro PCMs (MPCMs) by assuming the homogeneous temperature distribution in molten regimes in order to avoid solving more complicated diffusion equations. We purposely investigate the solidification process of spherical and ellipsoidal MPCMs and find that both the curvature and the surface tension promote the solidification process. Due to these outcomes, peanut-shaped MPCMs are proposed to demonstrate the usage of both effects for fabricating smart MPCMs. Finally, the validity of the model is confirmed by a Computational Fluid Dynamics (CFD) method.
A sharp numerical method for the simulation of Stefan problems with convective effects
2022, Journal of Computational Physics
We present a numerical method for the solution of interfacial growth governed by the Stefan model coupled with incompressible fluid flow. An algorithm is presented which takes special care to enforce sharp interfacial conditions on the temperature, the flow velocity and pressure, and the interfacial velocity. The approach utilizes level-set methods for sharp and implicit interface tracking, hybrid finite-difference/finite-volume discretizations on adaptive quadtree grids, and a pressure-free projection method for the solution of the incompressible Navier-Stokes equations. The method is first verified with numerical convergence tests using a synthetic solution. Then, computational studies of ice formation on a cylinder in crossflow are performed and provide good quantitative agreement with existing experimental results, reproducing qualitative phenomena that have been observed in past experiments. Finally, we investigate the role of varying Reynolds and Stefan numbers on the emerging interface morphologies and provide new insights around the time evolution of local and average heat transfer at the interface.
Error-correcting neural networks for semi-Lagrangian advection in the level-set method
2022, Journal of Computational Physics
Citation Excerpt :
In order to refine a grid efficiently, researchers have proposed nonuniform discretization schemes to increase the resolution only next to the interface [47,48]. Adaptive Cartesian grids, for instance, have served as the basis for high-order robust level-set algorithms and tools in both serial [49–55] and distributed computing systems [56,57]. Over the years, practitioners have also combined level-set technologies with other numerical frameworks to improve mass conservation.
We present a machine learning framework that blends image super-resolution technologies with passive, scalar transport in the level-set method. Here, we investigate whether we can compute on-the-fly, data-driven corrections to minimize numerical viscosity in the coarse-mesh evolution of an interface. The proposed system's starting point is the semi-Lagrangian formulation. And, to reduce numerical dissipation, we introduce an error-quantifying multilayer perceptron. The role of this neural network is to improve the numerically estimated surface trajectory. To do so, it processes localized level-set, velocity, and positional data in a single time frame for select vertices near the moving front. Our main contribution is thus a novel machine-learning-augmented transport algorithm that operates alongside selective redistancing and alternates with conventional advection to keep the adjusted interface trajectory smooth. Consequently, our procedure is more efficient than full-scan convolutional-based applications because it concentrates computational effort only around the free boundary. Also, we show through various tests that our strategy is effective at counteracting both numerical diffusion and mass loss. In simple advection problems, for example, our method can achieve the same precision as the baseline scheme at twice the resolution but at a fraction of the cost. Similarly, our hybrid technique can produce feasible solidification fronts for crystallization processes. On the other hand, tangential shear flows and highly deforming simulations can precipitate bias artifacts and inference deterioration. Likewise, stringent design velocity constraints can limit our solver's application to problems involving rapid interface changes. In the latter cases, we have identified several opportunities to enhance robustness without forgoing our approach's basic concept. Despite these circumstances, we believe all the above assets make our framework attractive to parallel level-set algorithms. Its appeal resides in the possibility of avoiding further mesh refinement and decreasing expensive communications between computing nodes.

View all citing articles on Scopus

View full text

Published by Elsevier Inc.