Optimal convergence of the arbitrary Lagrangian–Eulerian interface tracking method for two-phase Navier–Stokes flow without surface tension

Abstract

1. Introduction

Immiscible fluids mixture separated by a moving interface appears widely in nature and engineering applications. Current models for two-phase or multiphase flow can generally be categorized into two main types: diffuse interface models and sharp interface models. Diffuse interface models treat the interface as an extremely thin transition zone that separates the two fluids. In these models, physical quantities like density and viscosity change rapidly, but smoothly within this transition zone. On the other hand, sharp interface models consider the interface as a surface without any thickness. In such models, physical quantities exhibit discontinuities across this surface.

We consider the sharp interface model of two-phase Navier–Stokes (NS) flow Gross & Reusken (2011) without surface tension in a bounded domain |$\varOmega (t)=\varOmega _{+}(t)\cup \varOmega _{-}(t)\cup \varGamma (t)$| in |${{\mathbb{R}}}^{d}$|⁠, with |$d\in \{2,3\}$|⁠, where |$\varGamma (t)=\overline \varOmega _{+}(t)\cap \overline \varOmega _{-}(t)$| is a sharp interface which separates the two fluids occupying two subdomains |$\varOmega _{+}(t)$| and |$\varOmega _{-}(t)$|⁠, respectively; see Fig. 1. The fluid velocity |$u$| and pressure |$p$| are governed by NS equations in each subdomain |$\varOmega _\pm (t)$| and satisfy certain jump conditions on the interface |$\varGamma (t)$|⁠, i.e.,

Fig. 1.

Domain |$\varOmega $| in the case |$d = 2$|⁠.

Open in new tabDownload slide

where |$ \rho _{\pm } $| are the densities of the fluids in |$\varOmega _{\pm }(t)$|⁠, |$\nu $| is the unit normal vector on |$\varGamma (t)$| pointing to |$\varOmega _{+}(t)$|⁠; |$\sigma = (2\mu _{\pm } {{\mathbb{D}}}(u) -p{{\mathbb{I}}})$| is the stress tensor, in which |${{\mathbb{D}}}(u) = \frac 12(\nabla u + (\nabla u)^\top )$| and |${{\mathbb{I}}}$| denote the deformation matrix and identity matrix, respectively; |$\mu _{\pm }$| are the viscosities of the two fluids in |$\varOmega _{\pm }(t)$|⁠, respectively; |$[u]_{-}^{+} = u_{+} - u_{-}$| and |$[\sigma \nu ]_{-}^{+} = \sigma _{+}\nu - \sigma _{-}\nu $| denote the jumps of the quantities |$u$| and |$\sigma \nu $| across the interface, respectively. The gravitational force is given by |$f = (0,-g)^\top $| in two dimensions and |$f = (0,0, -g)^\top $| in three dimensions with |$g$| a constant (the gravitational acceleration).

The interface |$\varGamma (t)$| moves along with the fluid and therefore also has velocity |$u$|⁠. If |$\varGamma ^{0}$| is the initial interface, through a flow map |$X:\varGamma ^{0}\times [0,T]\rightarrow{{\mathbb{R}}}^{d}$|⁠, we can describe the evolution of interface |$\varGamma (t)=\{X(y,t):y\in \varGamma ^{0}\}$| by the following equation:

We do not take surface tension into account on the interface |$\varGamma (t)$|⁠, which differs from models that include surface tension, where the appearance of curvature leads to a coupling between the geometric properties of the interface and the fluid flow in the bulk. The fluid dynamics in (1.1)–(1.2) involve the movement of the interface |$\varGamma (t)$| and two subdomains |$\varOmega _{\pm }(t)$| over time, while the overall domain |$\varOmega = \varOmega (t)$| remains unchanged in shape.

Numerical methods for (1.1) can generally be classified into three categories, i.e., the Eulerian approach, the Lagrangian approach and the arbitrary Lagrangian–Eulerian approach (ALE). The ALE approach serves as a bridge between the Lagrangian and Eulerian approaches, by enabling the frame to move with an ‘arbitrary’ bulk velocity that fits the interface motion. At the discrete level, the mesh points can be displaced independently of the flow. Utilizing the ALE technique allows for flexible movement of the inner domain mesh, while the mesh on the interface can move alongside materials to precisely track the interfaces of a multi-material system.

In an early investigation of ALE methods, the stability analysis of the ALE finite element method (FEM) was firstly conducted by Formaggia and Nobile Formaggia & Nobile (1999) for a linear parabolic equation. Subsequently, Gastaldi Gastaldi (2001) developed a priori error estimates in the |$L^{2}$| norm for ALE-FEM when solving parabolic equations. Many studies on the numerical analysis of ALE-FEMs for parabolic moving boundary/interface problems can be found in works such as Nobile (2001), Boffi & Gastaldi (2004), Formaggia & Nobile (2004), Badia & Codina (2006), Kovács & Power Guerra (2018), and Lan & Sun (2020). Optimal-order error bounds of |$O(h^{k+1})$| in the |$L^\infty (0, T; L^{2})$| norm of ALE semidiscrete FEM with curved triangles/simplices for parabolic moving boundary problems were established by Elliott and Ranner Elliott & Ranner (2021). We also refer to Edelmann (2022) for an ALE method with harmonically evolving mesh. Optimal-order |$H^{1}$| convergence of the ALE-FEM for PDEs coupling boundary evolution arising from shape optimization problems was proved in Gong et al. (2023). These results were established for high-order curved evolving mesh. Optimal convergence of |$O(h^{k+1})$| in the |$L^\infty (0, T; L^{2})$| norm, with flat evolving simplices in the interior and curved simplices exclusively on the boundary, was proved in Li et al. (2023) for the ALE semidiscrete FEM utilizing a standard iso-parametric element of degree |$k$| in Lenoir (1986). The ALE methods for PDEs in bulk domains Gong et al. (2023) are also closely related to the evolving FEMs for PDEs on evolving surfaces, for which optimal-order convergence in the |$L^{2}$| and |$H^{1}$| norms were shown in Dziuk & Elliott (2007), Elliott & Venkataraman (2015) and Kovács et al. (2017). The above-mentioned research efforts have focused on diffusion equations.

The numerical analysis of ALE methods for the Stokes and NS equations has also noteworthy results, but remained suboptimal. Legendre and Takahashi Legendre & Takahashi (2008) established an |$L^{2}$| error estimate of |$O(\tau + h^{\frac 12})$| for the method of characteristics combined with a |$P_{1b}-P_{1}$| finite element approximation to the ALE formulation of 2D NS equations. In San Martín et al. (2009), by using the Taylor–Hood |$P_{2}$|-|$P_{1}$| element, the error estimates of |$O(h^{2}|\log h|)$| and |$O(\tau + h^{2}+ h^{2}/\tau )$| were proved for the semi-discrete and fully discrete ALE-FEM for the Stokes equations in a time-dependent domain, respectively. The error of ALE finite element solutions to the Stokes equations on a time-varying domain, with BDF-|$r$| in time (for |$1\leqslant r\leqslant 5$|⁠) and the Taylor–Hood |$P_{k}$|–|$P_{k-1}$| elements in space (with degree |$k\geqslant 2$|⁠), was proved to be |$O(\tau ^{r} + h^{k})$| in the |$L^{2}$| norm in Liu (2013). Optimal-order error bounds of |$O(h^{k+1})$| and |$O(\tau ^{2}+h^{k+1})$| of the semidiscrete and fully discrete ALE-FEMs in the |$L^{2}$| norm for the Stokes equations were proved recently in (Rao et al.).

