Addressing the checkerboard problem in an Eulerian meshless method for incompressible flows

The incompressible Navier-Stokes equations are considered.

where $u$,$v$ denote $x$ and $y$ components of velocity and $p$ denotes pressure. $\rho$ and $\nu$ denote density and kinematic viscosity.
The projection method [33] is used to solve the above equations. The momentum equations are first marched to solve for a provisional velocity field.

The convective and viscous terms are denoted as $C_{u}$, $C_{v}$ and $V_{u}$, $V_{v}$ respectively. The superscript ’n’ denotes the time level and ’*’ denotes the provisional values. The provisional velocity field does not satisfy the divergence-free condition. A correction pressure $p_{c}$ is computed from the Poisson equation.

Here, $\vec{V}^{*}=[u^{*},v^{*}]$. The velocity and pressure are corrected using the correction pressure as

The generalized finite difference method is a meshless method that estimates derivatives of flow variables from a set of neighbours that are a part of a point cloud using a weighted least squares error minimization procedure as discussed below [34, 35].
Consider a point $i$ which has a neighbourhood of points $j\in S_{i}$. $S_{i}$ denotes the support region around the point $i$. This is illustrated in Fig. 1.

For a given operator, there are consistency conditions that are required to be satisfied when applied on monomials up to a prescribed degree, which is usually 2. These monomials are

Here, $\Delta x=x-x_{i}$ and $\Delta y=y-y_{i}$. Now, we derive the procedure for the Laplacian operator. The process is identical for any other operator. Let us denote the Laplacian operator at a point $j$ in the neighbourhood of point $i$ as $C^{\Delta}_{ij}$. Applying Laplacian the operator to each of the monomials in Eq. 10, we get

This can be rewritten in matrix form as

where, $\text{V}=\begin{pmatrix}.&.&.&1&.&.&.\\ .&.&.&\Delta x_{j}&.&.&.\\ .&.&.&\Delta y_{j}&.&.&.\\ .&.&.&\Delta x_{j}^{2}&.&.&.\\ .&.&.&\Delta y_{j}^{2}&.&.&.\\ .&.&.&\Delta x_{j}\Delta y_{j}&.&.&.\\ \end{pmatrix}$ $\vec{C^{\Delta}}=\begin{pmatrix}C^{\Delta}_{i1}\\ .\\ .\\ C^{\Delta}_{ij}\\ .\\ .\\ C^{\Delta}_{iN}\\ \end{pmatrix}$ $\vec{b}=\begin{pmatrix}0\\ 0\\ 0\\ 2\\ 2\\ 0\\ \end{pmatrix}$. Eq. 18 is used in minimizing the functional

This gives

where W is a diagonal matrix with its entries given by the weight of the point $j\in S_{i}$ w.r.t. point $i$, $\textbf{W}=\begin{pmatrix}w_{i1}&0&0&.&.\\ 0&w_{i2}&0&.&.\\ .&.&.&.&.\\ &&.&.&w_{iN}\\ \end{pmatrix}$
Having derived the Laplacian operator from Eq. 20, the Laplacian of a general function $\phi$ at the point $i$ may be written as

For the first derivatives in $x$ and $y$, the operator $C^{\Delta}_{ij}$ in Eq. 21 is simply replaced with the corresponding operator - $C^{x}_{ij}$ and $C^{y}_{ij}$ respectively.
We, now, look at the discretization of the convective terms ($C_{u}$ and $C_{v}$) from the Navier-Stokes in Eqns. 4 and 5. From Eq. 2, the convective terms are

where $f=u^{2}$ and $g=uv$. The derivatives are computed as

The same is done for $C_{v}$. The viscous terms $V_{u}$ and $V_{v}$ involve the Laplacian of the velocity components.

For the pressure Poisson equation, the Laplacian operator is applied to the pressure correction field in the same manner as above. The zero-Neumann conditions are implemented at the boundaries

where $\hat{n}$ is the unit normal at the boundary, whith $n_{x}$ and $n_{y}$ being the individual components.

The test case considered is the lid driven cavity problem. A unit domain is considered in which the top lid slides with unit velocity ($u=1$), as shown in Fig. 2. The boundary conditions are no-slip at the left, bottom and right walls i.e. $u=0$ and $v=0$. The pressure follows the zero-Neumann boundary condition at all boundaries.

The pressure at the bottom, left corner is arbitrarily fixed at zero.

It is noted that the uniform point cloud for the simulation is derived from the nodes of a regular rectangular mesh, without considering the connectivity as given by the indices of the mesh.

