On Meshfree GFDM Solvers for the Incompressible Navier-Stokes Equations

TL;DR

This paper critically examines meshfree GFDM schemes for incompressible Navier-Stokes equations, identifies their limitations in local mass conservation, and proposes a modified scheme that improves accuracy and divergence approximation.

Contribution

It introduces a novel modification to a monolithic GFDM scheme that enhances local mass conservation and velocity divergence approximation.

Findings

The new scheme outperforms existing methods in accuracy.

Identifies key drawbacks of conventional meshfree schemes.

Provides numerical evidence of improved divergence approximation.

Abstract

Meshfree solution schemes for the incompressible Navier--Stokes equations are usually based on algorithms commonly used in finite volume methods, such as projection methods, SIMPLE and PISO algorithms. However, drawbacks of these algorithms that are specific to meshfree methods have often been overlooked. In this paper, we study the drawbacks of conventionally used meshfree Generalized Finite Difference Method~(GFDM) schemes for Lagrangian incompressible Navier-Stokes equations, both operator splitting schemes and monolithic schemes. The major drawback of most of these schemes is inaccurate local approximations to the mass conservation condition. Further, we propose a new modification of a commonly used monolithic scheme that overcomes these problems and shows a better approximation for the velocity divergence condition. We then perform a numerical comparison which shows the new monolithic scheme to be more accurate than existing schemes.

Tables

Table 1. Table 1: Errors, convergence orders and simulation times for the Taylor-Green vortices test case. h ℎ h is the smoothing length, N 𝑁 N is the number of points in the entire domain at the initial state, ϵ 2 subscript italic-ϵ 2 \epsilon_{2} is the relative error, r 𝑟 r is the order of convergence of ϵ 2 subscript italic-ϵ 2 \epsilon_{2} , and t 𝑡 t is the simulation time in seconds.

		Projection			Coupled:Penalty			Coupled:New
$h$	$N$	${ϵ}_2$	$r$	$t$	${ϵ}_2$	$r$	$t$	${ϵ}_2$	$r$	$t$
$1.0$	$293$	$3.1\times {10}^{-2}$	$-$	$6$	$2.4\times {10}^{-2}$	$-$	$7$	$1.2\times {10}^{-2}$	$-$	$9$
$0.5$	$\mathrm{1\hspace{0.17em}047}$	$9.7\times {10}^{-3}$	$1.67$	$28$	$6.5\times {10}^{-3}$	$1.88$	$29$	$3.1\times {10}^{-3}$	$1.95$	$32$
$0.25$	$\mathrm{3\hspace{0.17em}856}$	$3.2\times {10}^{-3}$	$1.59$	$195$	$2.5\times {10}^{-3}$	$1.37$	$202$	$1.0\times {10}^{-3}$	$1.63$	$207$
$0.125$	$\mathrm{14\hspace{0.17em}878}$	$1.1\times {10}^{-3}$	$1.54$	$1934$	$8.2\times {10}^{-4}$	$1.60$	$2010$	$3.3\times {10}^{-4}$	$1.59$	$2031$

Table 2. Table 2: Errors, convergence orders and simulation times for the bifurcated tube test case. h ℎ h is the smoothing length, N 𝑁 N is the number of points in the entire domain at the initial state, ϵ m ​ a ​ s ​ s subscript italic-ϵ 𝑚 𝑎 𝑠 𝑠 \epsilon_{mass} is the relative error in mass conservation, r 𝑟 r is the order of convergence of ϵ m ​ a ​ s ​ s subscript italic-ϵ 𝑚 𝑎 𝑠 𝑠 \epsilon_{mass} , and t 𝑡 t is the simulation time in seconds.

		Projection			Coupled:Penalty			Coupled:New
$h$	$N$	${ϵ}_{mass}$	$r$	$t$	${ϵ}_{mass}$	$r$	$t$	${ϵ}_{mass}$	$r$	$t$
$0.5$	$\mathrm{5\hspace{0.17em}805}$	$1.2\times {10}^{-1}$	$-$	$556$	$7.3\times {10}^{-2}$	$-$	$718$	$4.1\times {10}^{-3}$	$-$	$927$
$0.25$	$\mathrm{21\hspace{0.17em}294}$	$9.1\times {10}^{-2}$	$0.40$	$7378$	$5.6\times {10}^{-2}$	$0.38$	$8752$	$2.5\times {10}^{-3}$	$0.71$	$6056$
$0.125$	$\mathrm{81\hspace{0.17em}125}$	$7.1\times {10}^{-2}$	$0.36$	$76850$	$4.1\times {10}^{-2}$	$0.45$	$76936$	$1.7\times {10}^{-3}$	$0.56$	$53477$

Table 3. Table 3: Errors, convergence orders and simulation times for the sloshing test case. h ℎ h is the smoothing length, N 𝑁 N is the number of points in the entire domain at the initial state, ϵ V subscript italic-ϵ 𝑉 \epsilon_{V} is the error in volume conservation, r 𝑟 r is the order of convergence of ϵ V subscript italic-ϵ 𝑉 \epsilon_{V} , and t 𝑡 t is the simulation time in seconds.

		Projection			Coupled:Penalty			Coupled:New
$h$	$N$	${ϵ}_V$	$r$	$t$	${ϵ}_V$	$r$	$t$	${ϵ}_V$	$r$	$t$
$0.06$	$\mathrm{\hspace{0.17em}382}$	$9.59\times {10}^{-2}$	$-$	$10$	$7.09\times {10}^{-2}$	$-$	$12$	$5.73\times {10}^{-2}$	$-$	$15$
$0.03$	$\mathrm{1\hspace{0.17em}361}$	$3.27\times {10}^{-2}$	$1.56$	$51$	$2.26\times {10}^{-2}$	$1.65$	$59$	$1.47\times {10}^{-2}$	$1.96$	$68$
$0.015$	$\mathrm{5\hspace{0.17em}119}$	$2.42\times {10}^{-2}$	$0.43$	$723$	$1.67\times {10}^{-2}$	$0.44$	$602$	$1.01\times {10}^{-2}$	$0.54$	$746$
$0.0075$	$\mathrm{19\hspace{0.17em}704}$	$1.43\times {10}^{-2}$	$0.76$	$7871$	$1.23\times {10}^{-2}$	$0.44$	$6685$	$0.75\times {10}^{-2}$	$0.43$	$7268$