Many efforts have also been made in developing efficient ALE-FEMs for two-phase NS flow problems, primarily for models with surface tension. For example, Anjos et al. in Anjos et al. (2014) introduced a 3D ALE-FEM scheme with adaptive mesh to achieve good control of the mesh quality by using repeated interpolations of the velocity. Barrett et al. in Barrett et al. (2015a) also considered surface tension and designed a stable ALE-FEM in which the mesh points used to describe the interface are not mesh points of the underlying bulk finite element mesh. An implicit tangential velocity of the interface nodes is introduced in Barrett et al. (2015a) to acquire good mesh quality of the interface. Barrett, Garcke and Nürnberg introduced a novel weak formulation in Barrett et al. (2020), which is referred to as the BGN formulation from now on. The BGN formulation allows for tangential degrees of freedom to improve the mesh quality. This formulation was employed for two-phase Stokes flow in Barrett et al. (2013) and two-phase NS flow in Barrett et al. (2015b), both with surface tension and with unfitted finite element approximations, which leads to an unconditional stability bound. Moreover, the applications of BGN formulation to the two-phase NS flow with surface tension with moving fitted FEMs were also investigated in Agnese & Nürnberg (2020) in both Eulerian and ALE approaches, but no stability bound was available. Fu in Fu (2020) proposed an ALE-HDG FEM for two-phase NS flow with surface tension with an implicit tangential velocity to acquire good mesh quality of the interface, by using HDG to produce an exactly divergence-free velocity approximation. Recently, Duan et al. in Duan et al. (2022) proposed an energy-diminishing ALE scheme with surface tension by using Stokes extension of the velocity from the interface to the bulk domain, in order to produce an exactly divergence-free mesh velocity. Energy stability of the ALE-FEMs with surface tension was discussed in Barrett et al. (2013), Barrett et al. (2015a), Barrett et al. (2015b), Fu (2020) and Duan et al. (2022). We also refer to Garcke et al. (2023) for the structure-preserving ALE-FEMs for the two-phase NS flow with surface tension in both the conservative and non-conservative forms. The methods introduced in Garcke et al. (2023) were shown to satisfy unconditional stability and exact volume preservation on the fully discrete level. The numerical investigation of convergence behavior (without error analysis) with surface tension, can be found in Duan et al. (2022) and Fu (2020). As far as we know, there is no error analysis of ALE-FEMs for the sharp interface model of two-phase NS flow problems with surface tension. In this work, we focus on a simplified model of two-phase NS flow without surface tension. Two major difficulties are discussed below.

First, the error analysis of ALE-FEMs for parabolic and Stokes/NS equations in the above-mentioned articles all require the triangulation of the domain to fit the moving boundary or interface exactly, i.e., the triangulated domain |$\varOmega _{h}(t)$| for computation that interpolates the exact domain |$\varOmega (t)$|⁠. This is the case when the boundary or interface of domain |$\varOmega (t)$| is given, determined either by a given flow map or by a given velocity field. However, this is the not the case in two-phase NS flow where the velocity of the interface is unknown (determined by the solution of the PDE problem). Thus the location of the exact interface is unknown and needs to be approximated based on the numerical solution of the two-phase NS flow. As a result, the numerically computed interface at |$t>0$| generally has a distance from (therefore is not an interpolation of) the exact interface, and the triangulated subdomains |$\varOmega _{h,\pm }(t)$| are not interpolations of the exact subdomains |$\varOmega _{\pm }(t)$| at |$t>0$|⁠. Consequently, the error analysis of ALE-FEMs in the pre-existing articles for parabolic and Stokes/NS equations with a given moving boundary/interface cannot be applied to the two-phase NS flow without surface tension.

Second, the error analysis of ALE-FEMs for parabolic and Stokes/NS equations in the above-mentioned articles all require the solution to be sufficiently smooth in the whole domain, while the solution of the two-phase NS equation in (1.1) is generally not smooth in the whole domain due to the existence of the free interface, which separates the two phases of the fluid, e.g., |$\nabla u$| is generally discontinuous across the interface |$\varGamma (t)$|⁠. Since the numerical interface and the exact interface differ from each other by a certain distance at |$t>0$|⁠, the triangles of the mesh (which fit the numerical interface) generally do not fit the exact interface |$\varGamma (t)$|⁠. Thus the exact solution |$u$| may be nonsmooth inside a triangle of the mesh. This situation renders it difficult to prove high-order convergence using standard interpolation onto the mesh.

The aim of this article is to establish an optimal-order error bound in the |$L^\infty (0, T; H^{1}(\varOmega ))$| norm for the semidiscrete ALE-FEM for the two-phase NS flow problem without surface tension in (1.1) using high-order curved evolving mesh, where the mesh points on the numerical interface move with the fluid’s velocity and is harmonically extended to the bulk subdomains. We address the two major difficulties mentioned above by using the matrix-vector formulations of the ALE-FEM and the error equations. Such techniques were originally introduced in Kovács et al. (2017) for analyzing solution-driven surface evolutions. It was later developed for analyzing the convergence of evolving FEMs for surface evolution under geometric flows Kovács et al. (2019, 2021); Hu & Li (2022); Bai & Li (2024) and domain evolution under shape gradient flows in Gong et al. (2023). We adapt such matrix-vector formulation techniques here to address the mismatch between the numerical interface and the exact interface, using an intermediate auxiliary triangulated domain, as well as the missing of global regularity of the solution by using lift operations which map the numerical subdomains to the exact subdomains. This approach only requires using piecewise smoothness of the solution |$u$| on the subdomains |$\varOmega _\pm (t)$| instead of its global smoothness of |$u$| on the whole domain |$\varOmega (t)$|⁠. For a smoothly evolving sharp interface |$\varGamma (t)$| without topological changes, and a solution |$(u,p)\in C([0,T];H^{1}(\varOmega )^{d})\times L^{2}(0,T;L^{2}(\varOmega ))$|⁠, which is piecewise smooth in the subdomains |$\varOmega _\pm (t)$| for |$t\in [0,T]$|⁠, we prove optimal-order error bound of |$O(h^{k})$| in the |$C([0,T];H^{1}(\varOmega )^{d})\times L^{2}(0,T;L^{2}(\varOmega ))$| norm for the semi-discrete ALE-FEMs scheme with Taylor–Hood |$P_{k}$|–|$P_{k-1}$| elements of degree |$k \geqslant 2$|⁠; see Theorem 4.1.

The rest of this article is organized as follows: Section 2 introduces the ALE reformulation of (1.1) and the ALE weak formulation on the evolving domain. In Section 3, we introduce the ALE weak formulation and ALE evolving FEM on the interpolated evolving domain and provide a consistency analysis. In Section 3, we formulate the consistency error estimates based on interpolated evolving mesh and matrix-vector formulation, which serve, not only for implementation purposes, but also play a key role in the subsequent convergence analysis. In Section 4, we prove optimal-order convergence of the ALE evolving FEM considered in this article. Numerical experiments are presented in Section 5 to illustrate the convergence of the ALE evolving FEM and its performance on a benchmark example.

2. The ALE evolving FEM

2.1 The ALE reformulation of PDE

Let |$w(\cdot ,t): \varOmega (t) \rightarrow{{\mathbb{R}}}^{d}$| be the harmonic extension of fluid’s velocity |$u$| from the interface |$\varGamma (t)$| to the bulk subdomains |$\varOmega _\pm (t)$|⁠, i.e.,

Harmonic extension is a common method for generating mesh velocity |$w$| in the bulk subdomains; see Boffi & Gastaldi (2004), Fu (2020) and Edelmann (2022). Since |$w=u$| on |$\varGamma (t)$|⁠, it follows that |$w$| and |$u$| generate the same evolution of interface |$\varGamma (t)$| and subdomains |$\varOmega _\pm (t)$|⁠.

Denote the initial domain by |$ \varOmega ^{0} = \varOmega (0) $| and the initial subdomains by |$ \varOmega _\pm ^{0} = \varOmega _\pm (0) $|⁠. Let |$\phi :\varOmega ^{0}\times [0, T]\rightarrow{{\mathbb{R}}}^{d}$| be the flow map generated by the velocity field |$w$|⁠. Then |$\phi $| maps |$\varOmega _{\pm }^{0}$| to |$\varOmega _\pm (t)$| and |$\varGamma ^{0}$| to |$\varGamma (t)$|⁠, respectively, with

and satisfies the following ordinary differential equation:

We consider an ALE moving frame associated with the flow map |$\phi $|⁠, with frame velocity |$w$|⁠. Then the original problem in (1.1) can be written as

where |$\partial _{t}^{\bullet } u = \partial _{t} u + w \cdot \nabla u$| is the material derivative with respect to the frame velocity |$w$|⁠. The subdomains |$\varOmega _\pm (t)$| evolve with velocity |$w$|⁠.

We assume that problem (1.1), or equivalently (2.3), has a solution

which is piecewise smooth in subdomains |$\varOmega _\pm (t)$| for |$t \in (0, T]$|⁠, where

By piecewise smooth we mean that |$\sum _{j=0}^{2}\|(\partial _{t}^{\bullet })^{j} u\|_{H^{k+1}(\varOmega _\pm (t))}+\sum _{j=0}^{1}\|(\partial _{t}^{\bullet })^{j} p \|_{H^{k}(\varOmega _\pm (t))}$| is bounded uniformly with respect to |$t\in [0,T]$|⁠.

We assume that the flow map |$\phi (\cdot , t) \in C^{k}(\bar \varOmega _{\pm }^{0})$| is piecewise smooth in the subdomains |$\varOmega _{\pm }^{0}$|⁠. By pulling back to the initial domain |$\varOmega ^{0}$| after composition with |$\phi (\cdot , t)$| and utilizing the regularity estimate of the elliptic equation (Evans, 2022, Theorem 5, Section 6.3) on the fixed domain |$\varOmega ^{0}$|⁠, the frame velocity |$w$| defined in (2.1) satisfies

