A stable projection method for the incompressible Navier–Stokes equations on arbitrary geometries and adaptive Quad/Octrees

Abstract

We present a numerical method for solving the incompressible Navier–Stokes equations on non-graded quadtree and octree meshes and arbitrary geometries. The viscosity is treated implicitly through a finite volume approach based on Voronoi partitions, while the convective term is discretized with a semi-Lagrangian scheme, thus relaxing the restrictions on the time step. A novel stable implementation of the projection step is introduced, making use of the Marker And Cell layout for the data. The solver is validated numerically in two and three spatial dimensions.

Introduction

The prediction of fluid motion around structures is crucial in many important applications in science and engineering. Examples include the classical study of the aerodynamics of aircrafts or the hydrodynamics of ships, but also more modern applications such as the fluid dynamics occurring during materials processing such as the solidification of liquid metal alloys, the study of artificial swimmers or biological flows. One of the difficulties of solving the equations of fluid dynamics is in dealing with non-trivial geometries, which can be explicitly described (body-fitted approaches) or implicitly captured. We focus here on strategies on Cartesian grids, where the geometry is implicitly captured and refer the interested reader to the book by Peric and Ferziger [23] for discussions of body-fitted approaches.

Much of the early work was concerned with compressible flows and strategies based on Cartesian grids were first introduced on uniform grids by Purvis and Burkhalter [55], who used a finite volume approach to solve the two-dimensional potential equation. Later, researchers proposed solvers for the Euler equations in two [24], [28] and three [24] spatial dimensions. One of the difficulties inherent to fluid flows is their spatial multiscale nature. From boundary layers close to solid boundaries to the generation of vorticity and turbulence, specific regions require a very fine level of detail while other areas can be treated adequately with a coarser mesh. Researchers have proposed strategies to alleviate this problem by designing numerical methods on spatially adaptive grids, which include stretched grids (see e.g. [62], [19]), nested grids (see e.g. [8], [31], [60], [7]), or unstructured meshes (see e.g. [65], [33], [44], [42], [43]). Large parallel codes have also been written and used in commercial applications, e.g. TRANAIR, which is an adaptive Cartesian full potential solver coupled with a viscous boundary layer model [64], [15] or NASA Cart3D [1], which is a solver for the compressible Euler equations, with application to high-speed flows.

We are focusing in this paper on the incompressible Navier–Stokes equations on Octree data structures, which provide the ability to refine/coarsen continuously in space.¹ These approaches follow the general framework of the projection method of Chorin [13], which leverages the Hodge decomposition of vector fields: first an intermediate velocity field is computed, before applying a projection onto the divergence-free subspace (the interested reader is referred to the excellent paper by Brown, Cortez and Minion [10] for a review of different projection methods). In that vein, Popinet [54] introduced a solver on Octrees using finite volume discretizations on the Marker And Cell (MAC) configuration [29]. In this work, the size between adjacent cells is constrained by a 2:1 ratio, which reduces the number of local grid configurations. In turn, this can be exploited to construct second-order approximations of the projection step, albeit leading to a non-symmetric linear system. Later, Losasso et al.  introduced a solver for the incompressible Navier–Stokes equations and for free surface flows [39]. This approach considers octrees for which the ratio between adjacent cells is not constrained, allowing for increased adaptivity of the computational mesh. The projection step uses a finite volume approach and leads to a symmetric linear system. However, the Hodge variable is only first-order accurate, which impacts on the accuracy of the velocity field. Losasso et al. [38] also introduced a second-order solver for the Poisson equation, based on the work of Lipnikov [36], but did not use that solver for fluid simulations. Later, Min and Gibou introduced a solver that uses an approach based on finite differences instead of finite volumes [45]. The linear system for the Hodge variable is non-symmetric but the solution is second-order accurate with second-order accurate gradients, which in turn produces also second-order accurate velocity fields.

However, one of the main challenges when considering a finite difference approach on adaptive meshes is the potential loss of numerical stability. In [45], Min and Gibou showed that the standard projection method cannot be guaranteed to be stable in that framework; it was confirmed numerically and shown to be exacerbated by high size ratios between adjacent cells and high Reynolds numbers. They also introduced the so-called “orthogonal projection” method that guarantees numerical stability (even though their method is not conservative) in the case where Dirichlet boundary conditions for the velocity field are imposed. However, the numerical stability is not guaranteed for inflow/outflow boundary conditions, which limits the range of applications of that approach. An approach based on the MAC grid configuration is more amenable to designing stable projection solvers: following the standard proof of $L^2$-stability, one can show that a minus transpose relationship of the discrete gradient and divergence operators is enough to guarantee numerical stability in a weighted $L^2$-norm. Such a constraint can be enforced in an MAC grid sampling of the data, even if no constraint is imposed on the grid; a difference from a node-based approach. In this paper, we present a projection method that is stable, using the Poisson solver introduced in [38] and deriving the numerical approximations of gradient and divergence operators to ensure numerical stability. An MAC grid discretization also has the desirable property of being conservative.

Another challenge is the representation of the interface between the fluid and the irregular boundaries, and more specifically how to impose the boundary conditions in an implicit framework. A common approach is to use Peskin's immersed boundary method [53], [20]. However, we seek to avoid the smoothing of the solution induced by a delta formulation and the subsequent decrease in accuracy in the $L^{\infty }$-norm near the boundaries; thus we represent the location of the interface implicitly with a level-set function [52] and impose the boundary conditions sharply on the interface. Several strategies have been introduced, e.g. the rasterization approach used in Losasso et al. [39] or the Heaviside formulation of Batty et al. [6]. However, Ng et al.  showed that treatments such as these lead to a method that does not converge in the $L^{\infty }$-norm [51] while a finite-volume/cut-cell approach in a level-set framework produces accurate results. We use that approach.

The last challenge in solving the incompressible Navier–Stokes equations is the time step restriction usually resulting from the discretizations of the advection (CFL condition) and the viscous ($\mathrm{\Delta }t=O(\mathrm{\Delta }{x}^2)$) terms. We circumvent those issues by employing a BDF semi-Lagrangian scheme [63] for the advection term and by treating the viscous term implicitly. Since we opt for an MAC layout, producing a compact accurate implicit solver for the viscous term is not straightforward, but it can be done with a finite volume approach where the control volumes are the elements of a Voronoi partition, as discussed in [30].