Table 4. Table 4: Errors in Taylor expansions measured by the functional J = ∑ j ∈ S i W j ​ ( e j u ) 2 𝐽 subscript 𝑗 subscript 𝑆 𝑖 subscript 𝑊 𝑗 superscript superscript subscript 𝑒 𝑗 𝑢 2 J=\sum_{j\in S_{i}}W_{j}(e_{j}^{u})^{2} . The table shows the mean value of J 𝐽 J across all interior points at the end of the first time-step for the Taylor-Green vortices test case.

$h$	Classical GFDM: Projection	Modified GFDM: New Coupled Solver
$1.0$	$0.2976$	$0.3057$
$0.5$	$0.09579$	$0.09709$
$0.25$	$0.03255$	$0.03345$
$0.125$	$0.01153$	$0.01470$

Equations

$${\partial }^{\ast }u({\vec{x}}_i)\approx {\tilde{\partial }}_i^{\ast }u=\sum _{j\in S_i}{c}_{ij}^{\ast }{u}_j,$$

$$e_j+u({\vec{x}}_j)=u({\vec{x}}_i)+\nabla u\cdot ({\vec{x}}_j-{\vec{x}}_i)+\frac{1}{2}({\vec{x}}_j-{\vec{x}}_i{)}^{T}D({\vec{x}}_j-{\vec{x}}_i).$$

$$\text{min }J=\sum _{j\in S_i}{W}_j{e}_j^2,$$

$$W_j=\exp \left(-\alpha \frac{\Vert {\vec{x}}_j-{\vec{x}}_i{\Vert }^2}{h^2}\right),$$

$$\left(\begin{array}[]{c}e_{1}\\ \vdots\\ e_{n}\end{array}\right)=\left(\begin{array}[]{cc}\delta x_{1}&\frac{1}{2}\delta x_{1}^{2}\\ \vdots&\vdots\\ \delta x_{n}&\frac{1}{2}\delta x_{n}^{2}\\ \end{array}\right)\left(\begin{array}[]{c}u_{x}\\ u_{xx}\end{array}\right)-\left(\begin{array}[]{c}u_{1}-u_{i}\\ \vdots\\ u_{n}-u_{i}\end{array}\right)\,.$$

$$\text{min }J=(M\vec{a}-\vec{b}{)}^{T}W(M\vec{a}-\vec{b}),$$

$$\vec{a}=[(M^{T}WM{)}^{-1}{M}^{T}W]\vec{b}.$$

$$\displaystyle u_{x}\Bigr{|}_{i}$$

$$\displaystyle u_{xx}\Bigr{|}_{i}$$

$$au+b{u}_x+c{u}_{xx}=d.$$

$$a{u}_i+b\sum _{j\in S_i}{c}_{ij}^x(u_j-u_i)+c\sum _{j\in S_i}{c}_{ij}^{xx}(u_j-u_i)=d,i=1,\dots ,N$$

$$\sum _{j\in S_i}{c}_{ij}^x{m}_j$$

$$\text{min }J$$

$$\displaystyle\left(\begin{array}[]{c}e_{1}\\ \vdots\\ e_{n}\\ e_{n+1}\end{array}\right)$$

$$\text{min }J$$

$$u_i=\sum _{j\in S_i}{\alpha }_{ij}{u}_j+{\beta }_{i}d.$$

$$(1-{\alpha }_{ii})u_i-\sum _{\begin{array}{c}j\in S_i\\ j\ne i\end{array}}{\alpha }_{ij}{u}_j={\beta }_{i}d.i=1,\dots ,N.$$

$$\frac{D\vec{x}}{Dt}$$

$$\nabla \cdot \vec{v}$$

$$\frac{D\vec{v}}{Dt}$$

$${\vec{x}}_i^{(n+1)}={\vec{x}}_i^{(n)}+{\vec{v}}_i^{(n)}\Delta t+\frac{{\vec{v}}_i^{(n)}-{\vec{v}}_i^{(n-1)}}{\Delta t}(\Delta t{)}^2,$$

$$\frac{{\vec{v}}^{\ast }-{\vec{v}}^{(n)}}{\Delta t}=\frac{\eta }{\rho }\Delta {\vec{v}}^{\ast }-\frac{1}{\rho }\nabla p^{\ast }+\vec{g},$$

$${\vec{v}}^{(n+1)}={\vec{v}}^{\ast }-\frac{\Delta t}{\rho }\nabla p_{corr}.$$

$$\frac{\Delta t}{\rho }\Delta p_{corr}=\nabla \cdot {\vec{v}}^{\ast }.$$

$$p^{(n+1)}=p^{\ast }+p_{corr}.$$

$$\nabla \cdot \nabla p_{corr}=\Delta p_{corr}.$$

$$\left(\begin{array}[]{ccc}I-\frac{\Delta t}{\rho}\eta\mathbf{C}^{\Delta}&&\frac{\Delta t}{\rho}\mathbf{C}^{x}\\[5.0pt] &I-\frac{\Delta t}{\rho}\eta\mathbf{C}^{\Delta}&\frac{\Delta t}{\rho}\mathbf{C}^{y}\\[5.0pt] \mathbf{C}^{x}&\mathbf{C}^{y}&-A\frac{\Delta t}{\rho}\\ \end{array}\right)\left(\begin{array}[]{c}\vec{U}^{(n+1)}\\ \vec{V}^{(n+1)}\\ \vec{P}_{corr}\end{array}\right)=\left(\begin{array}[]{c}\vec{U}^{(n)}-\frac{\Delta t}{\rho}\mathbf{C}^{x}\vec{P^{*}}+\Delta t\vec{G_{x}}\\[5.0pt] \vec{V}^{(n)}-\frac{\Delta t}{\rho}\mathbf{C}^{y}\vec{P^{*}}+\Delta t\vec{G_{y}}\\[5.0pt] 0\\ \end{array}\right)\,,$$

$$\nabla \cdot {\vec{v}}^{(n+1)}-A\frac{\Delta t}{\rho }\Delta p_{corr}=0.$$

$$\frac{{\vec{v}}^{(n+1)}-{\vec{v}}^{(n)}}{\Delta t}$$

$$\nabla \cdot {\vec{v}}^{(n+1)}$$

$$e_1+u^{(n+1)}-\frac{\eta \Delta t}{\rho }\Delta u^{(n+1)}+\frac{\Delta t}{\rho }{\partial }^x{p}_{corr}$$

$$e_2+v^{(n+1)}-\frac{\eta \Delta t}{\rho }\Delta v^{(n+1)}+\frac{\Delta t}{\rho }{\partial }^y{p}_{corr}$$

$$e_3+\nabla \cdot {\vec{v}}^{(n+1)}$$