which is also piecewise smooth. Then, by using notations

testing (2.3) by |$v$| and using the interface and boundary conditions in (1.1), we obtain the following weak reformulation:

ALE weak formulation: find |$(u, p) \in C([0, T]; H_{0}^{1}(\varOmega )^{d}) \times C([0, T]; L_{0}^{2}(\varOmega )) $| with |$\partial ^{\bullet }_{t} u \in C([0, T]; L^{2}(\varOmega )^{d}) $|⁠, |$w\in C([0,T]; H^{1}(\varOmega _\pm )^{d})$| and |$\varOmega _\pm (t)=\phi (\varOmega _\pm ^{0},t)$|⁠, such that

with initial conditions |$\phi (0)=\mathrm{id}|_{\varOmega ^{0}}$| and |$u(0)=u^{0}$|⁠.

2.2 The ALE mapping and ALE evolving FEM

Assume that the given smooth initial domain |$\varOmega ^{0}\subset{{\mathbb{R}}}^{d}$|⁠, |$d\in \{2,3\}$|⁠, is divided into a set |${{\mathcal{K}}}_{h}^{\ \, 0}$| of shape-regular and quasi-uniform simplices with mesh size |$h$|⁠, using curved finite elements of degree |$k$| and fitting the interface |$\varGamma ^{0} = \varGamma (0)$|⁠. Every possibly curved closed simplex |$K \in{{\mathcal{K}}}_{h}^{\ \,0}$| is the image of a unique polynomial of degree |$k$| defined on the reference simplex |$\hat K$|⁠, denoted by |$F_{K}:\hat K\rightarrow K$|⁠, called the parametrization of |$K$|⁠; see (Elliott & Ranner, 2021, §8.6). Every boundary simplex |$K\in{{\mathcal{K}}}_{h}^{\ \,0}$| (with one face or edge attached to |$\varGamma ^{0}$|⁠) contains a possibly curved face or edge to interpolate |$\varGamma ^{0}$| with accuracy of |$O(h^{k+1})$|⁠. This can be obtained by using Lenoir’s isoparametric finite elements Lenoir (1986). Thus the initial domain |$\varOmega ^{0}$| is approximated by the triangulated domain

Let |$\mathbf{x}^{0}=(\xi _{1},\cdots ,\xi _{M})\in{{\mathbb{R}}}^{dM}$| be the vector that collects all the nodes |$\xi _{j}\in{{\mathbb{R}}}^{d}$|⁠, |$j=1,\ldots ,M$|⁠, (which define finite elements of degree |$k$|⁠) in the triangulation |${{\mathcal{K}}}_{h}^{0}$|⁠. We evolve the vector |$\mathbf{x}^{0}$| in time and denote its position at time |$t$| by

which determines the triangulation |${{\mathcal{K}}}_{h}[\mathbf{x}(t)]$| and domain |$\varOmega _{h}[\mathbf{x}(t)]$|⁠, as well as subdomains |$\varOmega _{h,\pm }[\mathbf{x}(t)]$| and interface |$\varGamma _{h}[\mathbf{x}(t)] $| such that

Similarly as the simplices on the initial domain, if |$K \in{{\mathcal{K}}}_{h}[\mathbf{x}(t)]$| is a simplex on the evolving domain |$\varOmega _{h}[\mathbf{x}(t)] $|⁠, then we denote by |$F_{K}:\hat K\rightarrow K$| the parametrization of |$K$|⁠. For the simplicity of notation, we denote

Correspondingly, the scalar-valued finite element space on the evolving domain |$\varOmega _{h}(t)$| is defined as

The ALE mapping |$\phi _{h}(\cdot , t)$| is defined as the unique finite element function |$\phi _{h}(\cdot ,t) \in S_{h}^{k}[\mathbf{x}^{0}]^{d}$| such that

which maps the initial domain |$\varOmega _{h, \pm }^{0}$| onto the evolving domain |$\varOmega _{h, \pm }(t)$|⁠, and the initial interface |$\varGamma _{h}^{0}$| onto the evolving interface |$\varGamma _{h}(t)$|⁠, i.e.,

The ALE mapping |$\phi _{h}(\cdot , t)$| is the flow map generated by a mesh velocity |$w_{h}(\cdot , t)$| and satisfies the following relation:

Based on the ALE mapping |$\phi _{h}(\cdot , t)$| and the mesh velocity |$w_{h}(\cdot , t)$|⁠, if |$v(\cdot , t)$| is a function defined on |$\varOmega _{h}(t)$|⁠, its discrete material derivative is defined as

Let |$X_{h}[\mathbf{x}(t)] \times M_{h}[\mathbf{x}(t)] \subset H_{0}^{1}(\varOmega _{h}(t))^{d} \times L_{0}^{2}(\varOmega _{h}(t)) $| denote the Taylor–Hood finite element space of degree |$k\geqslant 2$| subject to the triangulation |${{\mathcal{K}}}_{h}[\mathbf{x}(t)]$| of evolving domain |$\varOmega _{h}(t)$|⁠, i.e.,

We also denote by |$\overset{ \circ }{X}_{h}[\mathbf{x}(t)] $| the subspace of |$X_{h}[\mathbf{x}(t)]$| with zero interface and boundary conditions, i.e.,

It is known that the Taylor–Hood finite element space satisfies the discrete inf-sup condition (see (Boffi & Fortin, 2013, §8.5.3) on the domain decomposition techniques), i.e., there exists positive constant |$\kappa $| independent of |$h$| such that

where

Since the mesh velocity |$w_{h}$| should be an approximation of |$w$| defined in Section 2.1, we define |$w_{h}(\cdot , t)$| to be a discrete harmonic finite element function in the numerically computed subdomains |$\varOmega _{h, \pm }(t)$| with boundary condition matching the velocity of the fluid on the boundary |$\partial \varOmega _{h}(t)$| and the numerically computed interface |$\varGamma _{h}(t)$|⁠. This is equivalent to finding |$w_{h}\in u_{h} + \overset{ \circ }{X}_{h}[\mathbf{x}(t)]$| such that

We denote by |$\rho _{h}$| and |$\mu _{h}$| the piecewise constant functions in |$\varOmega _{h}(t)$|⁠, defined by

The ALE evolving FEM for (2.4) is defined as follows:

ALE evolving FEM: find a flow map |$\phi _{h}(\cdot ,t) \in S_{h}^{k}[\mathbf{x}^{0}]^{d} $|⁠, which determines nodal vector |$\mathbf{x}(t) = \phi _{h}(\mathbf{x}^{0}, t)$| and corresponding subdomains |$\varOmega _{h,\pm }(t)=\varOmega _{h,\pm }[\mathbf{x}(t)]$|⁠, as well as |$(u_{h}(\cdot ,t), p_{h}(\cdot ,t) ) \in X_{h}[\mathbf{x}(t)] \times M_{h}[\mathbf{x}(t)]$| and |$w_{h}(\cdot ,t) \in u_{h}(\cdot ,t) + \overset{ \circ }{X}_{h}[\mathbf{x}(t)]$| such that

with initial conditions |$\phi _{h}(0) = \mathrm{id}|_{\varOmega _{h}^{0}}$| and |$u_{h}(0)=u_{h}^{0}$|⁠, where |$u_{h}^{0}$| is the Lagrange interpolation of |$u^{0}$| and |$ f_{h} $| is the Lagrange interpolation of |$ f $| on |$\varOmega _{h}(t) $|⁠.

3. Interpolated evolving mesh and consistency analysis

In the error analysis, we need to compare functions on two different domains, i.e., the exact domain |$\varOmega (t)$|⁠, and the triangulated domain |$\varOmega _{h}(t)$| determined by the numerical solution. The consistency error is the defect upon inserting the Lagrange interpolation of the exact solution |$u$| and the interpolated interface into the discretized equations. To this end, we define the interpolated ALE mapping and interpolated evolving mesh below.

3.1 Interpolated ALE mapping and interpolated evolving mesh