In this section, we introduce the issue of pressure checkerboard pattern for the well-known lid driven cavity problem on a uniform point cloud of 1600 points. The checkerboard problem does not appear for low Reynolds numbers ($R_{e}$). For $R_{e}\geq 400$, pressure oscillations are noticeable, although the solution may still converge to the physically correct solution. For $R_{e}\geq 1000$, the issue of pressure checkerboard pattern becomes prominent. For $R_{e}=3200$, the numerics converge to a non-physical solution with checkerboard oscillations and does not allow the physically correct solution to evolve. The checkerboard pattern in pressure and velocity ($u$) is shown in Fig. 3. It can be seen at the top right corner of the box that adjacent points tend to take on oscillating values of pressure and velocity and these oscillations are retained as a permissible solution. This issue becomes more severe as $R_{e}$ is increased further, as seen in Fig. 4 in which $R_{e}=5000$.

In order to resolve the problem, multiple prospective solutions were attempted as a part of this work. However, they do not alleviate the problem. For the sake of completeness of this work, these attempts are documented below, along with a brief explanation.

The source of the checkerboard problem in mesh-based finite difference and finite volume methods is traced to the lack of participation of the point in consideration, in the evaluation of the pressure gradient at the same point. In other words, a hypothetical pressure field that oscillates between two distinct values ($p_{a}$ and $p_{b}$) from one point to the next, is still seen as a uniform pressure field since the pressure at the neighbouring points is the same.This is illustrated in Fig. 7.

In meshless methods, however, the neighbourhood of a point involves points that are not just in the immediate proximity, but also points that are farther away. Therefore, the hope is that the participation of these points in the estimation of the pressure gradient would make the checkerboard oscillations visible to the gradient operator and therefore, the oscillations would be suppressed. A regularly spaced neighbourhood of a point is shown in Fig. 8(a). The red and the blue colours indicate two distinct values between which pressure oscillates in the neighbourhood. The gradient operator is constructed at the point $P$. The numerical value of the differential operator, $C^{x}_{Pj}$,associated with each neighbour $j$, for the first derivative in $x$, is mentioned at corresponding locations in Fig. 8(b). It is evident from the values of the operator that neighbours that are positioned symmetrically w.r.t point $P$ have operators of the same magnitude and opposite sign. Evaluating the pressure gradient using these operators,

This results in the two summations in the above equation to be zero. Therefore, using these operators in the computation of the derivative of pressure at $P$, would result in a zero pressure gradient, just as seen in mesh-based methods. In summary, although the oscillations in pressure are present in the neighbourhood of the point $P$, the operators do not detect the presence of the oscillations.

The two most well-known solutions to the problem of checkerboard oscillations in mesh-based methods are – (a) Staggered meshes and (b) Rhie-Chow interpolation. In staggered meshes, the pressure and the velocities are stored at different locations such that grid the formed by the pressure nodes is staggered w.r.t to the grid formed by the velocity nodes. The momentum equations are solved at the velocity nodes and the pressure gradient at these nodes is estimated from the pressure nodes that are staggered w.r.t the velocity node. Due to this, oscillations in pressure are detected by the system of equations, thus, solving the problem of checkerboarding. In the case of Rhie-Chow interpolation, collocated grids are used, where the pressure and velocity are stored at the same locations. The face velocities are estimated as the sum of the average of the velocities of the nodes that straddle the face and a correction term that depends on the third derivative of pressure at the face. This correction acts as a damping mechanism for any high frequency oscillations that may develop in the numerical solution[3].

In the context of meshless methods, these solutions are not straightforward to apply. Both the above solutions are designed in the context of finite-volume methods where a particular variable is solved for its cell-centred value and the fluxes at the cell faces are summed for each equation. Staggered meshes could correspond to staggered point clouds i.e. two sets of point clouds, one for pressure and another for velocity. For each point, two sets of differential operators will have to be computed. Evidently, the difficulty with this approach is that it dramatically increases memory requirements, number of operations and the amount of book-keeping required in a simulation. It also involves the generation of a staggered point cloud w.r.t a base point cloud, which could be non-trivial. Rhie-Chow interpolation, on the other hand, would require a solution to the integral form of the Navier-Stokes equations, which is the operating framework of the finite-volume method. In a meshless method like GFDM, the equations are solved in the differential form. Although, flux-based formulations can be used in GFDM as well, it relies on error minimization in flux derivatives rather than summing of fluxes over a control volume.

An error indicator based on the function value (in this case, pressure) is proposed. Using the procedure outlined in Sec. 2.2, an operator for the function approximation can be developed by setting

We denote the operator at a point $i$ with neighbouring points $j\in S_{i}$ as $C^{0}_{ij}$ such that the pressure from the weighted least squares approximation is given by

The error ($e$) is defined as

In this equation, $p_{n}$ is the pressure obtained from the numerical solution that contains the checkerboard oscillations. On the other hand, $p_{i}$ comes from a polynomial approximation of GFDM as shown in Eq. 33 and hence, will be a smooth function. As a result, the error $e_{i}$ will be high in regions where the checkerboard oscillations are prominent and low in other regions. This is pictorially illustrated in Fig. 9(a).