$$e_4+\Delta p_{corr}+\rho \nabla \cdot \left[\left({\vec{v}}^{(n)}\cdot \nabla \right){\vec{v}}^{(n+1)}\right]=\frac{\rho }{\Delta t}\nabla \cdot {\vec{v}}^{(n)}-\Delta p^{\ast }+\rho \nabla \cdot \vec{g}.$$

Full text

††journal: Elsevier

1 Introduction

In hydrodynamic simulations, the computational domain is often very complex or is rapidly changing with time. For such simulation domains, the task of meshing can prove to be very challenging. Meshfree or meshless methods are becoming increasingly popular for such simulations to avoid the tasks of meshing and remeshing. These methods use the numerical basis of a set of arbitrarily distributed nodes without any underlying mesh to connect them. These nodes could either be mass carrying particles or numerical points. For each node, the only connectivity information required is a local set of neighbouring nodes over which approximations are carried out.

Meshfree Generalized Finite Difference Methods (GFDMs) are one such class of meshfree methods. Meshfree GFDMs are strong-form methods based on weighted least squares approximations. They have been widely used (for example, [13, 20, 36, 44]) and are referred under various names, including FPM, which stands for both the Finite Pointset Method [38] and the Finite Point Method [30]; the Kinetic Meshless Method (KMM) [32] and the Least Squares Kinetic Upwind Method (LSKUM) [14]. They provide several advantages over meshfree particle methods such as Smoothed Particle Hydrodynamics (SPH) [24]. They provide a framework to naturally incorporate boundary conditions which is not possible in SPH. They are more adaptive in terms of the spatial discretization since they use numerical approximation points instead of mass carrying particles. These advantages, however, come at the price of a lack of strict conservation.

Several mesh-based algorithms to solve the incompressible Navier–Stokes equations have been extended to meshfree methods. These include operator splitting methods such as the projection method [8], SIMPLE [31] and PISO [18]. However, drawbacks of these algorithms that are specific to meshfree methods have often been overlooked. One important example of such a drawback is that the discrete Laplace operator is not the same as the discrete divergence of the discrete gradient operator. This results in inaccurate approximations to the mass conservation equation. While this problem has been solved in the context of Finite Volume Methods (FVMs) with staggered grids, this problem persists in meshfree methods and introduces errors in most operator splitting algorithms.

In this paper, we present a few commonly used meshfree solvers, both operator splitting solvers and monolithic solvers, for the incompressible Navier–Stokes equations and study the drawbacks of these algorithms specific to meshfree GFDMs. We then propose a new monolithic solver which attempts to overcome these drawbacks. This is done by solving an over-determined problem. The mass and momentum conservation equations are solved together, and simultaneously with a pressure-Poisson equation which is needed to improve stability. All the methods considered are spatially second order accurate and use first order time integrations. However, the new monolithic method proposed here is numerically shown to be much more accurate than existing methods. This is as a consequence of the new method providing better local approximations to the mass conservation condition.

In our earlier paper [37], we presented a method to improve global conservative properties in meshfree GFDMs, with an application to the incompressible Navier–Stokes equations. That was done by introducing an approximate discrete divergence theorem, without changing the local accuracy of the conservation equations. In contrast, in this paper, we present a method to improve the local accuracy of the mass conservation equation.

The paper is organized as follows. In Section 2, we give an overview of two variants of meshfree GFDMs. In Section 3, we present existing meshfree GFDM algorithms used to numerically solve the incompressible Navier–Stokes equations and analyze their drawbacks. A new coupled velocity-pressure solver which addresses these drawbacks is presented in Section 4. Numerical results and comparisons of the new method with existing methods are presented in Section 5, and the paper is concluded with a discussion in Section 6.

2 Meshfree Finite Differences

In meshfree GFDMs, a cloud of $N$ numerical points is used to discretize the computational domain $\Omega$. These points, with positions $\vec{x}_{i}$, $i=1,\dots,N$, are usually irregularly spaced and include both points in the interior and on the boundary $\partial\Omega$ of the domain. The only connectivity information required for each point $i$ is a set of neighbouring points $S_{i}$. In most cases, this neighbourhood or support $S_{i}$ consists of the $n=n(i)$ points closest to it, $S_{i}=\{\vec{x}_{j}:\|\vec{x}_{j}-\vec{x}_{i}\|\leq h\}$, where $h=h(\vec{x},t)$ is the smoothing length, which determines the spatial discretization size. In addition to the smoothing length $h$, the point cloud is also described by the parameters $r_{max}$ and $r_{min}$. Initial set up of the point cloud, and local addition and deletion of points during a simulation is done such that there exists at least one point in every sphere of radius $r_{max}h$, and that no two points are closer than $r_{min}h$. Further details about point cloud organization and management can be found in [10, 19, 34].

There are two variations of meshfree finite differences that we consider in this paper, both of which are based on moving least squares approximations. The first is the widely used approach based on the work of Liszka and Orkisz [23], which we refer to as the classical GFDM. This method very closely resembles traditional finite differences, and is obtained by minimizing errors obtained from Taylor expansions. The second approach is a modification of the classical GFDM based on the work of Tiwari and Kuhnert [38] in which the errors in the PDE considered is minimized simultaneously with the errors obtained from Taylor expansions. A comparison between these approaches has been done by Illiev and Tiwari [17]. In both variations, we consider Taylor expansions up to second order terms. Higher order accuracy can be attained in both formulations and have been considered by Milewski [27].

2.1 Classical meshfree GFDM

This classical method is based on the work of Liszka and Orkisz [23]. For a function $u$ defined at each numerical point $i=1,2,\dots,N$, its derivatives are approximated as

[TABLE]

where ${}^{*}=x,y,xx,\Delta,etc.$ represents the differential operator being approximated, $\partial^{*}$ represents the continuous ∗-derivative, and $\tilde{\partial}^{*}_{i}$ represents the discrete derivative at point $i$. For a point $i$, consider Taylor expansions around it at each neighbouring point $j\in S_{i}$

[TABLE]

The unknown coefficients $[\nabla u,D]^{T}$ are computed by a weighted least squares method, by minimizing

[TABLE]

where $W$ is a weighting function used to make sure that the points closer to the central point $i$ have a larger impact than the points farther away. It is usually taken as a Gaussian distribution

[TABLE]

where $\alpha$ is a positive constant. Note that the weighting function is only defined on the local support $S_{i}$ consisting of $n(i)$ points. For the sake of brevity, we present only the case of one spatial dimension. Eq. (2) and Eq. (3) lead to the following system which is solved at each point $i=1,\dots,N$