We evolve the vector |$\mathbf{x}^{0}=(\xi _{1},\cdots ,\xi _{M})$| in time and denote |$\mathbf{x}^{*}(t)=(x_{1}^{*}(t),\dots ,x_{M}^{*}(t))$| the nodal vector which collects the points |$x_{j}^{*}(t) = \phi (\xi _{j}, t)$| on the exact domain |$\varOmega (t)$|⁠, which determines triangulation |${{\mathcal{K}}}_{h}[\mathbf{x}^{*}(t)]$| and interpolated domain |$\varOmega _{h}^{*}(t)$|⁠. Since the initial triangulation is shape-regular and quasi-uniform, the interpolated triangulations |${{\mathcal{K}}}_{h}[\mathbf{x}^{*}(t)]$| will keep shape-regular and quasi-uniform for |$t\in [0,T]$| when the flow map |$\phi (\cdot ,t):\varOmega _\pm (0)\rightarrow \varOmega _\pm (t)$| and its inverse are piecewise smooth on |$\varOmega _\pm (0)$| and |$\varOmega _\pm (t)$| for |$t\in [0,T]$|⁠, respectively.

Subject to the triangulation |${{\mathcal{K}}}_{h}[\mathbf{x}^{*}(t)]$| of domain |$\varOmega _{h}^{*}(t)$|⁠, the space of scalar-valued finite element functions on |$\varOmega _{h}^{*}(t)$| is defined as

We define the interpolated ALE mapping |$\phi _{h}^{*}(\cdot ,t):\varOmega _{h}^{0} \to \varOmega _{h}^{*}(t)$| to be the Lagrange interpolation of flow map |$\phi (\cdot , t): \varOmega (0) \to \varOmega (t)$|⁠, which has the following properties:

and

If we denote by |$w_{h}^{*}(\cdot , t)= I_{h}^{*}(t) w(\cdot , t)$| the Lagrange interpolations of the exact velocity |$w(\cdot , t)$|⁠, then the following relation holds:

This equality holds because both sides are finite element functions on |$\varOmega _{h}^{0}$| and they are equal on the nodes.

Let |$u_{h}^{*}$| be the Lagrange interpolation of the exact solution |$u$|⁠, and denote by |$w_{h,\pm }^{*}$| and |$u_{h,\pm }^{*}$| the restriction of |$w_{h}^{*}$| and |$u_{h}^{*}$| to subdomains |$\varOmega _{h,\pm }^{*}(t)$|⁠. Then the interpolated mesh velocity satisfies |$w_{h}^{*}=u_{h}^{*}$| on the interpolated interface |$\varGamma _{h}^{*}(t)$|⁠. Thus |$w_{h,\pm }^{*}-u_{h,\pm }^{*}\in \overset{ \circ }{X}_{h, \pm }[\mathbf{x}^{*}(t)] $|⁠.

For a function |$v(\cdot ,t)$| defined on |$\varOmega _{h}^{*}(t)$|⁠, we can define the discrete material derivative |$\partial _{t,*}^{\bullet }$| by using the interpolated ALE mapping |$\phi _{h}^{*}(\cdot ,t)$| and the interpolated mesh velocity |$w_{h}^{*}(\cdot , t)$| as

3.2 Lifting operator

The transition between the interpolated domain |$\varOmega _{h}^{*}(t)$| and the exact domain |$ \varOmega _{h}(t)$| is done by a lift operator. In Lenoir’s isoparametric finite element approximation to |$\varOmega ^{0}$|⁠, there exists a map |$\varPhi :\varOmega _{h}^{0}\rightarrow \varOmega ^{0}$| such that (cf. Lenoir (1986))

For the triangulation |${{\mathcal{K}}}_{h}[\mathbf{x}^{*}(t)]$| that fits the interface |$\varGamma (t)$|⁠, there exists a continuous mapping |$\varPhi _{h}^{*}(\cdot ,t): \varOmega _{h}^{*}(t)\rightarrow \varOmega (t)$| given by

that maps |$\varGamma _{h}^{*}(t)$| to |$\varGamma (t)$| and satisfies the following estimates:

Correspondingly, for |$v_{h}: \varOmega _{h}^{*}(t) \to{{\mathbb{R}}}$|⁠, we denote its lift by |$v_{h}^{\ell }$| given by

• (Lagrange interpolation) Define |$\hat I_{h, \hat K}$| be the Lagrange interpolation on the reference simplex |$\hat K$|⁠, then the following relation hold:

  $$\begin{align} (I_{h}^{*}(t) v) \circ F_{K} = \hat I_{h, \hat K} \big(v\circ \varPhi_{h}^{*}(\cdot, t) \circ F_{K} \big).\end{align}$$

  (3.6)
• (Commutation of material derivative and lift) (cf. (Elliott & Ranner, 2021, Lemma 3.5)):

  $$\begin{align} \partial^{\bullet}_{t} \big(v_{h}^{\ell}\big) = \big(\partial^{\bullet}_{t, *} v_{h} \big)^{\ell} \qquad \mbox{for} \quad v_{h}: \varOmega_{h}^{*}(t) \to{{\mathbb{R}}}.\end{align}$$

  (3.7)
• (Transport property of the lifted basis functions |$\varphi _{j} [\mathbf{x}^{*}]^{\ell }$|⁠) (cf. (Dziuk & Elliott, 2013, Lemma 4.1)):

  $$\begin{align} \partial^{\bullet}_{t} (\varphi_{j} [\mathbf{x}^{*}(t)]^{\ell}) = \big(\partial^{\bullet}_{t, *} \varphi_{j} [\mathbf{x}^{*}(t)] \big)^{\ell} = 0, \quad j = 1, \dots, M.\end{align}$$

  (3.8)
• (Commutation of material derivative of Lagrange interpolation): Utilizing the transport property of the lifted basis functions, we have

  $$\begin{align} \partial^{\bullet}_{t, *} I_{h}^{*}(t) v = I_{h}^{*}(t) \partial^{\bullet}_{t} v.\end{align}$$

  (3.9)

  Define |${{\mathcal{I}}}_{h}^{*} u:= (I_{h}^{*} u)^{\ell }: \varOmega (t) \to{{\mathbb{R}}}$|⁠. It follows from (3.8) and (3.7) that

  $$\begin{align} \partial^{\bullet}_{t} {{\mathcal{I}}}_{h}^{*} u = {{\mathcal{I}}}_{h}^{*} \partial^{\bullet}_{t} u.\end{align}$$

  (3.10)

  It plays a crucial role in proving the optimal-order in consistency error.
• (Approximation properties for the Lagrange interpolation) There exists a constant |$C> 0$| independent of |$h \leqslant h_{0}$|⁠, with a sufficiently small |$h_{0}> 0$|⁠, and |$t$| such that for |$u(\cdot ,t) \in H^{k+1}(\varOmega _{\pm }(t))$|⁠, for |$0\leqslant t \leqslant T$| (cf. (Kovács et al., 2017, Lemma 7.5)):

  $$\begin{align} \|u - {{\mathcal{I}}}_{h}^{*}(t) u \|_{L^{2}(\varOmega_{\pm}(t))} + h \|\nabla u - \nabla{{\mathcal{I}}}_{h}^{*}(t) u \|_{L^{2}(\varOmega_{\pm}(t))} \leqslant C h^{k+1} \|u\|_{H^{k+1} (\varOmega_{\pm}(t))}.\end{align}$$

  (3.11)
• (Equivalence of norms) The |$L^{p}$| and |$W^{1, p}$| norms on the discrete and continuous domains are equivalent for |$1 \leqslant p \leqslant \infty $|⁠, uniformly in the mesh size |$h \leqslant h_{0}$| (with sufficiently small |$h_{0}>0$|⁠) and in |$t \in [0,T]$|⁠. In particular, for |$v_{h}: \varOmega _{h}^{*}(t) \to{{\mathbb{R}}}$| with lift |$v_{h}^{\ell }: \varOmega (t) \to{{\mathbb{R}}}$|⁠, there is a constant |$C$| such that for |$h \leqslant h_{0}$| and |$0 \leqslant t \leqslant T$|⁠,

  $$\begin{align} &C^{-1}\|v_{h}\|_{L^{2}(\varOmega_{h}^{*}(t))} \leqslant \|v_{h}^{\ell}\|_{L^{2}(\varOmega(t))} \leqslant C\|v_{h}\|_{L^{2}(\varOmega_{h}^{*}(t))}, \end{align} $$

  (3.12)

   

  $$\begin{align} &C^{-1}\|v_{h}\|_{H^{1}(\varOmega_{h}^{*}(t))} \leqslant \|v_{h}^{\ell}\|_{H^{1}(\varOmega(t))} \leqslant C\|v_{h}\|_{H^{1}(\varOmega_{h}^{*}(t))}. \end{align} $$

  (3.13)

  The lift operator |$\ell $| maps a function on the interpolated domain |$\varOmega _{h}^{*}(t)$| to a function on the exact domain |$\varOmega (t)$|⁠, provided that |$\varOmega _{h}^{*}(t)$| is sufficiently close to |$\varOmega (t)$|⁠.

