High-order finite element methods for moving boundary problems with prescribed boundary evolution

Highlights

• A framework for building high-order methods for moving boundary problems is presented.

• Methods are built by blending standard elements with off-the-shelf time integrators.

• A universal mesh is used to maintain an exact representation of the moving domain.

• Unlike ALE schemes, the universal mesh can handle large domain deformations easily.

Abstract

We introduce a framework for the design of finite element methods for two-dimensional moving boundary problems with prescribed boundary evolution that have arbitrarily high order of accuracy, both in space and in time. At the core of our approach is the use of a universal mesh: a stationary background mesh containing the domain of interest for all times that adapts to the geometry of the immersed domain by adjusting a small number of mesh elements in the neighborhood of the moving boundary. The resulting method maintains an exact representation of the (prescribed) moving boundary at the discrete level, or an approximation of the appropriate order, yet is immune to large distortions of the mesh under large deformations of the domain. The framework is general, making it possible to achieve any desired order of accuracy in space and time by selecting a preferred and suitable finite-element space on the universal mesh for the problem at hand, and a preferred and suitable time integrator for ordinary differential equations. We illustrate our approach by constructing a particular class of methods, and apply them to a prescribed-boundary variant of the Stefan problem. We present numerical evidence for the order of accuracy of our schemes in one and two dimensions.

Introduction

Science and engineering are replete with instances of moving boundary problems: partial differential equations posed on domains that change with time. Problems of this type, which arise in areas as diverse as fluid–structure interaction, multiphase flow physics, and fracture mechanics, are inherently challenging to solve numerically.

Broadly speaking, computational methods for moving boundary problems generally adhere to one of two paradigms. Deforming-mesh methods employ a computational mesh that deforms in concert with the moving domain, whereas fixed-mesh methods employ a stationary background mesh in which the domain is immersed. While the former approach can require that efforts be made to avoid distortions of the mesh under large deformations  [1], the latter approach requires that special care be taken in order to account for any discrepancy between the exact boundary and element interfaces  [2], [3]. Figs. 1–2 illustrate these two paradigms schematically.

In this study, we eliminate these difficulties by employing a universal mesh: a stationary background mesh that adapts to the geometry of the immersed domain by adjusting a small number of mesh elements in the neighborhood of the moving boundary. An example is illustrated in Fig. 3. The resulting framework admits, in a general fashion, the construction of methods that are of arbitrarily high order of accuracy in space and time, without exhibiting the aforementioned drawbacks of deforming-mesh and fixed-mesh methods. This strategy was introduced for time-independent and quasi-steady problems in  [4], [5]. Here we present its extension to time-dependent problems posed on moving domains with prescribed evolution. We relegate a discussion of problems with unprescribed boundaries to future work, since the treatment of unprescribed boundaries introduces its own set of challenges–approximation of the boundary, discretization of the boundary evolution equations, and error analysis on approximate domains–that may have the undesired effect of blurring the focus of the present study.

In the process of deriving our method, we present a unified, geometric framework that puts our method and existing deforming-mesh methods on a common footing suitable for analysis. The main idea is to recast the governing equations on a sequence of cylindrical spacetime slabs that span short intervals of time. The clarity brought about by this geometric viewpoint renders the analysis of numerical methods for moving-boundary problems more tractable, as it reduces the task to a standard analysis of fixed-domain problems with time-dependent PDE coefficients.

Organization. This paper is organized as follows. We begin in Section  2 by giving an informal overview of our method, and illustrating the ideas by formulating the method for a moving boundary problem in one spatial dimension. We formulate a two-dimensional model moving boundary problem on a predefined, curved spacetime domain in Section  3, and proceed to derive its equivalent reformulation on cylindrical spacetime slabs. In Section  4 we present, in an abstract manner, the general form of a finite-element discretization of the same moving boundary problem, as well as its reformulation on cylindrical spacetime slabs. This formalism will lead to a statement of the general form of a numerical method for moving boundary problems with prescribed boundary evolution that includes our method and conventional deforming-mesh methods as special cases. We finish Section  4 by summarizing an error estimate for methods of this form, referring the reader to our companion paper  [6] for its proof. In Section  5, we present the key ingredient that distinguishes our proposed method from standard approaches: the use of a universal mesh. We specialize the aforementioned error estimate to this setting to deduce that the method’s convergence rate is suboptimal by half an order when the time step and mesh spacing scale proportionately. In Section  6 we demonstrate numerically our method’s convergence rate on a prescribed-boundary variant of a classic moving-boundary problem called the Stefan problem, which asks for the evolution of a solid–liquid interface during a melting process. Some concluding remarks are given in Section  7.