[TABLE]

where $\delta x_{j}=x_{j}-x_{i}$. Or, in short form $\vec{E}=M\vec{a}-\vec{b}$, with $\vec{a}=(u_{x},u_{xx})^{T}$ and $\vec{b}=(u_{1}-u_{i},\dots,u_{n}-u_{i})^{T}$. The minimization Eq. (3) can be rewritten as

[TABLE]

where $W$ is a diagonal matrix with entries $W_{1},\dots,W_{n}$. A formal minimization leads to

[TABLE]

Which leads to the differential operator stencils

[TABLE]

where $c_{ij}^{x}$ and $c_{ij}^{xx}$ represent the values in the first and second row respectively of the matrix $[(M^{T}WM)^{-1}M^{T}W]$ in Eq. (7). These stencils are then used to solve the PDE. For example, if we consider the PDE

[TABLE]

The derivative approximations Eq. (8), Eq. (9) are substituted into the PDE to obtain

[TABLE]

which forms a large sparse implicit system which is solved with an iterative method.

An alternate, but equivalent, way to arrive at the stencil coefficients in Eq. (8) is to ensure that the derivatives of monomials $m\in\mathcal{M}$, up to the order of accuracy desired, are exactly reproduced [34]

[TABLE]

To improve stability and to improve accuracy for Dirichlet boundary conditions, a Kronecker-delta property is often fulfilled, either by the addition of a Kronecker-delta function $\delta_{ij}$ to the monomial set $\mathcal{M}$ [22] or through the weighting function $W$ [28].

2.2 Tiwari Approach

This method is based on the work of Tiwari and Kuhnert [38]. In this method, stencil coefficients are found that include the PDE being solved as a local constraint. Then a function is found that satisfies the conditions on all the stencils computed. The least squares procedure for minimizing the errors includes not just the Taylor expansions, but also the error in solving the PDE itself. Consider the PDE Eq. (10) used in the previous section. We proceed in the same way as done in Section 2.1. The system Eq. (5) and Eq. (3) get extended to

[TABLE]

which is solved at each point $i=1,\dots,N$. The stencils to the function approximation of $u$ are the only ones of importance. Proceeding in the same way as done in Section 2.1 leads to a system similar to that obtained in Eq. (7), the first row of which gives the function approximation stencils

[TABLE]

We then find a function $u$ that satisfies all these stencils by solving a large sparse implicit system with an iterative method.

[TABLE]

This results in a diagonally dominant system. One of the advantages of this direct approach over the classical GFDM one is that it can be used to solve algebraically over-determined, but well-posed, systems. i.e. systems with more PDEs than variables. An important use of this is to solve PDE systems on boundary points along with the imposed boundary conditions, which is not possible by the classical GFDM. Tiwari and Kuhnert [38] used this framework for pressure-Poisson equations in Navier–Stokes solvers. There, an over-determined system was solved on boundary points, consisting of both the boundary conditions and the pressure-Poisson equation itself. The same idea was later extended for imposing free surface boundary conditions [39], for compressible flows [21], and recently for heat transfer boundary conditions by Reséndiz and Saucedo [33], among other applications. A similar idea of adding the PDE on boundaries was also done in context of the meshfree Radial Basis Functions (RBF), by using an additional set of nodes adjacent to the boundary [12], which was shown to increase accuracy over standard RBF. Using Tiwari’s approach, a similar increase in accuracy is obtained in meshfree GFDMs and can be done without the addition of extra nodes.

However, other than solving a PDE along with the boundary condition at boundary points, to the best of our knowledge, this framework has not been used to solve over-determined problems across all points in the computational domain. In this paper, we show how this framework could be used to devise a new algorithm which reduces some problems with existing meshfree GFDM solvers for the incompressible Navier–Stokes equations by solving an over-determined problem across the entire computational domain.

3 Existing meshfree GFDM algorithms for the incompressible Navier–Stokes equations

We consider the incompressible Navier–Stokes equations in Lagrangian form.

[TABLE]

where $\vec{v}$ is the fluid velocity, $p$ is the pressure, $\rho$ is the density, $\eta$ is the dynamic viscosity and $\vec{g}$ includes both gravitational acceleration and body forces. All the algorithms considered below start with an update of point locations by solving Eq. (33) according to

[TABLE]

for each point $i=1,\dots,N$, where the bracketed superscript refers to the time level. Following the movement of points, mass conservation Eq. (34) and momentum conservation Eq. (35) are solved according to one the methods mentioned below.

Two broad classes of algorithms are used to solve the incompressible Navier–Stokes equations. Namely, operator splitting methods, and monolithic methods. In operator splitting methods, also referred to as fractional step methods (FSM), partitioned methods, or pressure-segregation methods, the momentum conservation and mass conservation equations are solved in two separate steps. In monolithic methods, also referred to as coupled velocity-pressure methods, the mass conservation and momentum conservation equations are solved together in one large system. In the remainder of this section, we present existing meshfree methods in both these classes and study their drawbacks. In the next section, we present a new monolithic meshfree method which helps reduce these problems in existing solvers.

3.1 Operator Splitting Methods

Projection methods are a type of operator splitting methods, based on the Helmholtz-Hodge decomposition [5], which have been widely used for approximating incompressible and weakly compressible fluid flows in Finite Volume Methods [1, 7]. These method have also been commonly used in meshfree contexts (see, for example, [39, 42]). They consist of two steps. In the first step, an intermediate velocity is computed by solving the momentum conservation equation. This intermediate velocity is then projected to a divergence-free field with the help of a correction pressure to obtain the final velocity. The intermediate velocity $\vec{v}^{*}$ is obtained by solving

[TABLE]

where $p^{*}$ is a pressure guess. In the original projection method of Chorin [8], the pressure guess is taken to be zero. Non-zero values of this pressure guess have been shown to produce more accurate results [7]. Throughout this paper, we take the pressure guess to be the pressure at the previous time-step $p^{*}=p^{(n)}$. The step of computation of the intermediate velocity is often done explicitly, by replacing the $\Delta\vec{v}^{*}$ term in Eq. (37) with $\Delta\vec{v}^{(n)}$ [39]. However, the implicit way used in Eq. (37) produces more accurate, and often more stable, results.

In the second step, the velocity is corrected by projecting it to a divergence free space by

[TABLE]

The pressure correction $p_{corr}$ is computed by a pressure-Poisson equation obtained by applying the divergence operator to Eq. (38) and setting $\nabla\cdot\vec{v}^{(n+1)}=0$