In the following, we drop the argument |$t$| when it is not essential.

3.3 Weak formulation on the interpolated evolving mesh

For the simplicity of notation, we denote by |$\rho _{h}^{*}$| and |$\mu _{h}^{*}$| the piecewise constant functions in |$\varOmega _{h}^{*}(t)$| defined by

and define the following inner product:

After replacing the solution by its finite element interpolation in weak formulation (2.4), and replacing domains |$\varOmega (t)$| and |$\varOmega _{\pm }(t)$| by the interpolated domains |$\varOmega _{h}^{*}(t)$| and |$\varOmega _{h,\pm }^{*}(t)$|⁠, respectively, we obtain the following weak formulation satisfied by the interpolated functions |$u_{h}^{*}$|⁠, |$p_{h}^{*}$|⁠, |$w_{h}^{*}$| and |$\phi _{h}^{*}$|⁠:

with initial conditions |$\phi _{h}^{*}(0) =\mathrm{id}|_{\varOmega _{h}^{0}}$| and |$u_{h}^{*}(0)=u_{h}^{0}$|⁠, and |$ f_{h}^{*} $| is the Lagrange interpolation of |$ f $| on |$\varOmega _{h}^{*}(t) $|⁠, where |${{\mathcal{E}}}_{u1}(t, v_{h})$|⁠, |${{\mathcal{E}}}_{u2}(t, q_{h})$| and |${{\mathcal{E}}}_{w}(t, \chi _{h})$| are consistency errors arising from finite element approximations to the solution and domain. For example,

which contains errors arising from approximating domain |$\varOmega _\pm (t)$| by |$\varOmega _{h,\pm }^{*}(t)$|⁠. This can be controlled by using the following estimates for the error caused by domain perturbation (see Kovács et al. (2017)):

where |$\chi _{h},\psi _{h}\in S_{h}^{k}[\mathbf{x}^{*}(t)]$| and |$g\in L^{2}(\varOmega _{h}^{*}(t))$|⁠. By using this result and the estimates of the interpolation error shown in (3.11), the following estimates of the consistency errors can be derived when the domain |$\varOmega (t)$| and interface |$\varGamma (t)$| are smooth, and the solution is piecewise smooth in each subdomain |$\varOmega _\pm (t)$|⁠:

for |$v_{h}\in X_{h}[\mathbf{x}^{*}(t)] $|⁠, |$q_{h}\in M_{h}[\mathbf{x}^{*}(t)] $| and |$\chi _{h}\in \overset{ \circ }{X}_{h}[\mathbf{x}^{*}(t)] $|⁠. In the case |$\partial _{t,*}^{\bullet } v_{h}=0$| and |$\partial _{t,*}^{\bullet } q_{h}=0$| on |$\varOmega _{h}^{*}(t)$|⁠, the derivatives of these consistency errors also satisfy similar estimates, i.e.,

3.4 Matrix-vector form of the weak formulation

For both computation and error analysis it is convenient to use the matrix-vector form of the finite element weak formulation in (2.8). Through the utilization of finite element basis functions |$\varphi _{j}[\mathbf{x}(t)]$| on the domain |$\varOmega _{h}(t)$|⁠, the numerical solution |$u_{h}(t)$| can be expressed as a column vector |$\mathbf{u} = (u_{1}(t), \dots , u_{M}(t))$|⁠, where

Similarly, the solutions |$p_{h}(t)$| and |$w_{h}(t)$| on |$\varOmega _{h}(t)$| can also be represented as column vector |$\mathbf{p} $| and |$\mathbf{w}$|⁠, respectively. Then |$\mathbf{x}(t)=(x_{1}(t),\dots ,x_{M}(t)) =\phi _{h}(\mathbf{x}^{0},t)$| is the vector of positions of all finite element nodes at |$t$|⁠, satisfying

where the initial value |$\mathbf{x}^{0}$| is given by the triangulation of the domain |$\varOmega $| at time |$t=0$|⁠.

We introduce domain-dependent mass matrix |$\mathbf{M}(\mathbf{x})$|⁠, the stiffness matrix |$\mathbf{A}(\mathbf{x})$| and the matrix |$ \mathbf{M}_{\boldsymbol{\rho }} (\mathbf{x})$| and |$\mathbf{A}_{\boldsymbol{\mu }}(\mathbf{x})$| on the domain |$\varOmega _{h}(t)$| as follows:

With identity matrix |$I_{d} \in{{\mathbb{R}}}^{d\times d}$|⁠, we define |$\mathbf{K}(\mathbf{x}^{n})^{d}$| as the Kronecker product of |$I_{d}$| and |$\mathbf{K}(\mathbf{x}^{n}):=\mathbf{M}(\mathbf{x}^{n})+\mathbf{A}(\mathbf{x}^{n})$|⁠, i.e.,

To simplify the notation, we will use |$\mathbf{K}(\mathbf{x}^{n})$| to represent |$\mathbf{K}(\mathbf{x}^{n})^{d}$| when the dimension of the matrix is clear and therefore no confusion arises. With the matrices defined above, the |$L^{2}$| and |$H^{1}$| norm of finite element functions can be expressed as quadratic forms:

In the same way, we can define the matrices |$\mathbf{B}_{\boldsymbol{\rho }}(\mathbf{x}, \mathbf{u})$| and |$\mathbf{C}(\mathbf{x})$| such that

In situations where there is no ambiguity in the context, we omit the time dependency and rewrite the weak formulation (2.8) into the following matrix-vector form:

where the column vector |$ \mathbf{f} $| collects the nodal values of |$ f_{h} $|⁠. In the same way, we can define the corresponding matrices and vectors over the interpolated domain |$\varOmega _{h}^{*}(t)$|⁠. The matrices and vectors denoted with a superscript * collect the nodal values of the interpolation of the exact solutions. For example, the vector |$\mathbf{v}^{*}$| collects the nodal values of |$v_{h}^{*}$|⁠. With this convention, we can rewrite the weak formulation (3.16) in the following form:

where |$\boldsymbol{{\mathcal{E}}}_{\mathbf{u}1}$|⁠, |$\boldsymbol{{\mathcal{E}}}_{\mathbf{u}2}$|⁠, |$\boldsymbol{{\mathcal{E}}}_{\mathbf{w}}$| denote the vectorized forms of the defect terms |${{\mathcal{E}}}_{u1}(t, \cdot )$|⁠, |${{\mathcal{E}}}_{u2}(t, \cdot )$| and |${{\mathcal{E}}}_{w}(t, \cdot )$|⁠, which satisfy the following relations:

for |$v_{h}^{*}\in X_{h}[\mathbf{x}^{*}(t)]$|⁠, |$q_{h}^{*} \in M_{h}[\mathbf{x}^{*}(t)]$|⁠, |$\chi _{h}^{*} \in \overset{ \circ }{X}_{h}[\mathbf{x}^{*}(t)]$| and the associated nodal vectors |$\mathbf{v}$|⁠, |$\mathbf{q}$|⁠, |$\boldsymbol{\chi }$|⁠.

Convergence of the ALE evolving FEM is established in this article by comparing the matrix-vector formulations in (3.24) and (3.25). By subtracting (3.25) from (3.24), we obtain the following equations for the errors |$\mathbf{e}_{\mathbf{w}} = \mathbf{w} - \mathbf{w}^{*}$|⁠, |$\mathbf{e}_{\mathbf{u}} = \mathbf{u}- \mathbf{u}^{*}$|⁠, |$\mathbf{e}_{\mathbf{p}} = \mathbf{p} - \mathbf{p}^{*}$| and |$\mathbf{e}_{\mathbf{x}} = \mathbf{x} - \mathbf{x}^{*}$|⁠:

for all |${\boldsymbol \chi }$| associated to some |$\chi _{h}\in \overset{ \circ }{X}_{h}[\mathbf{x}^{*}(t)]$|⁠.

The estimates of the consistency errors in (3.17)–(3.19) can also be written into the matrix-vector form, as shown in the following proposition. This proposition generalizes the consistency estimates found in (Kovács et al., 2017, Lemma 8.1), as it directly follows from the estimates for the error caused by domain perturbation and the approximation properties of Lagrangian interpolation.

4. Convergence of the ALE evolving FEM

We are now in the position to formulate the main result of this paper, which yields optimal-order convergence in the |$H^{1}$| norm for the Taylor–Hood finite element semidiscretization of problem (3.24).

The proof of Theorem 4.1 will be presented in the next two subsections based on the technique of a homotopy map between two different domains |$\varOmega _{h}^{*}(t)$| and |$\varOmega _{h}(t)$|⁠.

