A multigrid ghost-point level-set method for incompressible Navier-Stokes equations on moving domains with curved boundaries

Highlights

• A finite-difference ghost-point method to solve Navier-Stokes equations in complex-shaped moving domains is proposed.

• The method is second order accurate.

• The accuracy is not affected by moving domains or complex-shaped objects.

• A proper geometric multigrid approach is presented, whose efficiency is not affected by complex-shaped moving objects.

Abstract

In this paper we present a numerical approach to solve the Navier-Stokes equations on moving domains with second-order accuracy. The space discretization is based on the ghost-point method, which falls under the category of unfitted boundary methods, since the mesh does not adapt to the moving boundary. The equations are advanced in time by using Crank-Nicholson. The momentum and continuity equations are solved simultaneously for the velocity and the pressure by adopting a proper multigrid approach. To avoid the checkerboard instability for the pressure, a staggered grid is adopted, where velocities are defined at the sides of the cell and the pressure is defined at the centre. The lack of uniqueness for the pressure is circumvented by the inclusion of an additional scalar unknown, representing the average divergence of the velocity, and an additional equation to set the average pressure to zero. Several tests are performed to simulate the motion of an incompressible fluid around a moving object, as well as the lid-driven cavity tests around steady objects. The object is implicitly defined by a level-set approach, that allows a natural computation of geometrical properties such as distance from the boundary, normal directions and curvature. Different shapes are tested: circle, ellipse and flower. Numerical results show the second order accuracy for the velocity and the divergence (that decays to zero with second order) and the efficiency of the multigrid, that is comparable with the tests available in literature for rectangular domains without objects, showing that the presence of a complex-shaped object does not degrade the performance.

Introduction

Navier-Stokes equations are a system of partial differential equations (PDEs) describing the motion of an incompressible viscous fluid subject to external forces. In the non-dimensional form, they read

$$\begin{gathered} \frac{\partial \mathbf{u}}{\partial t}+\mathbf{u}\cdot \mathrm{\nabla }\mathbf{u}+\mathrm{\nabla }p=\frac{1}{\text{Re}}{\mathrm{\nabla }}^2\mathbf{u}+\mathbf{f} \\ \mathrm{\nabla }\cdot \mathbf{u}=0 \end{gathered}$$

where u is the fluid velocity, p is the pressure, Re is the Reynolds number and f the external force. The first equation is obtained from the momentum equation, while the second one is the continuity equation and represents the incompressibility condition. Navier-Stokes equations have been central to innumerable scientific and engineering applications for several decades, as the necessity to model the motion of incompressible fluids has a tremendous impact on the real-life world. Several numerical approaches have been developed so far to provide efficient and accurate methods in different contexts, both for time and space discretization. An overview on traditional methods can be found in [43].

One of the first approaches for time discretization was the projection method, proposed by Chorin in [14]. It consists of a fractional step, where the velocity u is advanced in time by discretising the momentum equation and ignoring the pressure term, obtaining an intermediate velocity ${\mathbf{u}}^{⁎}$. This velocity is not divergence-free, and therefore is corrected by the pressure term, that is determined by enforcing the divergence-free condition. Projection methods have been very popular for Navier-Stokes equations and are widely adopted nowadays, as they are computationally efficient and easy to implement (they require the solution of a sequence of elliptic equations for the intermediate velocity and for the pressure). However, the definition of boundary conditions for ${\mathbf{u}}^{⁎}$ and p is not obvious and it is challenging to achieve an accuracy higher than first order, even if the internal equations are discretized with higher order. There are existing approaches that suggest appropriate boundary conditions to improve the accuracy, as described in [11]. However, these boundary conditions are fully coupled each other and the resulting elliptic equations for ${\mathbf{u}}^{⁎}$ and for p cannot be solved in sequence, as the boundary condition on p depends on ${\mathbf{u}}^{⁎}$.

Another approach consists of discretising both PDEs of Eq. (1) and solve them simultaneously in the unknowns $(\mathbf{u},p)$. This is the monolithic approach and has the advantage that no boundary conditions are required for the pressure, as described, for example, in [61, Sect. 8.1]. These methods can usually achieve higher order accuracy than projection methods, but they have the disadvantage that extra computational effort is needed to solve both PDEs simultaneously. This aspect can be alleviated if efficient approaches are adopted such as the monolithic solver for fluid/solid interactions proposed in [35], where the authors achieve second order accuracy by maintaining a symmetric positive definite discretization.