[TABLE]

Finally the pressure is updated by

[TABLE]

Different variations of such projection methods including higher order time integration have been studied, for example, by Brown et al. [7]. A variety of boundary conditions have been used in these methods. A discussion on boundary conditions for projection methods in the mesh-based context can be found in Denaro [25]. In meshfree projection methods, a common approach is to obtain the boundary condition by projecting the underlying equation solved at interior points on the outward facing unit normal. Such approaches and the required stabilization can be found in Fang and Parriaux [11] and Boroomand et al. [6]. In this paper, we use boundary conditions dependent on the physics of the underlying problem as done in Seibold [34]. Further, for the spatial discretization of Eq. (37) and Eq. (39) we use the classical meshfree GFDM presented in Section 2.1, as done, for example, by Drumm et al. [10]. We now consider a few drawbacks of such projection methods in the meshfree context.

Consistency of numerical differential operators:

While applying the divergence operator to Eq. (38) to obtain Eq. (39), the assumption is made that the Laplace operator is the same as the divergence of the gradient operator. Specifically,

[TABLE]

While this is certainly true for continuous operators, it does not hold for discrete differential operators. This leads to an error in the approximation of the numerical divergence of the new velocity. Thus, mass conservation is violated at the local level. Unlike the truncation error due to the order of spatial approximation, this error does not converge to zero with a decreasing spatial discretization.

In the context of Finite Volume Methods, this problem has been known for several decades. One proposed solution was to replace the classical discrete Laplace operator with wider stencils by taking the convolution of the divergence operator and the gradient operators, i.e. setting $\Delta:=\nabla\cdot\nabla$ at the discrete level [3]. Another, more widely used approach to overcome this problem is to use staggered grids in which the pressure and velocity fields are defined at different locations [15]. The staggered grid approach can not be generalized to meshfree methods since all properties, both scalars or vectors, are prescribed on the same nodes. Using the first approach of setting the numerical Laplace operator to be equal to the numerical divergence of the numerical gradient causes several issues in the meshfree context. Firstly, this would lead to different support sizes for the first and second derivatives. The second derivatives would be defined on a support double the size of that of the first derivative, which more than doubles the number of points in each support domain. Moreover, there is no control over the center stencil values, which makes diagonal dominance hard to achieve. As a result, convergence of the large linear systems can be troublesome. Thus, this problem still persists in meshfree projection methods. Not only for meshfree GFDMs, but also in various other meshfree approximation methods including SPH [9, 42].

Other operator splitting algorithms such as the PISO algorithm [18], SIMPLE algorithms [31] and their derivatives also rely on pressure-Poisson equations. The difference being that the Poisson equations are derived by applying the divergence operator to the momentum conservation equation. Their meshfree equivalents possess the same drawback of $\Delta\neq\nabla\cdot\nabla$ at the discrete level, which leads to the same inaccuracies in the approximation of the mass conservation condition.

Compressible boundary layer:

When using the classical GFDM approach of Section 2.1 the pressure-Poisson equation Eq. (39) is solved at interior points with appropriate boundary conditions on the boundary points. The divergence of the velocity at boundary points depends on the pressure-Poisson equation at boundary points which is not solved at all. This results in the formation of a numerical boundary layer of compressible fluid, with non-zero divergence of velocity, during the simulation of incompressible fluids. This problem has been alleviated by Tiwari and Kuhnert [39] by solving a system where an over-determined system is solved at boundary points, which considers both the relevant boundary conditions and the pressure-Poisson equation, by using the modified framework of Section 2.2. This problem of numerical boundary layers is also avoided by several solvers that solve the pressure-Poisson system before the implicit velocity system [45]. More recently, alternate solutions to this same problem were also considered by Idelsohn and Oñate [16], but specific to free surface boundaries. However, the earlier issue of $\Delta\neq\nabla\cdot\nabla$ at the discrete level is present even at the boundaries, and is a larger source of error than the compressible boundary layer.

Poor accuracy at low Reynolds flow:

Projection methods usually suffer from low accuracy, and often instability, for fluid flows at low Reynolds numbers. This problem is reduced, but not overcome, by the implicit nature of Eq. (37). As a result, projection methods do not provide good approximations for fluids with high viscosity such as molten glass [29].

3.2 Monolithic Methods

Coupled velocity-pressure solvers usually have the advantage of being more stable, especially for larger time-step sizes, but often come at the disadvantage of ill-conditioned systems. To avoid the problem of ill-conditioned systems, Kuhnert [22] developed a penalty formulation based on the classical GFDM approach of Section 2.1. In two spatial dimensions, for $\vec{v}=(u,v)$ and $\vec{g}=(g_{x},g_{y})$, it can be written in matrix form as

[TABLE]

where $\mathbf{C}^{x}$ is the matrix formed by the stencil coefficients for the numerical differential operators for the $x$ derivative, $c_{ij}^{x}$ (see Eq. (7) and Eq. (8)). Similarly for $\mathbf{C}^{y}$ and $\mathbf{C}^{\Delta}$. The vector $\vec{U}^{(n+1)}$ is formed by $u^{(n+1)}$ at all points in the computational domain, $\vec{U}^{(n+1)}=(u^{(n+1)}_{1},\dots,u^{(n+1)}_{N})^{T}$ and similarly for the other upper case vectors. Thus, the first two blocks of rows in Eq. (42) represents the momentum conservation equation at all points. The last block of rows in Eq. (42) is a penalty formulation for the conservation of mass equation

[TABLE]

Rows corresponding to boundary points in Eq. (42) are replaced with the appropriate problem specific boundary conditions. The ideal scenario would be to set $A=0$, to get an exact conservation of mass. However, that leads to an ill-conditioned system. Further, typical iterative solvers do not converge for such a system due to a complete lack of diagonal dominance in the last block of rows in Eq. (42). While specialized solvers for saddle-point problems [4, 26], could possibly be used to obtain solutions to these sparse linear systems, simulations are still usually unstable when $A=0$. Thus, non-zero values of $A$ need to be used. Kuhnert [22] takes this parameter to be in the range of $A\in(0,0.3)$. Lower values of $A$ are preferred for greater accuracy, however, higher values are needed for better conditioning and faster convergence of the resulting linear system. We note that using $A=1$ would result in an implicit and coupled projection method, with an additional velocity correction step, Eq. (38), required. The final pressure is given as done before $p^{(n+1)}=p^{*}+p_{corr}$.