Previous work. In what follows, we review some of the existing numerical methods for moving-boundary problems, beginning with deforming-mesh methods and finishing with fixed-mesh methods.

Deforming-mesh methods have enjoyed widespread success in the scientific and engineering communities, where they are best known as Arbitrary Lagrangian–Eulerian (ALE) methods. The appellation refers to the fact that in prescribing a motion of the mesh, a kinematic description of the physics is introduced that is neither Eulerian (in which the domain moves over a fixed mesh) nor Lagrangian (in which the domain does not move with respect to the mesh). The resulting formalism leads to governing equations that contain a term involving the velocity of the prescribed mesh motion that is otherwise absent in schemes on a fixed mesh  [7], [8]. Early appearances of the ALE framework date back to the works of Hirt et al.  [9], Hughes et al.  [10], and Donea et al.  [11]. ALE methods have seen use in fluid–structure interaction  [12], [13], [14], [15], [16], [17], solid mechanics  [18], [19], [20], [21], thermodynamics  [22], [23], [24], [25], and other applications.

Relative to methods for problems with fixed domains, less attention has been directed towards the development of ALE methods of high order of accuracy and the associated error analysis. Schemes of second-order in time are well-studied  [15], [16], [13], [26], [27], [28], [29], [30], though the analysis of higher-order schemes has only recently been addressed by Bonito and co-authors  [31], [32], who study the spatially continuous setting with discontinuous Galerkin temporal discretizations.

One of the key challenges that ALE methods face is the maintenance of a good-quality mesh during large deformations of the domain  [33], [34]. Fig. 1 illustrates a case where, using an intentionally naive choice of nodal motions, a domain deformation can lead to triangles with poor aspect ratios. In more severe cases, element inversions can occur. Such distortions are detrimental both to the accuracy of the spatial discretization and to the conditioning of the discrete governing equations  [35]. For this reason, it is common to use sophisticated mesh motion strategies that involve solving systems of equations (such as those of linear elasticity) for the positions of mesh nodes  [36], [37], [38], [39].

A related class of methods are spacetime methods (e.g.,  [40]), where the spacetime domain swept out by the moving spatial domain is discretized with straight or curved elements. These methods resemble deforming-mesh methods in the sense that spatial slices of the spacetime mesh at fixed temporal nodes constitute a mesh of the moving domain at those times. Bonnerot and Jamet  [41], [42] have used a spacetime framework to construct high-order methods for the Stefan problem in one dimension. They require the use of curved elements along the moving boundary to achieve the desired temporal accuracy. Jamet  [43] provides a generalization of these high-order methods to dimensions greater than one in the case that the boundary evolution is prescribed in advance. More recently, Rhebergen and Cockburn  [44], [45] created hybridizable-discontinuous-Galerkin-based spacetime methods for advection–diffusion and incompressible flow problems with moving domains.

At the other extreme are fixed-mesh methods, which cover a sufficiently large domain with a mesh and evolve a numerical representation of the boundary, holding the background mesh fixed  [46], [47], [48], [49]. A variety of techniques can be used to represent the boundary, including level sets  [50], [2], marker particles  [51], and splines  [52]. Fixed mesh methods require that special care be taken in constructing the numerical partial differential operators in the neighborhood of the moving boundary, so as to avoid losses in accuracy arising from the disagreement between the moving boundary and element interfaces. Some authors  [53], [54] propose adaptively refining the mesh in the neighborhood of the moving boundary to mitigate these losses. In the special case of a cartesian mesh, Gibou and Fedkiw  [2] have developed a third-order method for the Stefan problem in two dimensions using extrapolation to allow finite-difference stencils to extend beyond the moving boundary.

The method presented in this paper classifies neither as a deforming-mesh method nor as a fixed-mesh method, though it shares attractive features from both categories. It exhibits the immunity to large mesh distortions enjoyed by fixed-mesh methods without sacrificing the geometric conformity offered by deforming-mesh methods. Despite its conceptual simplicity, the method has not been proposed in the literature. An idea similar to ours, dubbed a “fixed-mesh ALE” method, has recently been proposed by Baiges and Codina  [55], [56], though there are several important differences. In particular, their method uses element splitting to define intermediate meshes during temporal integration, whereas our method leaves the connectivity of the mesh intact. Second, they advocate imposing boundary conditions approximately to improve efficiency; our method imposes boundary conditions exactly without extra computational effort. Finally, they focus only on low-order schemes with piecewise linear approximations to the domain deformation, while we derive schemes of arbitrarily high order.