4.1 Comparison of integrals on two different domains

To establish the convergence of the numerical schemes, it is necessary to compare the integrals over two different domains, i.e, the triangulated domain |$\varOmega _{h}^{*}(t)$| obtained from interpolating the exact domain |$\varOmega (t)$|⁠, and the triangulated domain |$\varOmega _{h}(t)$| determined by the numerical solution. In the spirit of the techniques in Kovács et al. (2017); Edelmann (2022), we can obtain a similar sequence of results employing shape derivatives by constructing a homotopy map between |$\varOmega _{h}^{*}(t)$| and |$\varOmega _{h}(t)$|⁠.

To handle the rate of change of an integral over a moving domain, we will frequently make use of the following lemma (cf. (Walker, 2015, Lemma 5.7)):

We establish a linear continuous deformation from |$\varOmega _{h}^{*}(t)$| to |$\varOmega _{h}(t)$|⁠. Using the basis functions |$\varphi _{j}[\mathbf{x}(0)]$| of |$S^{k}_{h}(\varOmega _{h}^{0})^{d}$|⁠, we uniquely determine the domains |$\varOmega _{h}^{*}(t)$| and |$\varOmega _{h}(t)$| by vectors |$\mathbf{x}^{*}$| and |$\mathbf{x}$|⁠, respectively. The vector

then defines an intermediate domain |$\varOmega _{h}^{\theta }(t)$| that changes continuously from |$\varOmega _{h}^{*}(t)$| to |$\varOmega _{h}(t)$| when the parameter |$\theta $| varies over the interval |$[0,1]$|⁠. For the vector |$\mathbf{u}^{*} = (u_{1}^{*}(t), \dots , u_{M}^{*}(t))$|⁠, we define the finite element function on |$\varOmega _{h}^\theta (t)$| as

Similarly, we can define functions |$e_{x}^{\theta }$|⁠, |$e_{u}^{\theta }, e_{w}^\theta $|⁠, |$w_{h}^{*, \theta }$| and |$(\partial _{t,*}^{\bullet } u_{h}^{*})^\theta $|⁠, which are defined on intermediate domain |$\varOmega _{h}^\theta (t)$| and share the nodal vectors |$\mathbf{e}_{\mathbf{x}}, \mathbf{e}_{\mathbf{u}}, \mathbf{e}_{\mathbf{w}}$|⁠, |$\mathbf{w}^{*}$| and |$\dot{\mathbf{u}}^{*}$|⁠, respectively. These notations are consistently applied to other functions without redundant introductions.

Let |$\partial _{\theta }^{\bullet }$| denote the material derivative with respect to the velocity field |$e_{x}^{\theta }$|⁠. By using the transport property, we have

Then the following relations follow from (4.8) and (4.5):

We denote by |$\rho _{h}^\theta $| and |$\mu _{h}^\theta $| the piecewise constant functions in |$\varOmega _{h}^\theta (t)$| defined by

In combination with the fundamental theorem of calculus, Lemma 4.1 and the transport property (2.9), the following lemma was proved in (Edelmann, 2022, Lemma 5.1).

A direct consequence of Lemma (4.2) is the following conditional equivalence of norms, i.e., if |$\nabla e_{x}^{\theta }$| is small, the norms of finite element functions on the two domains |$\varOmega _{h}^{*}(t)$| and |$\varOmega _{h}(t)$| with same nodal vectors are equivalent. This lemma was proved in (Edelmann, 2022, Lemma 5.2) by using the Gronwall’s inequality.

The following lemma was proved in (Edelmann, 2022, Lemma 5.3), which says that the condition for |$e_{x}^\theta $| in Lemma 4.3 can be reduced to |$\theta = 0$|⁠.

We also need results that bound the time derivatives of the mass and stiffness matrices. The following results can be shown in a similar way as the analogous results on a surface (Kovács et al., 2017, Lemma 4.6).

4.2 The proof of Theorem 4.1

Let |$ e_{x}^{*}\in S_{h}^{k}[\mathbf{x}^{*}(t)]^{d}$|⁠, |$e_{u}^{*} \in X_{h}[\mathbf{x}^{*}(t)]$|⁠, |$e_{p}^{*} \in M_{h}[\mathbf{x}^{*}(t)]$| and |$ e_w^* \in e_u^* + \mathring{X}_h[\mathbf{x}^*(t)] $| be the finite element error functions on |$\varOmega _{h}^{*}(t)$| with nodal vectors |$\mathbf{e}_{\mathbf{x}}$|⁠, |$\mathbf{e}_{\mathbf{u}}$|⁠, |$\mathbf{e}_{\mathbf{p}}$| and |$\mathbf{e}_{\mathbf{w}}$|⁠, respectively.

Let |$t^{*} \in (0,T]$| be the maximal time such that the following inequalities hold for |$t \in [0, t^{*}]$|⁠:

A key issue in the proof is to ensure that the |$W^{1, \infty }$| norm of the position error |$e_{x}^{*}(t)$| remains small. Since (4.22a) is small only for |$k \geqslant 2$|⁠, the estimates in this proof cannot be derived for the linear case |$k = 1$|⁠. In fact, |$t^{*}$| is positive because the inequalities above hold for |$t = 0$| when |$h$| is smaller than some sufficiently small constant.

• (1) Since |$\mathbf{x}^{*}(0) = \mathbf{x}^{0}$| and |$\mathbf{u}^{*}(0) = \mathbf{u}^{0}$|⁠, it follows that |$\mathbf{e}_{\mathbf{x}}(0) = 0$| and |$\mathbf{e}_{\mathbf{u}}(0) = 0$|⁠, which gives |$e_{x}^{*}(0) = 0$| and |$e_{u}^{*}(0) = 0$|⁠, respectively.

• (2) Testing (3.27c) with |$\mathbf{e}_{\mathbf{w}}-\mathbf{e}_{\mathbf{u}}$| at |$t=0$| and using the fact that |$\mathbf{e}_{\mathbf{x}}(0)=\mathbf{e}_{\mathbf{u}}(0)=0$| yields the following relation:

  $$ \begin{align*} \|\mathbf{e}_{\mathbf{w}}(0) \|_{\mathbf{A}(\mathbf{x}^{*}(0))}^{2} = - \mathbf{e}_{\mathbf{w}}(0)^\top \boldsymbol{{\mathcal{E}}}_{\mathbf{w}}(0). \end{align*} $$

The relation above, together with Lemma 3.1, implies that |$\|\mathbf{e}_{\mathbf{w}}(0) \|_{\mathbf{K}(\mathbf{x}^{*})} \leqslant Ch^{k}$|⁠. By the inverse inequality of finite element function, we have

Since |$k \geqslant 2> d/2$|⁠, it follows that (4.22b)–(4.22c) hold for |$t = 0$| when |$h$| is sufficiently small.

We first prove that error estimate (4.3) holds for |$t\in [0,t^{*}]$| under condition (4.22). Then we complete the proof by showing that |$t^{*}=T$|⁠.

(A) Estimates for |$\mathbf{e}_{\mathbf{x}}$|: Multiplying (3.27d) by matrix |$\mathbf{A}(\mathbf{x}^{*})$|⁠, testing with |$\mathbf{e}_{\mathbf{x}}$| and dropping the omnipresent argument |$ t \in [0, t^{*}]$|⁠, we have

In order to apply Gronwall’s inequality, we relate |$\frac{{\mathrm d}}{{\mathrm d} t}\| \mathbf{e}_{\mathbf{x}}\|_{\mathbf{A}(\mathbf{x}^{*})}^{2}$| to |$(\mathbf{e}_{\mathbf{x}})^\top \mathbf{A}(\mathbf{x}^{*})\dot{\mathbf{e}}_{\mathbf{x}}$| as follows:

By using Lemma 4.5, we have

(B) Estimates for |$\mathbf{e}_{\mathbf{w}}$| under condition (4.22): Since |$e_{w}^{*}- e_{u}^{*} = 0$| on |$\varGamma _{h}^{*}(t)$|⁠, i.e., |$e_{w}^{*}- e_{u}^{*}\in \overset{ \circ }{X}_{h}[\mathbf{x}^{*}(t)]$|⁠, we can choose |${\boldsymbol \chi }=\mathbf{e}_{\mathbf{w}} - \mathbf{e}_{\mathbf{u}}$| in (3.27c). This yields the following result:

The first term on the right-hand side of (4.26) can be estimated by using Lemma 4.2 and Hölder’s inequality, i.e.,