This penalty approach has been shown to be more stable than the meshfree projection method [19, 22], especially at low Reynolds flows. However, setting $A\neq 0$ in Eq. (43) leads to an artificial compressibility that is numerically observed to be similar to that introduced by the inconsistency between the $\Delta$ and the $\nabla\cdot\nabla$ discrete operators in operator splitting methods. Further, this approach has the same issue of a compressible boundary layer as that in classical operator splitting methods explained earlier.

This penalty approach coupled solver and the meshfree projection method have both been widely used and have shown to be robust methods with a wide variety of applications. Under the name of the Finite Pointset Method (FPM), they have also been used as the numerical basis of two commercially used meshfree simulation tools: NOGRID [29] and the meshfree module of VPS-PAMCRASH [40].

Other approaches to coupled solvers include solving the momentum conservation as in the first two blocks of rows in Eq. (42), with a pressure-Poisson equation replacing the third block of rows in Eq. (42). Like most fractional step methods, this solves the mass conservation indirectly, and introduces the same error as mentioned earlier of $\Delta\neq\nabla\cdot\nabla$ at the discrete level. In mesh-based contexts, especially for Finite Element Methods, the equivalent systems of Eq. (42) with $A=0$ are often solved by algebraically decomposing the system which is essentially equivalent to a fractional step method [35]; or with the help of stability conditions [2] which introduce errors similar to that in the method described in this section.

4 A New Monolithic Solver

We wish to solve the momentum and mass conservation equations as

[TABLE]

We emphasize that the desire is to solve the mass conservation directly as in Eq. (45) and not indirectly via a pressure-Poisson equation. Using the classical GFDM approach of Section 2.1, Eq. (44) and Eq. (45) lead to Eq. (42) with $A=0$. As mentioned earlier, this results in an ill-conditioned system that is very hard to solve with typical iterative procedures. This problem of ill-conditioned systems can be avoided by using Tiwari’s approach presented in Section 2.2. For the same, we start by rewriting Eq. (44) and Eq. (45), in two spatial dimensions, to obtain

[TABLE]

where $e_{1}$, $e_{2}$ and $e_{3}$ are the errors in the discretizations of the respective PDEs. A quadratic minimization of $e_{1}$, $e_{2}$ and $e_{3}$, along with the errors from the Taylor expansions in all three variables, $u$, $v$ and $p_{corr}$ can be done, similar to that done in Section 2.2. This results in a system of equations similar to Eq. (32) and provides better conditioned systems than Eq. (42) with $A=0$. However, resulting simulations are mostly unstable. Thus, such an approach of using Tiwari’s framework of GFDM discretization to coupled solvers has not been used successfully in the past.

A possible explanation of the instability is the lack of information about the pressure correction $p_{corr}$. To correct this and to introduce a further coupling condition between the velocity and pressure, we make use of the fact that this modified GFDM framework can be used to solve algebraically over-determined problems. The system of Eq. (46) – Eq. (48) are augmented with a pressure-Poisson equation to form an over-determined system. The Poisson equation is obtained in a manner similar to that done by SIMPLE and PISO methods, by applying the divergence operator to the conservation of momentum equation Eq. (35) to obtain

[TABLE]

Eq. (46) – Eq. (49) are solved at all interior points. Note that the mass conservation condition is solved directly in Eq. (48) and indirectly in Eq. (49). At boundary points, the mass conservation condition Eq. (48) is solved in addition to the relevant boundary conditions, once again leading to more PDEs than variables. The errors in each of these PDEs or boundary conditions are minimized simultaneously with the errors in the Taylor expansions for each velocity component and the pressure.

[TABLE]

where $e_{j}^{u}$, $e_{j}^{v}$ and $e_{j}^{p}$ are the errors in the Taylor expansions around point $i$ of $u$, $v$ and $p_{corr}$ respectively. For all simulations, $W_{j}$ are taken to as in Eq. (4), and $W_{\text{PDE},k}=2$. Thus, the following least squares system is solved at each interior point to obtain the function approximation stencil coefficients : $\vec{E}=M\vec{a}-\vec{b}$ with

[TABLE]

where $p$ is used as shorthand for $p_{corr}$, and $r_{k}$ are the right hand sides of equations Eq. (46) – Eq. (49). Further,

[TABLE]

with $M_{T}$ being the parts coming from the Taylor expansions

[TABLE]

and the part coming from the PDEs is given by

[TABLE]

The first two rows in Eq. (62) represent the momentum conservation in $x$ and $y$ directions respectively, the third row represents the mass conservation equation, and the last row represents the pressure-Poisson equation. Similar systems are obtained for boundary points, with the relevant rows in Eq. (62) replaced with the boundary conditions, as done in [39]. As done earlier, a formal minimization leads to $\vec{a}=[(M^{T}WM)^{-1}M^{T}W]\vec{b}$. Only the function approximation stencils are of interest to us, i.e. the rows of $u$, $v$ and $p_{corr}$. They can be written as

[TABLE]

for $i=1,\dots,N$. The coefficients $\alpha$, $\beta$, $\gamma$ and $\zeta$ represent the values in the relevant row of the matrix $[(M^{T}WM)^{-1}M^{T}W]$. These can be rearranged to obtain the sparse linear system

[TABLE]

This resulting sparse linear system is solved using an iterative solver. Eq. (66) – Eq. (68) produce a diagonally dominant system. In Eq. (66), for example, the magnitude of the diagonal values of $1-\alpha_{ii}^{u}$ are significantly larger than that of the off-diagonal values of $\alpha^{u}_{ij}\;j\neq i$, $\beta^{u}_{ij}$ and $\gamma^{u}_{ij}$. Thus, the resulting linear system converges well with typical iterative procedures.

A short comparison between the meshfree projection method presented in Section 3.1, the penalty approach coupled solver of Section 3.2 and the new coupled solver presented in this section is listed below

All three methods have the same theoretical convergence rate with the spatial discretization as all use Taylor expansions up to the same order of accuracy. All three methods also use similar first order temporal discretizations. Higher order methods have been studied extensively, especially for mesh-based fractional step methods (for example, [43]). Such higher order approximations could be applied to all three methods considered here.

• 2.

All three methods are “approximate” methods as apposed to “exact” methods in the sense that they only solve the mass conservation equation up to the truncation error. However, the new method does not contain the additional sources of error present in the other two methods. The projection method solves the mass conservation indirectly which leads to errors due to a lack of consistency between the first and second order derivatives. The penalty approach coupled solver attempts to solve the mass conservation directly, but introduces an artificial compressibility to improve conditioning of the resulting system. On the other hand, the new coupled solver solves the mass conservation equation directly and without introducing an artificial compressibility, and thus provides a much better approximation to the mass conservation equation than both other methods.