Projection and monolithic methods refer only to time discretization and then spatial discretization methods must be addressed regardless of the choice between the methods described above. Standard problems on rectangular domains are usually discretized in space by using a finite volume approach on a staggered grid (MAC grid, see for example [59]), where the velocities are defined at the edges of the cell and the pressure is defined at the centre. The MAC grid is used for example to avoid the checkerboard instability for the pressure term observed in non-staggered grids due to the fact that p appears in the equations only in the form of ∇p. The extension of MAC grids to non-rectangular domains (curved boundaries) is not straightforward, since a particular treatment is needed close to the boundary. Moreover, there is an increasing scientific interest to Navier-Stokes equations for complex-shaped (possibly moving) domains, for example to model the motion of incompressible fluids around a deforming object, that may be central to more sophisticated applications such as fluid-structure interactions. Numerical methods for steady complex-shaped objects are mainly divided between fitted and unfitted boundary methods. In fitted boundary methods, the grid, either structured or unstructured (such as Finite Element Methods (FEM)), is adapted to the curved boundary [58], [38], [36], [47], [46], [25], [24]. Although they provide an extremely accurate representation of the domain, the generation of the mesh may become cumbersome and significantly expensive from a computational point of view, especially in moving domain problems where a new mesh generation is needed at each time step. In unfitted boundary methods the moving domain is instead embedded in a fixed grid. The methods of this class are based either on FEM, such as XFEM [40], [55] or GFEM [15], on meshless methods [62], or Cartesian meshes, such as cut cell methods [37].

An alternative approach of finite–difference methods for curved boundaries is represented by the Ghost–Point method. A pioneering work of this class of methods is the Immersed Boundary Method [54] to model blood flows in the heart, where Peskin presented a first-order accurate method, later extended to higher order by the Immersed Interface Methods proposed by LeVeque and Li in [44]. More recent finite-difference methods to discretize PDEs on complex domains are the Ghost-Fluid Methods proposed in [30], [33], [34], [50], where the grid points external to the domain are called ghost points and a fictitious (ghost) value of the numerical solution is formally extrapolated onto these grid points in order to maintain a stable algorithm and achieve the desired accuracy order. The main feature of these methods is that the internal equations are firstly formally solved for the ghost values, resulting in a final linear system with eliminated boundary conditions. This class of methods may fail to obtain high accuracy for Neumann boundary condition, and anyway it remains first order accurate in the post-processing reconstruction of the gradient of the solution.

A new finite-difference method with an improved accuracy was proposed by the author and Russo in [20], [22], [19] for elliptic equations and extended to thermo-poroelastic Cauchy-Navier equations for elastic deformation [17], Euler equations for gas dynamics [13], [21], electron transport in semiconductors [48], [49], fluid flow in porous media for water-rock interaction [18], [16]. In this approach, the ghost point values are eventually coupled each other, resulting in a bigger linear system with non-eliminated boundary conditions, solved by a multigrid approach suitably designed for ghost-point methods. The domain is implicitly described by a level-set function.

Other numerical methods to solve sharp-edge boundary problems are the finite volume methods, the non-symmetric positive definite finite element method, the matched interface and boundary (MIB) method [63], the arbitrary Lagrangian Eulerian method (ALE) [31], [27], the penalization methods [42], [5], [12], and the class of Immersed Finite Volume Methods (IFVM) [29]. Other recent advances have been obtained in [41] to achieve higher accuracy in the presence of Neumann boundary conditions, in [9] to solve elliptic interface problems with discontinuous coefficients, in [28] to solve sharp interface conditions and in [8], [6], [53] to obtain higher accuracy in the gradient of the solution.

Time and space discretizations of Navier-Stokes equations lead to a linear system that must be solved efficiently. Multigrid method is one of the most efficient iterative approaches to solve the linear systems arising from the discretization of a class of PDEs. Initially designed for elliptic problems, multigrid methods have been widely adopted to solve countless problems [61], [39], [10]. Several approaches have been proposed for boundaries aligned with grid lines. As example, we mention [4], [26], [3], [45], [1]. These approaches fall under the class of geometric multigrid, since the operators are defined from the geometrical background of the problem. An extension of geometric multigrid is represented by the class of algebraic multigrid methods [60], where the operators are defined merely from the information gathered from the linear systems, and therefore it is more practical for complex problems such as curved boundaries. However, the inclusion of geometrical information to the definition of the operators usually results in a more efficient solver, as highlighted for example in [2], where a comparison between geometric and algebraic multigrid methods for curved boundaries is performed.

In this paper, we provide a finite-difference method to solve the Navier-Stokes equations on moving domains following a monolithic approach for the time discretization and a finite-difference ghost-point method for the space discretization. The singularity of the problem (due to the fact that p is not uniquely defined, since it appears only in terms of ∇p) is circumvented by an additional equation for the pressure that requires that its average is zero. This additional equation is balanced with an additional unknown ξ that measures the average divergence of the velocity (that is therefore not necessarily zero, as in traditional methods, but decays with the same order of accuracy of the method). The linear system arising from the discretization is solved by a geometric multigrid approach, where the relaxation, restriction and interpolation operators are modified close to the boundary in order to maintain the same efficiency that standard multigrid techniques achieve away from the boundary.

The paper is structured as follows. In Sect. 2 and 3 we describe the time and space discretizations, respectively. The multigrid technique is detailed in Sect. 4, while numerical tests are provided in Sect. 5 and conclusions in Sect. 6.