where we have used |$ \| \nabla w_{h}^\theta \|_{L^\infty (\varOmega _{h}^\theta (t))} = \| \nabla e_{w}^\theta - \nabla w_{h}^{*, \theta } \|_{L^\infty (\varOmega _{h}^\theta (t))} \leqslant C \big (1 + \|\nabla e_{w}^{*}\|_{L^\infty (\varOmega _{h}^{*}(t))}\big )$|⁠.

The second term on the right-hand side of (4.26) can be bounded by

The last term on the right-hand side of (4.26) can be estimated similarly, i.e.,

where we have converted |$ \|\mathbf{e}_{\mathbf{w}} - \mathbf{e}_{\mathbf{u}}\|_{\mathbf{K} (\mathbf{x}^{*})} $| (i.e., the full |$H^{1}$| norm) to |$ \|\mathbf{e}_{\mathbf{w}} - \mathbf{e}_{\mathbf{u}}\|_{\mathbf{A} (\mathbf{x}^{*})} $| (i.e., the |$H^{1}$| semi-norm) because the corresponding function |$e_{w}^{*}-e_{u}^{*}$| satisfies the zero boundary condition in |$\varOmega _{h,\pm }^{*}(t)$|⁠. Altogether we obtain the bound as

Then, substituting (4.28) into (4.25), we obtain

By integrating (4.29) on both sides with respect to time from |$0$| to |$t$|⁠, we obtain

(C) Estimates for |$\mathbf{e}_{\mathbf{u}}$| under condition (4.22): Testing (3.27a) by |$\dot{\mathbf{e}}_{\mathbf{u}}$|⁠, we obtain

The last term on the left-hand side of (4.31) can be estimated by using Lemma 4.5, i.e.,

In the following, we present estimates for each term on the right-hand side of (4.31).

The first term on the right-hand side of (4.31) can be estimated by using Lemma 4.2 and Hölder’s inequality, and the norm equivalence in Lemma 4.3, i.e.,

where the positive constant |$\varepsilon $| can be arbitrarily small and we have used |$\| \partial _{t, *}^{\bullet } u_{h}^{*}\|_{L^\infty (\varOmega _{h}^{*}(t))} \leqslant C$|⁠.

The second term on the right-hand side of (4.31) can be estimated by using integration by parts in time, i.e.,

where the first term on the right-hand side of (4.34) can be bounded after integrating it over |$[0, t]$| for |$0 \leqslant t \leqslant t^{*}$|⁠, i.e.,

In order to estimate |$S_{22}$|⁠, we define |$w_{h}^\theta = w_{h}^{*, \theta } + \theta e_{w}^{\theta }$| to be the finite element function on |$\varOmega _{h}^\theta (t)$| with nodal vector |$ \mathbf{w}^\theta = \mathbf{w}^{*} + \theta \mathbf{e}_{\mathbf{w}}$|⁠. Note that |$\frac{{\mathrm d}}{{\mathrm d} t} \mathbf{x}^\theta = \dot{\mathbf{x}} = \mathbf{w}^\theta $| and |$\frac{{\mathrm d}}{{\mathrm d} \theta } \mathbf{x}^\theta = \mathbf{e}_{\mathbf{x}}$|⁠. By using the transport formula, Lemma 4.1 and (4.5), we have

where |$\partial _{\theta }^{\bullet }$| denotes the material derivative with respect to the velocity field |$e_{x}^{\theta }$| with

and |$\partial _{t, \theta }^{\bullet }$| denotes the material derivative with respect to the velocity field |$w_{h}^{\theta }$| with

where we note |$\partial _{t, \theta }^{\bullet } e_{x}^\theta = e_{w}^\theta $| and |$\eta _{h}^\theta $| is defined from (4.7), which is the finite element function on |$\varOmega _{h}^\theta (t)$| and shares the nodal vector |${\boldsymbol{\eta }} \in{{\mathbb{R}}}^{dM}$| with |$ \dot{\boldsymbol{\eta }} = 0$|⁠.

The second term on the right-hand side of (4.34) can be estimated by using the above-mentioned setting and Lemma 4.2, we have

Similarly, as the estimates of (4.35), the third term on the right-hand side of (4.34) can be estimated as

The estimates of |$\int _{0}^{t} S_{21}(s)\, {\mathrm d} s$|⁠, |$S_{22}$| and |$S_{23}$| show that

The third term on the right-hand side of (4.31) can be written as

where |$S_{31}$| and |$S_{32}$| can be treated similarly as |$S_{21}$| and |$S_{22}$|⁠, with the following estimates:

where we have used the results |$\|e_{x}^{*}\|_{W^{1,\infty }(\varOmega _{h}^{*}(t))} \leqslant h^{-\frac{d}{4}}h^{\frac{k}{2}}$| and |$\|e_{u}^{*}\|_{W^{1,\infty }(\varOmega _{h}^{*}(t))}\leqslant 1$| in (4.22). Therefore,

The fourth term on the right-hand side of (4.31) can be estimated similarly, i.e.,

The fifth term on the right-hand side of (4.31) can be written as

where the first term on the right-hand side of (4.47) can be bounded after integrating it over |$[0, t]$| for |$0 \leqslant t \leqslant t^{*}$|⁠, i.e.,

The second and third terms on the right-hand side of (4.47) can be treated similarly as |$S_{22}$| and |$S_{23}$|⁠, with the following estimates:

The sixth term on the right-hand side of (4.31) can be estimated by differentiating (3.27b) in time. This yields the following relation:

Testing this equation with |$\mathbf{e}_{\mathbf{p}}$| and using relation |$\mathbf{q}^\top \mathbf{C}(\mathbf{x})\mathbf{v} = (\mathbf{C}(\mathbf{x})\mathbf{v})^\top \mathbf{q}$|⁠, we have

which can be estimated similarly as |$S_{22}$|⁠, |$S_{23}$| and so on. In particular, we have

The seventh term on the right-hand side of (4.31) can be written as

By using the approximation properties of the Lagrange interpolations |$ f_{h} $| and |$ f_{h}^{*} $| for |$ f $| on the domains |$ \varOmega _{h}(t) $| and |$ \varOmega _{h}^{*}(t) $|⁠, respectively, the terms |$ S_{71} $| and |$ S_{72} $| can be estimates as

By applying |$ \partial _\theta ^{\bullet } f = \nabla f \cdot e_{x}^{\theta }$| on |$ \varOmega _{h}^\theta (t) $|⁠, the term |$ S_{73} $| can be bounded similarly to the estimates for (4.33) as

The last term on the right-hand side of (4.31) can also be estimated by using integration by parts in time, i.e.,

Then, using the consistency estimates in Proposition 3.1, we have

Integrating the inequality (4.31) in time from |$0$| to |$t$| and using the estimates of |$S_{j}$|⁠, |$j=1,\dots ,6$|⁠, as well as the estimates in (4.58), with a sufficiently small |$\varepsilon $| in these estimates, the terms |$\varepsilon \|\mathbf{e}_{\mathbf{u}}(t)\|_{\mathbf{A}_{\boldsymbol{\mu }} (\mathbf{x}^{*})}^{2}$| and |$\varepsilon \int _{0}^{t} \|\dot{\mathbf{e}}_{\mathbf{u}}(s)\|_{\mathbf{M}_{\boldsymbol{\rho }}(\mathbf{x})}^{2}\, {\mathrm d} s$| on the right-hand side can be absorbed by the left-hand side. For |$k \geqslant 2$| and sufficiently small |$h$|⁠, the term |$Ch^{\frac{k}{2}-\frac{d}{4}} \|\mathbf{e}_{\mathbf{u}}(t)\|_{\mathbf{A}_{\boldsymbol{\mu }} (\mathbf{x}^{*})}^{2}$| on the right-hand side can also be absorbed by the left-hand side. Then we obtain the following result:

The first term on the right-hand side of (4.59) is estimated in the next part of the proof.

(D) Estimates for |$\mathbf{e}_{\mathbf{p}}$| under condition (4.22): To estimate |$\|\mathbf{e}_{\mathbf{p}}\|_{\mathbf{M}(\mathbf{x}^{*})}^{2}$|⁠, we test (3.27a) by a vector |$\mathbf{v}$| associated to a finite element function |$v_{h}$| on |$\varOmega _{h}(t)$| to get

The estimates of the terms on the right-hand side of (4.60) are similar as that in part (C) of this proof and therefore omitted. Specifically, we have

By using the inf-sup condition (2.6) in the expression for |$\|\mathbf{e}_{\mathbf{p}}\|_{\mathbf{M}(\mathbf{x})}$|⁠, we can estimate |$\|\mathbf{e}_{\mathbf{p}}\|_{\mathbf{M}(\mathbf{x})}$| as follows:

which together with the norm equivalence in Remark 4.2 gives

By integrating the above inequality with respect to time from |$0$| to |$t$|⁠, we have