Fig. 9(b) shows the error distribution in the lid driven cavity problem for $R_{e}=3200$; the corresponding pressure field is shown in Fig 3(a). It is seen that the error is high in regions of the checkerboard oscillations.
As a criteria for point cloud adaptation, we define

A threshold value of $\psi\geq 0.1$ is set in this work for point cloud adaptation.

The error indicator, as described in the previous section, is used as the criteria for adding points in the neighbourhood of a given point where the error is higher than the threshold. Consider the central point $i$ in Fig. 10.

The set of blue points are a part of the original point cloud. A new point is added at the mid-point between the point $i$ and each of its neighbours $j\in S_{i}$. Before the addition of a new point to the point cloud, the algorithm checks if there already exists a point at that location either from the original neighbourhood $S_{i}$ or an earlier adaptation at one of its neighbours $j$. A minimum distance criteria $d_{\text{min}}$ is set for the insertion of the new point such that the new point lies at least at $d_{\text{min}}$ from any of the existing points in the point cloud. This is done to ensure that the operators of GFDM, described in Sec. 2.2, are not unstable or singular.
The new point is initialized with the mean of the solution variables at the two points between which it is being inserted.

After the insertion of new points, the recalculation of the neighbourhood of points and differential operators is necessary. The adaptation of the point cloud would, in general, affect both the neighbourhood of points that have and have not participated in the adaptation procedure. Any change in neighbourhood would also result in the change in the differential operators. The advantage of using an Eulerian method is that estimation of neighbourhoods and differential operators need to be done only at the beginning of the simulation and whenever point cloud adaptation occurs. In Lagrangian methods, the points are continuously moving and would therefore, require these recalculations every time-step. The point cloud adaptation algorithm is summarized in Algorithm 1.

In this test case, the Reynolds number is set to $3200$. A test case with 1600 points (extracted from a regular $40\times 40$ grid) which exhibits checkerboard patterns, is used as reference. After the simulation is run for sometime on this base point cloud, the adaptation procedure is used to add points in regions of checkerboarding. The adapted point cloud is shown in Fig. 11.

The adapted point cloud contains 2585 points. For the sake of comparison of the extent of checkerboarding, we also use a finer uniform point cloud of 3600 points without adaptation(extracted from a regular $60\times 60$ grid). It is noted that the numerical results such as in Figs. 12 and 13 are presented here with the ’Delauney2D’ feature of Paraview that triangulates the point cloud before rendering the results. This method was preferred since it provides clearer visualization of the flow-field than using just points.
Figs. 12 and 13 show $p$ and $u$ contours on the three point clouds. While in the subfigures (a) and (b), the solution is on a uniform point cloud of 1600 and 3600 points, the subfigure (c) is on an adapted point cloud of 2585 points. It is noted that the checkerboard patterns in $p$ and $u$ are prominent in both the uniform point clouds, with the finer point cloud showing less oscillations. In comparison to both the uniform point clouds, the solution on the adapted point cloud is relatively smooth and manages to counter the checkerboard oscillations satisfactorily. It is also important to note that the adapted point cloud although has only 2585 points and it still performs better than a fine uniform point cloud of 3600 points.
In the cases with the checkerboard patterns, the solution stops changing after a certain point and converges to a non-physical solution. This can be observed in the streamline patterns of the converged solution in the three cases, as shown in Fig. 14. In Fig. 14(a) and (b), the solution on the uniform point cloud does not converge to form the secondary vortices despite sufficient number of points being present to resolve them as flow features. In Fig. 14(c), however, the adaptation of the point cloud suppresses the formation of spurious checkerboard oscillations and allows the solution to converge to the physically-correct solution. Fig. 15(a) reveals a good match in the $u$ velocity profile at the vertical mid-section with that of Ghia et al.[36].

In this test case, the Reynolds number is set to $5000$. A test case with 3600 points (extracted from a regular $60\times 60$ grid) which exhibits checkerboard patterns, is used as reference. The adapted point cloud is shown in Fig. 16.

The adapted point cloud contains 4900 points. For the sake of comparison of the extent of checkerboarding, we also use a finer uniform point cloud of 4900 points without adaptation.
Figs. 17 and 18 show $p$ and $u$ contours on the three point clouds. While in the subfigures (a) and (b), the solution is on a uniform point cloud of 3600 and 4900 points, the subfigure (c) is on an adapted point cloud of 4900 points. It is noted that the checkerboard patterns in $p$ and $u$ are prominent in both the uniform point clouds. In comparison to the case of $R_{e}=3200$, the oscillations have a more severe effect on the convergence. This is evident in the $u$ contours in Fig. 18. The solutions on the uniform point clouds are far from the correct solution, which is obtained on the adapted point cloud (Fig. 18(c)).
In Fig. 19, the streamlines on the adapted point cloud show the formation of the secondary vortices. These features do not form on the uniform point clouds. Fig. 15(b) reveals a good match in the $u$ velocity profile at the vertical mid-section with that of Ghia et al.[36].