• 3.

While the penalty approach coupled solver results in a compressible boundary layer, this is not present in the projection method of Tiwari and Kuhnert [39] or the new coupled solver presented here, due to the addition of the mass balance equation on boundary points.

• 4.

Numerically it is observed that the new coupled solver has stability comparable to the penalty approach coupled solver, and thus, much better than the projection method.

• 5.

Both coupled solvers solve one large implicit linear system while the projection method solves two smaller implicit linear systems. In both coupled solvers, the linear systems are of the same size, but the system is denser in the new coupled solver. The resulting system of Eq. (66) – Eq. (68) can be written in matrix form as

[TABLE]

where $\boldsymbol{\alpha}$ , $\boldsymbol{\beta}$ and $\boldsymbol{\gamma}$ are matrices formed from the coefficients in Eq. (66) – Eq. (68) and $\vec{R}_{1}$, $\vec{R}_{2}$ and $\vec{R}_{3}$ are vectors formed from the right hand sides of Eq. (66) – Eq. (68) respectively. While the sparsity pattern of each block of rows in Eq. (69) and Eq. (42) are identical, as they are dependent only on the support domains for each point, several zero blocks are present in Eq. (42), while all blocks are non-zero in Eq. (69). This is no longer true for cases of spatially varying viscosity $\eta$, which are very common in fluid flow applications, for which all blocks are non-zero in the systems arising from both the coupled solvers.

• 6.

The larger size of the implicit linear systems means that both coupled solvers have higher memory requirements than the projection method.

• 7.

The fastest simulation times between the three methods vary on a case by case basis, but the overall simulation time is similar for all three methods.

5 Numerical Results

The explicit time-integration for the movement of points according to Eq. (36) results in a CFL-like condition on the time-step size

[TABLE]

Thus, a varying time-step size according to Eq. (70) is used in all performed simulations. To determine the numerical order of accuracy of the three methods, a small time-step is used in the first test case below by using a small value of $C_{\Delta t}$. The remaining test cases use a much larger time-step which results in lower experimental convergence rates. In the comparison of simulation times, the total time of the simulation refers to the total clock time (in seconds) of the simulation, including the initial point cloud setup and the post-processing error calculations. To ensure that these comparisons are realistic, all three methods were implemented by the same programmer and use identical memory management. All simulations were carried out in Fortran and were run serially on an Intel XeonE5-2670 CPU rated at 2.60GHz. All sparse linear systems are solved using the BiCGSTAB iterative solver [41] without the use of any preconditioner. For the penalty approach coupled solver, the penalty coefficient is taken to be $A=0.1$ in all simulations, since lower values result in much larger simulation times due to poor convergence of the linear systems. In all the figures below, the projection method presented in Section 3.1 is referred to by ‘Projection’, the coupled solver with the penalty formulation presented in Section 3.2 is referred to by ‘Coupled: Penalty’ and the new coupled solver which directly solves an algebraically over-determined system, as done in Section 4, is referred to as ‘Coupled: New’.

5.1 Taylor-Green Vortices

As a validation case, we consider the two-dimensional decaying vortices referred to as the Taylor-Green vortices on $[0,2\pi]\times[0,2\pi]$. The analytical solution is given by

[TABLE]

where $\vec{v}_{a}=(u_{a},v_{a})$ is the analytical solution for the velocity, $p_{a}$ is the analytical solution for the pressure, and $Re=\frac{\rho UL}{\eta}$ is the Reynolds number, with the characteristic velocity $U=1$ and the characteristic length $L=2\pi$. Error in the numerical solution $\vec{v}$ is measured by

[TABLE]

where $V_{i}$ is a volume associated with point $i$. For an error $\epsilon$ and consecutive smoothing lengths $h_{1}$ and $h_{2}$, the rate of convergence of the solution with changing smoothing length is measured as

[TABLE]

The initial condition for the velocity is taken in accordance with the exact solution, and is shown in Figure 1. Dirichlet boundary conditions are used on all boundaries. The simulations are done up to an ending time of $t_{end}=1s$, for $\eta=1Pa\,s$ and $\rho=1kg/m^{3}$, which gives $Re=2\pi$. A small time-step is used according to Eq. (70) with $C_{\Delta t}=0.005$. The convergence of the errors with the smoothing length $h$ are shown in Figure 2 and are tabulated in Table 1. All three methods match the analytical solution well. All three methods exhibit a similar convergence rate, which matches the theoretical expectation. Errors in the new coupled solver are smaller than the other two methods, while the errors in the projection method are the largest, but are only slightly larger than those in the penalty approach coupled solver. Total simulation times for each case are also shown in Table 1. All three methods take approximately the same time, with the new coupled solver being the slowest and the projection method being the fastest. Similar results were obtained for larger Reynolds numbers.

5.2 Bifurcated Tube

We consider flow of a fluid through the bifurcated tube shown in Figure 3. The length of the tube is $60m$, the width is $4m$ in the thick region and $1m$ in the thin region.

Simulation parameters are set as $t_{end}=1s$, $\rho=10^{3}kg/m^{3}$ and $\eta=2Pa\,s$. The velocity at the inflow, on the left of the tube, is kept constant at $\vec{v}_{in}=(2m/s,0)$. This results in a Reynolds number of the order of $10^{3}$. A varying time-step is used according to Eq. (70) with $C_{\Delta t}=0.05$. Homogeneous Neumann boundary conditions for the velocity are used at the outflow and no-slip conditions on the walls. The pressure is kept constant at atmospheric pressure at the outflow and homogeneous Neumann boundary conditions are considered elsewhere for the pressure. The error in mass conservation is measured as the difference between the total volume of fluid flowing in and that flowing out, throughout the entire simulation. The mass conservation error is given by

[TABLE]

where $\vec{n}$ is the outward pointing unit normal and $\partial\Omega_{in}$ and $\partial\Omega_{out}$ are the inflow and outflow boundaries respectively. Note that $\epsilon_{mass}$ measures the errors during transient states too, and not just the errors in the steady state solution. A measure for the velocity divergence throughout the domain is taken as the integral of the divergence of velocity scaled by the total volume

[TABLE]