(E) Combination of the estimates: Summing up inequalities (4.30), (4.59) and |$\varepsilon \times $|(4.64), substituting estimate (4.28) into the resulted inequality and choosing a sufficiently small |$\varepsilon $|⁠, we obtain the following combined estimate:

By applying Gronwall’s inequality, we obtain

This result can be substituted into (4.28) and (4.64) to yield an error estimate for |$\mathbf{e}_{\mathbf{w}}$| and |$\mathbf{e}_{\mathbf{p}}$|⁠, i.e.,

It remains to show that |$t^{*} = T$| when |$h$| is sufficiently small (smaller than some constant). In fact, for |$0 \leqslant t \leqslant t^{*}$|⁠, we can use the inverse inequality and (4.66) to bound the left-hand side in (4.22a) as

where |$C$| is a constant dependent on |$T$|⁠, but independent of |$t^{*}$|⁠. Since |$k \geqslant 2$|⁠, we have the following inequality for |$0 \leqslant t \leqslant t^{*}$| when |$h$| is sufficiently small:

For the same reason, for |$k \geqslant 2$| and sufficiently small |$h$|⁠, we have

Hence, we can extend the bound (4.22a)–(4.22c) beyond |$t^{*}$|⁠, which contradicts the maximality of |$t^{*}$| unless |$t^{*} = T$|⁠. This proves that the error estimates in (4.66)–(4.67) actually hold for |$t^{*}=T$|⁠. Since the flow map |$\phi _{h}^{*}(\cdot , t)$| and its inverse are bounded in the |$W^{1,\infty }$| norm, it follows that error estimates can be pulled back to the initial domain |$\varOmega _{h}^{0}$| after composition with |$\phi _{h}^{*}(\cdot , t)$|⁠. This proves the desired error estimates in Theorem 4.1.

5. Numerical results

In this section, we present numerical results for a benchmark example without surface tension. One of the most well-known benchmark examples for the two-phase NS flow problem was proposed in (Hysing et al., 2009, Table I) for simulating a moving bubble in fluid. This example includes surface tension and has been widely adopted to test the convergence and robustness of numerical methods for two-phase NS flows; see Aland & Voigt (2012); Barrett et al. (2015b); Fu (2020). In our study, we modify this model by considering a simplified version without surface tension and with no-slip boundary conditions. All the computations are performed using the software package NGSolve (https://ngsolve.org), and the mesh is generated by Gmsh (https://gmsh.info/).

The computational domain is |$\varOmega = (0, 1)\times (0, 2)$| and the no-slip pieces of the boundary are given by

The initial interface is a circle with radius 0.25, i.e.,

The gravitational force is |$f = (0, - 0.98)^\top $| and the physical parameters are shown in Table 1.

We first present numerical results to illustrate the convergence rate of the proposed method (2.8) for BP-1 in Table 1. We apply the time discretization of the proposed method (2.8) with a linearly semi-implicit Euler method, which only requires solving several decoupled linear systems at every time level. Since the exact solution of the problem is unknown, we evaluate the convergence rate in space by the formula

based on the finest three meshes, where |$u_{\tau , h}$| denotes the numerical solution at |$t_{N}=T$| computed by using a stepsize |$\tau $| and mesh size |$h$|⁠.

The spatial discretization errors of the numerical solutions at |$T$| = 0.1 and 0.5 are presented in Figs 2 and 3, respectively, where we have used the following finite elements:

Fig. 2.

Spatial discretization errors of the numerical solutions at |$T = 0.1$|⁠.

Open in new tabDownload slide

Fig. 3.

Spatial discretization errors of the numerical solutions at |$T = 0.5$|⁠.

Open in new tabDownload slide

with sufficiently small time step size |$\tau = 10^{-3}$| so that the temporal discretization errors are negligibly small compared to the spatial errors. From Figs 2–3, we see that the error of space discretization is |$O(h^{k})$|⁠, which is consistent with the theoretical results proved in Theorem 4.1.

Next, we present simulations of the bubble’s circularity, center of mass, rise velocity and the total energy of the two-phase NS system with parameters BP-1 in Table 1. We focus on the following benchmark quantities in Hysing et al. (2009) with respect to time.

• (a) Center of mass:

  $$ \begin{align*} X_{d} = \frac{1}{|\varOmega_{-}|} \int_{\varOmega_{-}} x_{d} \, {\mathrm d} x. \end{align*} $$
• (b) Circularity (referred as sphericity in 3D): In two dimensions it is defined as the quotient between perimeter of area-equivalent circle and perimeter of the bubble. In three dimensions it is defined as the quotient between surface area of volume-equivalent ball and surface area of the bubble, i.e.,

  $$ \begin{align*} Cir = \frac{2\sqrt{|\varOmega_{-}| \pi}}{|\varGamma|} \quad(\mbox{in 2-D}) \qquad \mbox{and} \qquad Cir = \frac{ \pi^{1/3} |\varOmega_{-}|^{2/3}}{|\varGamma|} \quad(\mbox{in 3-D}). \end{align*} $$
• (c) Rise velocity:

  $$ \begin{align*} V_{d} = \frac{1}{|\varOmega_{-}|} \int_{\varOmega_{-}} u_{h} \cdot e_{d} \, {\mathrm d} x, \end{align*} $$

  where the integrand |$ u_{h} \cdot e_{d}$| refers to the vertical component of the velocity |$u_{h}$|⁠.

We consider the proposed method (2.8) with polynomial degree |$k = 2$| and quadratic isoparametric triangular elements on the two different meshes as depicted in Fig. 4, where the mesh |$M_{1}$| has mesh size |$h_{1} = 0.04$| and the mesh |$M_{2}$| has mesh size |$h_{2} = 0.02$|⁠. Zoom-in of the deformed mesh around |$\varOmega _{-}(t)$| at |$t = 1$|⁠, |$1.5$| and |$2.5$| are shown in Fig. 5.

Fig. 4.

Left: mesh |$M_{1}$|⁠. Right: mesh |$M_{2}$|⁠. Dark color: domain |$\varOmega _{-}(0)$|⁠.

Open in new tabDownload slide

Fig. 5.

Zoom-in of the deformed mesh around |$\varOmega _{-}(t)$| with |$ k=2 $|⁠. Top: |$M_{1}$|⁠. Bottom: |$M_{2}$|⁠. Left: |$t_{1}$| = 1. Middle: |$t_{2}$| = 1.5. Right: |$t_{3}$| = 2.5.

Open in new tabDownload slide

In the computation the mesh may degenerate and in this case one has to regenerate the mesh. In our computation we regenerate the mesh when the following angle criterion is satisfied:

where |$\angle K$| denotes the set of interior angles of a triangle |$K$|⁠. The remeshing is performed using the software Gmsh. If the mesh is regenerated then we take the new mesh as the initial mesh to compute the above benchmark quantities. Fig. 6 shows the number of remeshing steps over time for two different meshes, using a fixed stepsize |$\tau = 1/200$|⁠. The meshes before and after the 3rd and the 6th remeshing steps are shown in Figs 7–8. In Fig. 9, we present the numerical results of circularity, center of mass, rise velocity and energy on the two different meshes with |$\tau = 1/200$|⁠. From Figs 4–9 we see that the method with the coarse mesh captures the shape of the interface very well. The geometry of the interface at |$t = 1$|⁠, |$1.5$| and |$2.5$|⁠, as well as the distribution of interface-velocity with the mesh |$M_{2}$|⁠, are presented in Fig. 10.

Fig. 6.

The number of remeshing steps with respect to time.

Open in new tabDownload slide

Fig. 7.

Meshes before and after a remeshing step for mesh |$M_{1}$| with |$k=2$|⁠.

Open in new tabDownload slide

Fig. 8.

Meshes before and after a remeshing step for mesh |$M_{2}$| with |$k=2$|⁠.

Open in new tabDownload slide

Fig. 9.

Benchmark quantities with respect to time.

Open in new tabDownload slide

Fig. 10.

Velocity. Left: |$t_{1}$| = 1. Middle: |$t_{2}$| = 1.5. Right: |$t_{3}$| = 2.5.

Open in new tabDownload slide

In summary, with the numerical comparison of results between coarse and fine meshes in Figs 4–9, we have good reasons to trust the numerical simulation results for this benchmark example.

Acknowledgements

All authors contributed equally.

Funding

National Natural Science Foundation of China (grant No. 12201418 to B.L.) and a fellowship award from the Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. PolyU/RFS2324-5S03 to B.L.); Research Grants Council of the Hong Kong Special Administrative Region, China (Project No. CityU 11300621 to S.M. and W.Q.).

References