This can be interpreted as the average value of local error in the mass conservation equation. Note that there is no direct correlation between the measures of divergence and mass conservation, $D(\vec{v})$ and $\epsilon_{mass}$, because the absolute value of velocity divergence is taken in $D(\vec{v})$. The presence of a numerical source and a sink of equal magnitudes for $\nabla\cdot\vec{v}$ would cancel out while measuring the mass conservation, but they would add up while measuring $D(\vec{v})$. Figure 4 shows the convergence of the velocity divergence averaged over all time-steps, and the convergence of the error in mass conservation with respect to changing smoothing length $h$. The same are also tabulated in Table 2. The errors follow a similar pattern to those in the Taylor-Green vortices test case. The new coupled solver produces significantly smaller errors than the other two methods in both the average velocity divergence and mass conservation. The difference between the projection method and the coupled penalty approach is quite small. We note that convergence orders are observed to be much smaller than in the earlier test case due to the use of larger time-steps, Neumann boundary conditions and a non-standard domain geometry.

Simulation times for the three methods are also shown in Table 3. While the new coupled solver is the slowest by a large margin on the coarsest point cloud, it is the fastest on the finest point cloud. This could be due to slower convergence of linear systems in the other two methods.

To illustrate that the difference between the new coupled solver presented in this paper and the penalty formulation coupled solver goes beyond the addition of the zero divergence equation at the boundary points, we split the average divergence of velocity, Eq. (78), along the interior and boundary points $D(\vec{v})=D_{int}(\vec{v})+D_{bnd}(\vec{v})$, with

[TABLE]

where $\Omega\backslash\partial\Omega$ represents only the interior points. Figure 5 shows $D_{int}(\vec{v})$ and $D_{bnd}(\vec{v})$ averaged over all time-steps. It illustrates that the new coupled solver improves the accuracy of the mass conservation condition across both interior and boundary points.

5.3 Sloshing

The most common area of application of Lagrangian meshfree methods is for flows with moving free surfaces. We consider the sloshing of a fluid contained in a constantly moving rectangular box as shown in Figure 6. The dimensions of the initial state of the fluid are $1.2m\times 0.12m$, and that of the box containing the fluid are $1.2m\times 0.6m$.

The initial state is taken to be at rest. Slip boundary conditions are used at the walls for the velocity. The free surface boundary conditions are given by

[TABLE]

where $\vec{t}$ is a unit vector in the tangential direction to the free surface, $p_{0}$ is the atmospheric pressure, and $\mathbf{S}=\mathbf{S}(\vec{v})$ is the stress tensor. Homogeneous Neumann boundary conditions are used for the pressure at the walls. The movement of the box is represented in the gravitational and body forces term by setting $\vec{g}=(2\cos(10t),-10)$. The simulation parameters are set as $t_{end}=3s$, $\rho=10^{3}kg/m^{3}$, and $\eta=0.1Pa\,s$, which leads to a Reynolds number of the order of $10^{4}$. A varying time-step is used according to Eq. (70) with $C_{\Delta t}=0.3$. The error in mass conservation is measured by the change in total volume occupied by all points, since the density $\rho$ is fixed and constant throughout the domain

[TABLE]

where $\Omega_{0}$ is the initial domain and $\Omega_{end}$ is the domain at $t_{end}$. The average velocity divergence is measured as done earlier in Eq. (78). The convergence of volume conservation error and average velocity divergence $D(\vec{v})$ with respect to a changing smoothing length $h$ for all three methods are shown in Figure 7 and are tabulated in Table 3. The results follow a similar trend to that in the earlier two test cases. The new coupled solver shows the highest accuracy in both velocity divergence and mass conservation while the projection method shows the least accuracy in both. However, the difference between the new method and the two older methods is not as large as in the previous test case. The sharp increase in accuracy between $h=0.06$ and $h=0.03$ is due to a bad approximation of the free surface at the coarsest grid. For the finer point clouds, a small convergence rate is observed once again due to the boundary conditions used and the use of large time-steps.

Simulation times for the three methods are also present in Table 3. All three methods take almost the same simulation time, with the new coupled solver being the slowest for 3 out of the 4 point clouds considered.

We note that more turbulent sloshing problems than the one considered here produce larger errors in volume conservation, but for such problems, the largest source of error is in the management of the point cloud and not in the approximation of the PDEs.

5.4 Accuracy of Gradient Reconstruction

In the classical meshfree GFDM approach of Section 2.1, the errors in Taylor expansions are minimized directly. As a result, first order monomials are differentiated exactly, which is evident from the formulation in Eq. (12). However, the same need not be true for the modified approach of Section 2.2 used for the new coupled solver presented in this paper. The addition of the PDE error terms to the functional minimization in Eq. (30) and Eq. (50) could result in larger errors in the Taylor expansions. This would introduce errors in gradient reconstruction.

For $\vec{v}=(u,v)$, we examine this difference numerically by looking at the errors in the Taylor expansions for $u$. Formally, the functional being compared is $J=\sum_{j\in S_{i}}W_{j}(e_{j}^{u})^{2}$ for each point $i$. After the velocity at the new time-level is computed, the obtained velocity is checked for errors in the Taylor expansions. Thus, these errors include not just the truncation error in the discretization, but also the error due to the tolerance of the sparse iterative solver. This comparison is done using the Taylor-Green vortices test case considered earlier in Section 5.1, and for velocities after the completion of the first time step of the simulation. These errors are tabulated in Table 4 for the classical GFDM (with the projection method) and the modified framework used for the new coupled solver. The errors are larger in the new coupled solver, but not significantly.

6 Conclusion

We presented a new monolithic algorithm for the incompressible Navier–Stokes equations, solved using a meshfree Generalized Finite Difference Method (GFDM). While existing algorithms either solve the mass conservation indirectly via a pressure-Poisson equation or introduce an artificial compressibility, in the new method presented here, the mass conservation is solved directly without the introduction of any artificial compressibility. This results in improved local approximations for the mass conservation equation for both interior and boundary points, which results in better accuracy overall. In this method, the momentum and mass conservation equations are solved together, and simultaneously with a pressure-Poisson equation. The addition of this Poisson equation is essential for stabilizing the scheme and results in an over-determined system of PDEs. Accuracy is further improved at boundary points by solving the mass conservation equation in addition to the usual boundary conditions.

This new scheme was compared with two existing methods: one fractional step method and one monolithic method. All three methods have the same spatial and temporal order of accuracy. Numerical comparisons were done for different Reynolds flows, and the new method was shown to produce more accurate results.

Our future work in this direction will be devoted to extending this algorithm to three spatial dimensions. A further interesting point of study is the impact of the ratio of weights in the minimization of the functional in Eq. (50). Higher weights could be used, for example, for the minimization of the error in the mass conservation equation.