An explicit meshless point collocation solver for incompressible Navier-Stokes equations

We present a strong form, meshless point collocation explicit solver for the numerical solution of the transient, incompressible, viscous Navier-Stokes (N-S) equations in two dimensions. We numerically solve the governing flow equations in their stream function-vorticity formulation. We use a uniform Cartesian embedded grid to represent the flow domain. We compute the spatial derivatives using the Meshless Point Collocation (MPC) method. We verify the accuracy of the numerical scheme for commonly-used benchmark problems including lid-driven cavity flow, flow over a backward-facing step and vortex shedding behind a cylinder. We have examined the applicability of the proposed scheme by considering flow cases with complex geometries, such as flow in a duct with cylindrical obstacles, flow in a bifurcated geometry, and flow past complex-shaped obstacles. Our method offers high accuracy and excellent computational efficiency as demonstrated by the verification examples, while maintaining a stable time step comparable to that used in unconditionally stable implicit methods.


Introduction
In this paper, we describe the development of a strong form meshless point collocation method for numerically solving the two-dimensional, non-stationary, incompressible Navier-Stokes (N-S) equations in their stream function-vorticity formulation. Our method relies on the ability of meshless methods to accurately compute spatial derivatives, and on their flexibility to solve partial differential equations (PDEs) in complex geometries. The proposed numerical method can efficiently handle uniform Cartesian point clouds embedded in a complex geometry, and/or irregular sets of nodes that directly discretize the flow domain. To compute the spatial derivatives that appear in the governing flow equations, we apply an interpolation meshless scheme, which computes spatial derivatives on both Cartesian embedded grids and irregular nodal distributions.
The mathematical formulation of the N-S equations, depending on the choice of dependent variables, can be classified as (i) primitive variables, (ii) stream function-vorticity, and (iii) velocity-vorticity. The majority of numerical methods are expressed in the primitive variables The governing equations for the non-stationary, viscous, laminar flow of an incompressible fluid are based on conservation of mass and momentum. The non-dimensional form of the stream function-vorticity (ψ-ω) formulation is: where Re = UL v is the Reynolds number, with U and L being a characteristic velocity and length, respectively, and v the kinematic viscosity of the fluid. Taking the spatial derivative of the stream function ψ yields the velocity components u and v, as follows: We require a solution to Equations (1) and (2), defined in the spatial domain Ω with boundary ∂Ω, given the initial conditions at time t = 0 u = u 0 (5) and boundary conditions u = u ∂Ω (7) ω = (∇ × u ∂Ω )·k (8) withk being the unit vector in the direction of the z axis of a three-dimensional Cartesian coordinate system. We solve the governing equations using an explicit Euler time integration scheme for the transient terms and a strong form meshless collocation method to discretize the flow equations.

Discretization Corrected Particle Strength Exchange (DC PSE) Method
The accuracy of the proposed scheme relies on the accurate computation of spatial derivatives for the unknown field functions in Equations (11)- (14). In the present study, we use the Discretization Corrected Particle Strength Exchange (DC PSE) method, which, to our knowledge, is one of the most accurate meshless methods for computing spatial derivatives. The DC PSE method originated as a Lagrangian particle-based numerical method [34] and is based on Particle Strength Exchange (PSE) operators. To solve PDEs using the DC PSE meshless method, the authors in [35] reformulated the Lagrangian DC PSE method to work in the Eulerian framework. For completion, the PSE operators and the DC PSE method are described below.

Particle Strength Exchange (PSE) Operators
The Particle Strength Exchange (PSE) method utilizes kernels to approximate differential operators that conserve particle strength in particle-particle interactions. The PSE method was proposed by Degond and Mas-Gallic [36] for diffusion and convection-diffusion problems. Eldredge et al. [37] then developed a framework for approximating arbitrary derivatives.
In general, a PSE operator Q β f (x) for approximating the derivative D β f (x) has the form with η β ε = η β (x/ε)/ε d a scaled kernel of radius ε, and d the number of dimensions. The sign in Equation (15) is negative when |β| is even and positive when |β| is odd, with β a multi-index [34]. Now, let β = (β 1 , β 2 , . . . , β n ), where β i , i = 1, 2, . . . , n is a non-negative integer. Then the partial differential operator D β can be expressed as D β = . The challenge is to find a kernel η β ε that leads to good approximations for D β . To find such kernels for arbitrary derivatives we adopt the idea in Eldredge et al. [37] and start from the Taylor expansion of a function f (y) about a point x: Fluids 2019, 4, 164 5 of 32 Next, we subtract or add f (x), depending on whether |β| is odd or even, to both sides and convolute the equation with the unknown kernel η β ε . Considering Equation (15), this leads to: Introducing the continuous α-moments M a = (x − y) α η β (x − y)dy = z α η β (z)dz (18) and isolating the derivatives D β on the right-hand side of Equation (17) we obtain: Finally, to approximate Q β f (x) with order of accuracy r, the following set of conditions is imposed for the moments M a : In addition, if we impose the mollification error ε (x) = D β f (x) is bounded [34]. The procedure to construct a kernel that satisfies the conditions in Equation (16) has been described in [34]. Once the kernel is defined, the operator in Equation (17) can be discretized using a midpoint quadrature over the nodes as where N(x) is the number of all nodes in a neighborhood around x, which can be defined by a cut-off radius r c , usually chosen such that N(x) coincides to a certain level of accuracy with the kernel support, and V p is the volume associated with each particle. Given such a discretization, the discretization error is also bounded [34].

The Discretization Corrected PSE Operators
The Discretization Corrected PSE operators were introduced by Schrader et al. [34] to reduce the discretization error h (x) in the PSE operator approximation. We can derive the DC PSE approximation from Equation (22), aiming of finding a kernel function which minimizes the difference between this discrete operator and the actual derivative. We can achieve this, by replacing each term f (x p ) in Equation (22) with its Taylor expansion about x. This leads to the following expression for the derivative approximation:  (24) and the discrete moments defined as: Therefore, the set of moment conditions becomes: for all β 0, where a min is one when β is odd and zero when β is even. The kernel η β is chosen as: The kernel function consists of a polynomial correction function P(x, z) and the weight function W(z). Different types of weight functions can be applied with the possible choices described in [38].
The unknown coefficients a γ (x) are obtained by requiring the kernel given by Equation (28) to satisfy the conditions in Equation (27), resulting in the following linear system of equations: with weights To compute the approximated derivative Q β h f (x) at node x p , the coefficients are found by solving the linear system of Equation (28) for x = x p . Given our choice of kernel function, the DC PSE derivative approximation becomes: where p(x) = [p 1 (x), p 2 (x), p l (x)] and a(x) are the vectors of terms in the monomial basis and their coefficients, respectively. By using the DC PSE method, the spatial derivatives Q β up to second order are given as: and

Vorticity Boundary Conditions
In finite difference methods, a number of formulae can be used to impose the vorticity boundary conditions [33,39]. These formulae, generally referred to as local, compute vorticity values on a given node on the boundary from the vorticity values in the interior domain, without involving other nodes on the boundary. The most widely used method is Thom's formula [39]. Unfortunately, applying these formulae in real applications had limited success. As stated in [33], vorticity boundary conditions should be global, in the sense that computing the boundary vorticity values should involve vorticity values at nodes in the interior and the boundary. Herein, to achieve global vorticity boundary conditions, we impose vorticity boundary conditions through strong form meshless differential operators, used to compute spatial derivatives. This way, imposing vorticity boundary conditions is consistent with the rest of the algorithm (described in Section 2.2) and becomes quite straightforward.
DC PSE methods are recognized as one of the most accurate and efficient numerical methods to compute derivatives on an irregularly distributed set (cloud) of nodes [34,35,38]. For the stream function-vorticity N-S Equations, given the velocity field values computed previously by the updated stream function values ψ n+1 (Step 2, Section 2.2), we can compute the updated vorticity values ω n+1 for the entire spatial domain (including boundaries) using the strong form meshless operators for first order spatial derivative as We have already computed the updated values ω n+1 at the interior nodes (Step 1 in the Euler explicit solver). The updated vorticity values on the boundary nodes are computed using the updated velocity values and the DC PSE operators (described in Section 2.3) for spatial derivatives (Equation (27)).

Poisson Solver
To compute the stream function field values (Equation (2)), we numerically solve a Poisson-type equation. Depending on the flow problem and spatial discretization, a direct or iterative solver is chosen accordingly. Iterative and direct solvers for elliptic type PDEs are well established. Direct solvers are robust and easy to use but can be computationally costly and memory intensive when the systems of equations to be solved are very large (e.g., a few million degrees of freedom). On the other hand, a decisive advantage of iterative solvers is their low memory usage, which is significantly less than a direct solver for the same sized problems. Unfortunately, iterative solvers are less robust than direct solvers and may fail to compute the numerical solution.
It is computationally more efficient to use iterative solvers to solve the Poisson equation when a large number of nodes is used. There are many iterative solvers for Poisson type equations. For example, Gauss-Seidel or successive over relaxation (SOR) algorithms have been widely used. These solvers are accurate but relatively computationally demanding because they need O (N 2 ) iterations to converge (N being the total number of grid points). Multi-grid algorithms have been developed to accelerate convergence and the computational effort has been significantly reduced [40]. For equally-spaced grids, fast Fourier transform (FFT) solvers are the fastest available algorithms for solving Poisson equations [41].
In the present study, we use a direct solver. The discrete Laplace operator defined in the Poisson type equation is accurately discretized using DC PSE method, which performs accurately on uniform and locally refined Cartesian grids. The grid is fixed in time so a LU factorization can be precomputed, which reduces the computational cost significantly. For up to 4 million nodes, the factorization is fast and requires low memory. For larger number of nodes, we utilize iterative solvers, such as conjugate gradient (CG) and biconjugate gradient stabilized (BiCGSTAB), since the linear systems are always positive definite, diagonally dominant and symmetric.

Critical Time Step
To compute the critical time step needed to obtain stable results, we rewrite Equation (1) as which can be written in matrix form as where ω n+1 is a column vector of dimension N × 1 (N is the number of nodes). The N × N matrix A is defined as the DC PSE approximation of the convection and diffusion terms, viz. where and with diag Q 1,0 ψ n and diag Q 1,0 ψ n being diagonal matrices of dimension N × N. The explicit method defined by Equation (34) where λ A are the eigenvalues of the matrix A. The above stability condition requires the eigenvalues λ A to be either real or negative, leading to or complex with negative real parts, leading to When the problem discretization includes a large number of nodes, matrix A becomes very large, and calculating the eigenvalues is not practical. The terms of matrix A are determined by the Laplacian operator L, which consists of second order derivatives, and by the advection operator K, which includes only first order derivatives. The relative weight of the two operators on the structure of matrix A is dictated by discretisation and velocity field values (velocity is related to stream function and vector potential field values); a more refined discretization leads to a higher weight of operator L (diffusion) and lower influence of operator K (convection), while higher Reynolds number and higher velocity field values leads to a higher influence of operator K. When the nodal discretization is not adequately refined, the higher weight of the operator K can lead to eigenvalues with a positive real part; in such case explicit time integration is not possible and discretization has to be refined.
For refined discretizations, the resulting matrix is diagonally dominant, having its eigenvalues distributed close to the real axis. We can use the Gershgorin circle theorem [31] to estimate the maximum time step using Equations (36) and (37) with an upper bound of the maximum eigenvalue of matrix A. Given the composition of matrix A, the upper bound of the maximum eigenvalue can be computed without assembling the matrix: Equation (42) clearly shows that the eigenvalues of matrix A depend on the spatial derivatives of the field variables, computed using the DC PSE method. Therefore, how derivatives are computed reflects on the critical time step. We notice that for the same Reynolds number and the same nodal discretization, the critical time step increases when spatial derivatives are computed using a large number of nodes (more than 40) in the support domain of each node, compared with the critical time step computed using a small number of nodes (around 15). In some cases, this increase in the critical time step can be up to two orders of magnitude, which decreases the computational cost dramatically.
The reason for the increased critical time step is the upwind inherent nature of meshless methods. Meshless methods are locally based methods and use neighboring nodes to compute spatial derivatives. Furthermore, since matrix A involves derivatives of the stream function, by computing the critical time step using the Gershgorin circle theorem we ensure that the computation resembles the estimation of the critical time step based on the nodal spacing ∝ ∆x −1 and the Courant−Friedrichs−Lewy (CFL) condition ∝ uRe∆x −1 . Therefore, by adjusting the number of support domain nodes, we can adjust the critical time step to be comparable to that used by implicit time integration schemes.

Algorithm Verification
To demonstrate the efficiency and accuracy of the proposed numerical scheme, we first apply our method to well-established benchmark problems in computational fluid dynamics, including lid-driven cavity flow, the backward-facing step (BFS), and the unbounded flow past a cylinder.
For all the examples presented here, we examine the accuracy and the efficiency of the proposed scheme for each benchmark problem. Furthermore, we examine the dependence of the critical time step on the grid resolution, Reynolds number and nodes in the support domain. Computations were conducted using an Intel i7 quad core processor with 16 GB RAM.

Lid-Driven Cavity
Our first example involves incompressible, non-stationary flow in a square cavity ( Figure 1). Despite its geometrical simplicity, the lid-driven cavity flow problem exhibits a complex flow regime, mainly due to the vortices formed in the center and at the corners (bottom left and right, and upper left) of the square domain. We impose no-slip boundary conditions at the bottom and vertical walls (left and right wall), while the top wall slides with unit velocity (U = 1), generating the interior viscous flow. For all the examples presented here, we examine the accuracy and the efficiency of the proposed scheme for each benchmark problem. Furthermore, we examine the dependence of the critical time step on the grid resolution, Reynolds number and nodes in the support domain. Computations were conducted using an Intel i7 quad core processor with 16 GB RAM.

Lid-Driven Cavity
Our first example involves incompressible, non-stationary flow in a square cavity ( Figure 1). Despite its geometrical simplicity, the lid-driven cavity flow problem exhibits a complex flow regime, mainly due to the vortices formed in the center and at the corners (bottom left and right, and upper left) of the square domain. We impose no-slip boundary conditions at the bottom and vertical walls (left and right wall), while the top wall slides with unit velocity ( = 1), generating the interior viscous flow. In most numerical studies on the lid-driven cavity flow problem, the flow regime reaches a steady state solution for Reynolds (Re) numbers lower than Re = 10,000. In our study, we report on the results obtained by the proposed scheme for Reynolds numbers Re = 7500 and 10,000. We consider both uniform Cartesian and irregular nodal distributions. For uniform nodal distributions, our previous studies [27][28][29]35] show that for Reynolds numbers up to Re = 10,000, a uniform Cartesian grid with resolution of 361 × 361 captures the forming vortices and provides a grid-independent In most numerical studies on the lid-driven cavity flow problem, the flow regime reaches a steady state solution for Reynolds (Re) numbers lower than Re = 10,000. In our study, we report on the results obtained by the proposed scheme for Reynolds numbers Re = 7500 and 10,000. We consider both uniform Cartesian and irregular nodal distributions. For uniform nodal distributions, our previous studies [27][28][29]35] show that for Reynolds numbers up to Re = 10,000, a uniform Cartesian grid with resolution of 361 × 361 captures the forming vortices and provides a grid-independent and accurate numerical solution. In the present study, to highlight the efficiency of the proposed scheme, we use successively finer uniform Cartesian grids, from 601 × 601 to 2048 × 2408 nodes. We notice that for Re = 10,000 the 1024 × 1024 grid resolution can capture all the flow scales and can be used for Direct Numerical Simulation (DNS) studies.
We use the Gershgorin circle theorem to define the critical time step for different nodal distributions (Tables 1 and 2) and for various Reynolds numbers. Table 1 lists the critical time step computed using the Gershgorin circle theorem for several grid resolutions and different Reynolds numbers (up to Re = 10,000). For the results reported in Table 1, 20 nodes were used in the support domain. We further examine the dependence of the critical time step on the number of neighboring nodes in the support domain. Table 2 lists the critical time step computed using the Gershgorin circle theorem with 40 nodes in the support domain. We can observe that the critical time step increases with the number of nodes in the support domain.

Grid Resolution
Re = 5000 Re = 7500 Re = 10,000 Table 3 lists the computational time (in s) for computing the spatial derivatives for various grid resolutions, and for the numerical solution of N-S equations (for each time iteration) in the case of Re = 10,000. Recall that we precompute the LU factorization for the discrete Laplacian operator defined for the Poisson type equation for the stream function. Hence, steps 1 and 2 must be performed only once at the beginning of the simulation.  Figure 2 shows the uvelocity profiles along the vertical line passing through the geometric center of the cavity, and the vvelocity profiles along the horizontal line passing through the geometric center of the cavity, for Re = 7500 and 10,000, obtained using the proposed scheme and compared to those reported by Ghia et al. [42], which are considered as benchmark solution.  Figure 2 shows the u-velocity profiles along the vertical line passing through the geometric center of the cavity, and the v-velocity profiles along the horizontal line passing through the geometric center of the cavity, for Re = 7500 and 10,000, obtained using the proposed scheme and compared to those reported by Ghia et al. [42], which are considered as benchmark solution.   Table 3 demonstrates the efficiency of the proposes scheme and its ability to solve relatively large algebraic systems (Poisson equation for stream function) with low computational cost.
We set the Reynolds number to Re = 10,000 and the total time for the flow simulation to T final = 250. At that time, the flow regime has all the characteristic features of steady state [42]. We use a 1024 × 1024 grid resolution and a time step for the simulations of dt = 5 × 10 −4 .
Figures 3-5 show stream function and vorticity contours and demonstrate the formation of secondary vortices (Figures 3 and 4) which appear as the Reynolds number increases. There are often used as qualitative results, to highlight the accuracy of the numerical method. In detail, Figure 3 shows the stream function contour plots for Re = 7500, while Figure 4 shows the stream function contour plots for Re = 10,000. Magnified views of the secondary vortices are also included. The stream function values for these contours are listed in Table 4.  (Figures 3 and 4) which appear as the Reynolds number increases. There are often used as qualitative results, to highlight the accuracy of the numerical method. In detail, Figure 3 shows the stream function contour plots for Re = 7500, while Figure 4 shows the stream function contour plots for Re = 10,000. Magnified views of the secondary vortices are also included. The stream function values for these contours are listed in Table 4.     Figure 5 shows the vorticity contours for Re = 7500 and 10,000. The vorticity values for the contours are listed in Table 4.  To verify the applicability of the proposed scheme in irregular nodal configurations, we compute the flow regime for Re = 10,000 using randomly distributed nodes. We use 496,987 nodes to discretize the spatial domain, which provides a grid-independent numerical solution (for the uniform Cartesian grid, a grid independent solution was obtained using a 401 × 401 grid resolution with 160,801 nodes). We obtain the irregular point cloud through a 2D triangular mesh generator (MESH2D-Delaunay-based unstructured mesh-generation). To compare the results with those obtained using the uniform nodal distribution (regular Cartesian grid), we interpolate the velocity, stream function and vorticity values computed using the irregular point cloud onto the uniform nodal distribution using the moving least squares (MLS) approximation method [43] and compute the maximum   Figure 5 shows the vorticity contours for Re = 7500 and 10,000. The vorticity values for the contours are listed in Table 4.

Stream Function Vorticity Contour Level Value of ψ Contour Number Value of ψ Contour Number Value of
To verify the applicability of the proposed scheme in irregular nodal configurations, we compute the flow regime for Re = 10,000 using randomly distributed nodes. We use 496,987 nodes to discretize the spatial domain, which provides a grid-independent numerical solution (for the uniform Cartesian grid, a grid independent solution was obtained using a 401 × 401 grid resolution with 160,801 nodes). We obtain the irregular point cloud through a 2D triangular mesh generator (MESH2D-Delaunay-based unstructured mesh-generation). To compare the results with those obtained using the uniform nodal distribution (regular Cartesian grid), we interpolate the velocity, stream function and vorticity values computed using the irregular point cloud onto the uniform nodal distribution using the moving least squares (MLS) approximation method [43] and compute the maximum normalized root mean square error (NRMSE), defined as [44], For all field variables considered, the NRMSE was less than 10 −4 , which highlights the accuracy of the proposed method for irregular nodal distributions.
We further demonstrate the accuracy of the proposed scheme, providing quantitative results. We use a uniform Cartesian grid of 1024 × 1024 resolution and we set the time step to dt = 5 × 10 −4 . In Table 5, we present the strengths and positions of the primary (located in the center of the cavity), secondary (bottom right corner) and ternary vortexes (bottom and top left corners). The numerical results are in excellent agreement with the results reported in Ghia et al. [42]. Table 6 lists the uvelocity values along vertical line through the center of cavity, and Table 7 the vvelocity values along horizontal line through the center of cavity.

Backward-Facing Step
The second benchmark problem considered is the backward-facing step [45]. This problem has been studied by several researchers using different numerical methods and is considered to be a demanding flow problem to solve, mainly due to the flow separation that occurs when the fluid passes over a sharp corner and re-attaches downstream [1,45,46]. Figure 6 shows the spatial domain for the backward facing step. The coordinate system is centered at the step corner, with the x-coordinate being positive in the downstream direction, and the ycoordinate across the flow channel. The height H of the channel is set to H = 1 (ranging from (0, −0.5) to (0, 0.5)), while the step height and upstream inlet region are set to H/2. To ensure fully developed flow, the downstream channel length L is set to L = 30 H. The inlet velocity has a parabolic profile, with horizontal component u(y) = 12y − 24y 2 for 0 ≤ y ≤ 0.5, which gives a maximum inflow velocity of u max = 1.5 and average velocity of u avg = 1. At the outlet, we assume fully developed flow (du/dx = 0, v = 0). The Reynolds number is defined as Re = u avg H/v f , with v f being the kinematic viscosity.  We first discretize the flow domain with a uniform Cartesian grid. In our previous studies [29,35], we have shown that grid with resolution 31 × 901 is sufficient to compute a gridindependent numerical solution for Reynolds numbers up to Re = 800. As in the previous section, to highlight the efficiency of our scheme, we use a denser grid with resolution 121 × 3630, resulting in 439,230 nodes. We set the total time Ttotal = 250 (to ensure that the solution will reach steady state), and the time step dt = 10 −4 . The simulation terminates when the NRMSE of the time derivative of the stream function and vorticity field values in two successive time steps is less than 10 −6 . The time needed to create the grid was 0.023 s, and it takes 0.3 s to update the solution for each time step. We obtained a steady state solution after 100,000 time steps. Figure 7 shows the stream function and vorticity contours for Re = 800. The flow separates at the step corner and vortices are formed downstream. We can observe the two vortices formed at the lower and upper wall. After reattachment of the upper wall eddy, the flow in the duct slowly recovers towards a fully developed flow. We compute the separation and reattachment points at Llower ≈ 6.1 for the lower wall separation zone, Lupper ≈ 5.11 for the upper separation zone, where the separation begins at x ≈ 5.19. Our numerical findings show good agreement with other numerical methods for 2D computations [1,46]. In [1], the authors used a finite difference method and predicted separation lengths of Llower ≈ 6.0 and Lupper ≈ 5.75, while [46] using the A Fluid Dynamics Analysis Program (FIDAP) code (Fluid dynamics International Inc, Evanston, IL, USA) predicted Llower ≈ 5.8 and the upper Lupper ≈ 4.7.  We first discretize the flow domain with a uniform Cartesian grid. In our previous studies [29,35], we have shown that grid with resolution 31 × 901 is sufficient to compute a grid-independent numerical solution for Reynolds numbers up to Re = 800. As in the previous section, to highlight the efficiency of our scheme, we use a denser grid with resolution 121 × 3630, resulting in 439,230 nodes. We set the total time T total = 250 (to ensure that the solution will reach steady state), and the time step dt = 10 −4 . The simulation terminates when the NRMSE of the time derivative of the stream function and vorticity field values in two successive time steps is less than 10 −6 . The time needed to create the grid was 0.023 s, and it takes 0.3 s to update the solution for each time step. We obtained a steady state solution after 100,000 time steps. Figure 7 shows the stream function and vorticity contours for Re = 800. The flow separates at the step corner and vortices are formed downstream. We can observe the two vortices formed at the lower and upper wall. After reattachment of the upper wall eddy, the flow in the duct slowly recovers towards a fully developed flow. We compute the separation and reattachment points at L lower ≈ 6.1 for the lower wall separation zone, L upper ≈ 5.11 for the upper separation zone, where the separation begins at x ≈ 5.19. Our numerical findings show good agreement with other numerical methods for 2D computations [1,46]. In [1], the authors used a finite difference method and predicted separation lengths of L lower ≈ 6.0 and L upper ≈ 5.75, while [46] using the A Fluid Dynamics Analysis Program (FIDAP) code (Fluid dynamics International Inc., Evanston, IL, USA) predicted L lower ≈ 5.8 and the upper L upper ≈ 4.7.  We first discretize the flow domain with a uniform Cartesian grid. In our previous studies [29,35], we have shown that grid with resolution 31 × 901 is sufficient to compute a gridindependent numerical solution for Reynolds numbers up to Re = 800. As in the previous section, to highlight the efficiency of our scheme, we use a denser grid with resolution 121 × 3630, resulting in 439,230 nodes. We set the total time Ttotal = 250 (to ensure that the solution will reach steady state), and the time step dt = 10 −4 . The simulation terminates when the NRMSE of the time derivative of the stream function and vorticity field values in two successive time steps is less than 10 −6 . The time needed to create the grid was 0.023 s, and it takes 0.3 s to update the solution for each time step. We obtained a steady state solution after 100,000 time steps. Figure 7 shows the stream function and vorticity contours for Re = 800. The flow separates at the step corner and vortices are formed downstream. We can observe the two vortices formed at the lower and upper wall. After reattachment of the upper wall eddy, the flow in the duct slowly recovers towards a fully developed flow. We compute the separation and reattachment points at Llower ≈ 6.1 for the lower wall separation zone, Lupper ≈ 5.11 for the upper separation zone, where the separation begins at x ≈ 5.19. Our numerical findings show good agreement with other numerical methods for 2D computations [1,46]. In [1], the authors used a finite difference method and predicted separation lengths of Llower ≈ 6.0 and Lupper ≈ 5.75, while [46] using the A Fluid Dynamics Analysis Program (FIDAP) code (Fluid dynamics International Inc, Evanston, IL, USA) predicted Llower ≈ 5.8 and the upper Lupper ≈ 4.7.   Figure 8 shows the comparisons of the uand vvelocity components and vorticity values, between the present scheme and those obtained using the finite element method (FEM) [45], along the line x = 7 and x = 15 for Re = 800. Additionally, we computed the shear stress, defined as ∂u ∂y on the lower and upper wall, for increasingly refined nodal resolution. Figure 9a shows plots of the shear stress for the upper and lower wall, while Figure 9b displays the convergence of the upper wall shear stress when successively denser point clouds are used. We can see that shear stress converges and is accurate when compared to results of [45]. The accurate computation of spatial derivatives is a particular advantage of the strong form meshless method. In contrast, convergence of many finite element methods can be obtained only in terms of integral norms, and the point-wise convergence for the velocity gradient cannot be guaranteed even in the case of a smooth numerical solution [47].  Figure 8 shows the comparisons of the u-and v-velocity components and vorticity values, between the present scheme and those obtained using the finite element method (FEM) [45], along the line x = 7 and x = 15 for Re = 800. Additionally, we computed the shear stress, defined as on the lower and upper wall, for increasingly refined nodal resolution. Figure 9a shows plots of the shear stress for the upper and lower wall, while Figure 9b displays the convergence of the upper wall shear stress when successively denser point clouds are used. We can see that shear stress converges and is accurate when compared to results of [45]. The accurate computation of spatial derivatives is a particular advantage of the strong form meshless method. In contrast, convergence of many finite element methods can be obtained only in terms of integral norms, and the point-wise convergence for the velocity gradient cannot be guaranteed even in the case of a smooth numerical solution [47].     Figure 8 shows the comparisons of the u-and v-velocity components and vorticity values, between the present scheme and those obtained using the finite element method (FEM) [45], along the line x = 7 and x = 15 for Re = 800. Additionally, we computed the shear stress, defined as on the lower and upper wall, for increasingly refined nodal resolution. Figure 9a shows plots of the shear stress for the upper and lower wall, while Figure 9b displays the convergence of the upper wall shear stress when successively denser point clouds are used. We can see that shear stress converges and is accurate when compared to results of [45]. The accurate computation of spatial derivatives is a particular advantage of the strong form meshless method. In contrast, convergence of many finite element methods can be obtained only in terms of integral norms, and the point-wise convergence for the velocity gradient cannot be guaranteed even in the case of a smooth numerical solution [47].    Tables 8 and 9 list the velocity and vorticity field values at the at x = 7 and x = 15, respectively, for Re = 800.

Unbounded Flow Past a Cylinder
We examine the case of external flow past a circular cylinder in an unbounded domain [48]. The flow domain is large compared to the dimensions of the cylinder, as shown in Figure 10. The cylinder cross-section has a radius R c = 0.5 and is located at the origin O (0,0) of a square domain with dimensions −10 ≤ x ≤ 30 and −20 ≤ y ≤ 20.
The uniform flow is not perturbed near the inlet, which is far from the body. Therefore, we assume that the inflow velocity u inlet behaves like the potential flow u potential given as: (44) or in terms of the stream function: The uniform flow is not perturbed near the inlet, which is far from the body. Therefore, we assume that the inflow velocity uinlet behaves like the potential flow upotential given as: , −2 ( + ) (44) or in terms of the stream function: We represent the flow domain with a uniform Cartesian embedded grid, locally refined in the vicinity of the cylinder (Figure 11a). We use 402,068 nodes (0.5 s to create the nodal distribution) and we set Reynolds number to Re = 40. We use a time step of = 10 for the entire flow simulation.
It takes 0.3 s to update the solution for each time step, and we obtained a steady state solution after 20,000 time steps. The spatial derivatives (up to 2nd order) are computed in 4.24 s using a C++ code. Additionally, we solve the flow problem using an irregular nodal distribution with 221,211 nodes, locally refined in the vicinity of the cylinder We represent the flow domain with a uniform Cartesian embedded grid, locally refined in the vicinity of the cylinder (Figure 11a). We use 402,068 nodes (0.5 s to create the nodal distribution) and we set Reynolds number to Re = 40. We use a time step of dt = 10 −4 for the entire flow simulation. It takes 0.3 s to update the solution for each time step, and we obtained a steady state solution after 20,000 time steps. The spatial derivatives (up to 2nd order) are computed in 4.24 s using a C++ code. Additionally, we solve the flow problem using an irregular nodal distribution with 221,211 nodes, locally refined in the vicinity of the cylinder. We compute the following flow parameters: the pressure coefficient (Cp) on the body surface, the length (L) of the wake behind the body, the separation angle (θs), and the drag coefficient (CD) of the body. The drag and lift coefficient of the body are given by: where D is the characteristic length of the body, and are the unit normal vector in x-and ydirection, respectively, and is the total drag force acting on the body. The total drag force is given as: where , the pressure drag, can be determined from the flux of vorticity on the surface of the cylinder as: The friction drag may be computed from the vorticity on the surface of the body as: Table 5 lists the recirculation length (Lrec) of the wake behind the body, the separation angle (θs), and the drag coefficient (CD) for Re = 40.
The numerical findings obtained by the proposed scheme are listed in Table 10, and are in good agreement with those in the literature [48][49][50][51]. The stream function and vorticity contours around the cylinder are illustrated in Figure 12 (as in [48]).  We compute the following flow parameters: the pressure coefficient (C p ) on the body surface, the length (L) of the wake behind the body, the separation angle (θ s ), and the drag coefficient (C D ) of the body. The drag and lift coefficient of the body are given by: where D is the characteristic length of the body,ê x andê y are the unit normal vector in xand ydirection, respectively, and F b is the total drag force acting on the body. The total drag force is given as: where F p , the pressure drag, can be determined from the flux of vorticity on the surface of the cylinder as: The friction drag F f may be computed from the vorticity on the surface of the body as: withê θ = − sin θi + cos θj. Table 5 lists the recirculation length (L rec ) of the wake behind the body, the separation angle (θ s ), and the drag coefficient (C D ) for Re = 40. The numerical findings obtained by the proposed scheme are listed in Table 10, and are in good agreement with those in the literature [48][49][50][51]. The stream function and vorticity contours around the cylinder are illustrated in Figure 12 (as in [48]).

Numerical Results
In this section, we examine the applicability and reproducibility of the proposed scheme for several flow cases involving complex geometries.

Vortex Shedding Behind Cylinders
In the previous example (unbounded flow past a cylinder), we considered the flow past a cylinder in an unbounded domain. In this example, we study the vortex shedding behind a cylinder in a rectangular duct [52] (internal flow), as shown in Figure 13. The duct has length L = 2.  First, we consider a uniform Cartesian embedded grid with grid spacing h, defined by the average distance of the boundary nodes on the cylinder circumference. We conducted a grid independence analysis, applying successively refined Cartesian grids, starting from 34,999 (90 nodes on the circular boundary) up to 406,678 (720 nodes on the circular boundary) nodes.
The time needed to create the grid of 406,678 nodes is 0.8 s. We delete the nodes located close to the boundary (with distance less than 0.25 h), since they would increase the condition number of the

Numerical Results
In this section, we examine the applicability and reproducibility of the proposed scheme for several flow cases involving complex geometries.

Vortex Shedding Behind Cylinders
In the previous example (unbounded flow past a cylinder), we considered the flow past a cylinder in an unbounded domain. In this example, we study the vortex shedding behind a cylinder in a rectangular duct [52] (internal flow), as shown in Figure 13. The duct has length L = 2.

Numerical Results
In this section, we examine the applicability and reproducibility of the proposed scheme for several flow cases involving complex geometries.

Vortex Shedding Behind Cylinders
In the previous example (unbounded flow past a cylinder), we considered the flow past a cylinder in an unbounded domain. In this example, we study the vortex shedding behind a cylinder in a rectangular duct [52] (internal flow), as shown in Figure 13. The duct has length L = 2.  First, we consider a uniform Cartesian embedded grid with grid spacing h, defined by the average distance of the boundary nodes on the cylinder circumference. We conducted a grid independence analysis, applying successively refined Cartesian grids, starting from 34,999 (90 nodes on the circular boundary) up to 406,678 (720 nodes on the circular boundary) nodes.
The time needed to create the grid of 406,678 nodes is 0.8 s. We delete the nodes located close to the boundary (with distance less than 0.25 h), since they would increase the condition number of the First, we consider a uniform Cartesian embedded grid with grid spacing h, defined by the average distance of the boundary nodes on the cylinder circumference. We conducted a grid independence analysis, applying successively refined Cartesian grids, starting from 34,999 (90 nodes on the circular boundary) up to 406,678 (720 nodes on the circular boundary) nodes.
The time needed to create the grid of 406,678 nodes is 0.8 s. We delete the nodes located close to the boundary (with distance less than 0.25 h), since they would increase the condition number of the Vandermonde matrix. A Vandermonde matrix with too large condition number would affect the accuracy of the spatial derivative computation and consequently the overall precision of the numerical simulation. A grid independent solution is obtained with 203,890 nodes (360 nodes on the circular boundary). The critical time step computed for this grid resolution is δt critical = 2 × 10 −3 . In the simulations, we set the time step to dt = 10 −3 s. For each time step, it takes 0.38 s to update the solution. Figure 14a shows the stream function contours at time instances t = 2, 4 and 8 s, while Figure 14b plots the vorticity isocontours. Figure 15 plots the drag C D = computed with our meshless explicit scheme (as in the previous example). To highlight the applicability of our code for locally refined nodal distributions, we use a uniform Cartesian embedded grid, locally refined in the vicinity of the cylinder. As before, the average distance of the cylinder nodes dictates the nodal spacing h d of the Cartesian grid in the refined part of the domain. Nodes of the coarse grid have spacing h c = 2h, and nodes located close to the cylinder (with distance less than 0.25 h d ) are deleted since they affect the condition number of the Vandermonde matrix and decrease the accuracy of the numerical results. The flow domain is represented by 321,897 nodes, which ensure a grid independent solution. For the simulations, we set the time step to dt = 10 −3 s. Vandermonde matrix. A Vandermonde matrix with too large condition number would affect the accuracy of the spatial derivative computation and consequently the overall precision of the numerical simulation. A grid independent solution is obtained with 203,890 nodes (360 nodes on the circular boundary). The critical time step computed for this grid resolution is = 2 × 10 . In the simulations, we set the time step to dt = 10 −3 s. For each time step, it takes 0.38 s to update the solution. Figure 14a shows the stream function contours at time instances t = 2, 4 and 8 s, while Figure 14b plots the vorticity isocontours. Figure 15 plots the drag = • and lift = • coefficients, computed with our meshless explicit scheme (as in the previous example). To highlight the applicability of our code for locally refined nodal distributions, we use a uniform Cartesian embedded grid, locally refined in the vicinity of the cylinder. As before, the average distance of the cylinder nodes dictates the nodal spacing hd of the Cartesian grid in the refined part of the domain. Nodes of the coarse grid have spacing hc = 2h, and nodes located close to the cylinder (with distance less than 0.25 hd) are deleted since they affect the condition number of the Vandermonde matrix and decrease the accuracy of the numerical results. The flow domain is represented by 321,897 nodes, which ensure a grid independent solution. For the simulations, we set the time step to dt = 10 −3 s.   Additionally, to highlight the versatility of the proposed scheme, we examined the flow in the rectangular duct with multiple (seven in total) cylindrical obstacles. The length of the duct was increased to L = 4.4 m in order to ensure fully developed flow at the outlet. The boundary conditions applied are the same as before. We use Cartesian embedded and irregular nodal distributions ( Figure  16) to represent the flow domain. We conducted a comprehensive analysis in order to obtain a grid independent numerical solution. For the Cartesian grid we use 279,178 nodes (0.42 s to create the nodal distribution), while for the irregular point cloud we use 353,617 nodes (1.2 s to create the irregular nodal distribution). Both nodal distributions provided a grid independent numerical solution. The total time was set to Ttot = 8 s, the time step was set to dt = 10 −4 s. For each time step, the time needed to update the solution was 0.32 s. Figure 17 displays the stream function and vorticity contours for Re = 100 at time t = 1 and t = 2 s. Additionally, to highlight the versatility of the proposed scheme, we examined the flow in the rectangular duct with multiple (seven in total) cylindrical obstacles. The length of the duct was increased to L = 4.4 m in order to ensure fully developed flow at the outlet. The boundary conditions applied are the same as before. We use Cartesian embedded and irregular nodal distributions ( Figure 16) to represent the flow domain. Additionally, to highlight the versatility of the proposed scheme, we examined the flow in the rectangular duct with multiple (seven in total) cylindrical obstacles. The length of the duct was increased to L = 4.4 m in order to ensure fully developed flow at the outlet. The boundary conditions applied are the same as before. We use Cartesian embedded and irregular nodal distributions ( Figure  16) to represent the flow domain. We conducted a comprehensive analysis in order to obtain a grid independent numerical solution. For the Cartesian grid we use 279,178 nodes (0.42 s to create the nodal distribution), while for the irregular point cloud we use 353,617 nodes (1.2 s to create the irregular nodal distribution). Both nodal distributions provided a grid independent numerical solution. The total time was set to Ttot = 8 s, the time step was set to dt = 10 −4 s. For each time step, the time needed to update the solution was 0.32 s. Figure 17 displays the stream function and vorticity contours for Re = 100 at time t = 1 and t = 2 s. We conducted a comprehensive analysis in order to obtain a grid independent numerical solution.
For the Cartesian grid we use 279,178 nodes (0.42 s to create the nodal distribution), while for the irregular point cloud we use 353,617 nodes (1.2 s to create the irregular nodal distribution). Both nodal distributions provided a grid independent numerical solution. The total time was set to T tot = 8 s, the time step was set to dt = 10 −4 s. For each time step, the time needed to update the solution was 0.32 s. Figure 17 displays the stream function and vorticity contours for Re = 100 at time t = 1 and t = 2 s.
We further demonstrate the accuracy of the proposed scheme, by comparing our results with those computed using the finite element method (FEM). Flow equations were solved using the Incremental Pressure Correction Scheme (IPCS), which is an operator splitting method [53]. For the meshless solver, we use 353,617 nodes, generated using a triangular mesh generator, and we set the time step to dt = 10 −3 s. For the finite element model, we solve the flow equations using the generated mesh and FEniCS (https://fenicsproject.org). We compute the Normalized Root Mean Square Error (NRMSE), as defined in Equation (43), for both velocity components (u x , u y ), considering as gold standard the velocity field values calculated using FEniCS. Figure 18 shows the NRMSE for the velocity components for the time interval between t = 0 and t = 1 s. We further demonstrate the accuracy of the proposed scheme, by comparing our results with those computed using the finite element method (FEM). Flow equations were solved using the Incremental Pressure Correction Scheme (IPCS), which is an operator splitting method [53]. For the meshless solver, we use 353,617 nodes, generated using a triangular mesh generator, and we set the time step to = 10 s. For the finite element model, we solve the flow equations using the generated mesh and FEniCS (https://fenicsproject.org). We compute the Normalized Root Mean Square Error (NRMSE), as defined in Equation (43), for both velocity components (ux, uy), considering as gold standard the velocity field values calculated using FEniCS. Figure 18 shows the NRMSE for the velocity components for the time interval between t = 0 and t = 1 s.   We further demonstrate the accuracy of the proposed scheme, by comparing our results with those computed using the finite element method (FEM). Flow equations were solved using the Incremental Pressure Correction Scheme (IPCS), which is an operator splitting method [53]. For the meshless solver, we use 353,617 nodes, generated using a triangular mesh generator, and we set the time step to = 10 s. For the finite element model, we solve the flow equations using the generated mesh and FEniCS (https://fenicsproject.org). We compute the Normalized Root Mean Square Error (NRMSE), as defined in Equation (43), for both velocity components (ux, uy), considering as gold standard the velocity field values calculated using FEniCS. Figure 18 shows the NRMSE for the velocity components for the time interval between t = 0 and t = 1 s.

Flow in a Duct with Irregular Geometry
Finally, we consider two flow cases in the flow domains shown in Figure 19. The first one has an irregularly shaped obstacle downstream, while in the second the flow domain is split into two branches, forming a bifurcation.

Flow in a Duct with Irregular Geometry
Finally, we consider two flow cases in the flow domains shown in Figure 19. The first one has an irregularly shaped obstacle downstream, while in the second the flow domain is split into two branches, forming a bifurcation.

Flow in a Duct with Irregular Geometry
Finally, we consider two flow cases in the flow domains shown in Figure 19. The first one has an irregularly shaped obstacle downstream, while in the second the flow domain is split into two branches, forming a bifurcation.   First, we represent the flow domain with a uniform Cartesian embedded grid. For the irregularly shaped obstacle (bypass), we consider successively denser nodal distributions to obtain a grid independent numerical solution. We start with a grid spacing of h = 1/8200, h = 1/12,300 and h = 1/16,400 for the coarse, moderate and dense Cartesian embedded grids, resulting in 51,468, 114,387 and 201,468 nodes, respectively (for the denser grid (201,468 nodes) it takes 0.25 s to create the nodal distribution). The total time for the simulation was set to T tot = 30 s and the time step used was dt = 10 -4 s (for each time step the solution is computed in 0.35 s). Figure 21 shows the iso-contours for the stream function and vorticity field functions, at different time instances: t = 0.5 and t = 1 s, and for steady state (steady state is reached when, for two successive time instances t n+1 and t n for both vorticity and stream function field values, NRMSE/dt < 10 −8 ). First, we represent the flow domain with a uniform Cartesian embedded grid. For the irregularly shaped obstacle (bypass), we consider successively denser nodal distributions to obtain a grid independent numerical solution. We start with a grid spacing of h = 1/8200, h = 1/12,300 and h = 1/16,400 for the coarse, moderate and dense Cartesian embedded grids, resulting in 51,468, 114,387 and 201,468 nodes, respectively (for the denser grid (201,468 nodes) it takes 0.25 s to create the nodal distribution). The total time for the simulation was set to Ttot = 30 s and the time step used was dt = 10 -4 s (for each time step the solution is computed in 0.35 s). Figure 21 shows the iso-contours for the stream function and vorticity field functions, at different time instances: t = 0.5 and t = 1 s, and for steady state (steady state is reached when, for two successive time instances and for both vorticity and stream function field values, / < 10 ). We observe that two vortices are formed; the first is located at the inlet of the bypass, the second at the upper domain of the bypass which is moving towards the outlet of the flow domain. Figure  22a   We observe that two vortices are formed; the first is located at the inlet of the bypass, the second at the upper domain of the bypass which is moving towards the outlet of the flow domain. Figure 22a plots the uvelocity profile, computed at x = 0.07 m and x = 0.08 m at time t = 3 s for the coarse, moderate and dense nodal distributions, while Figure 22b plots the uvelocity profile, computed at x = 0.07 m and x = 0.08 m at time t = 3 s for the denser grid used.
Additionally, to highlight the applicability of the method to irregular nodal distributions, we solved the flow equations by representing the flow domain with nodes generated by a 2D triangular mesh generator (we use only the coordinates and not their connectivity). We used 392,122 nodes, which is higher than the number of nodes used in the Cartesian embedded grid. The numerical results obtained were compared against those obtained by the Cartesian grid and they were in excellent agreement. Additionally, to highlight the applicability of the method to irregular nodal distributions, we solved the flow equations by representing the flow domain with nodes generated by a 2D triangular mesh generator (we use only the coordinates and not their connectivity). We used 392,122 nodes, which is higher than the number of nodes used in the Cartesian embedded grid. The numerical results obtained were compared against those obtained by the Cartesian grid and they were in excellent agreement.
For the bifurcated geometry, we use both uniform Cartesian and irregular nodal distributions. In the case of Cartesian embedded grid, we obtain a grid independent numerical solution by using successively denser nodal distributions, with the denser one consisting of 678,967 nodes (for the fine grid it takes 0.68 s to create the nodal distribution). The total time for the simulation was set to Ttot = 30 s and the time step used is dt = 10 −4 s, which is smaller than the critical time that ensures stability for the explicit solver (the solution in each time step is computed in 0.52 s). The inflow velocity = 4 , with Um = 0.015 m/s results in Reynolds number Re = 212. Figure 23 shows the iso-contours for the stream function and vorticity field functions, at time t = 0.5 and t = 1 s and for steady state (steady state is reached when, for two successive time instances and for both vorticity and stream function field values, / < 10 ).  For the bifurcated geometry, we use both uniform Cartesian and irregular nodal distributions. In the case of Cartesian embedded grid, we obtain a grid independent numerical solution by using successively denser nodal distributions, with the denser one consisting of 678,967 nodes (for the fine grid it takes 0.68 s to create the nodal distribution). The total time for the simulation was set to T tot = 30 s and the time step used is dt = 10 −4 s, which is smaller than the critical time that ensures stability for the explicit solver (the solution in each time step is computed in 0.52 s). The inflow velocity U inlet = 4U m y H−y H 2 , with U m = 0.015 m/s results in Reynolds number Re = 212. Figure 23 shows the iso-contours for the stream function and vorticity field functions, at time t = 0.5 and t = 1 s and for steady state (steady state is reached when, for two successive time instances t n+1 and t n for both vorticity and stream function field values, NRMSE/dt < 10 −8 ). Additionally, to highlight the applicability of the method to irregular nodal distributions, we solved the flow equations by representing the flow domain with nodes generated by a 2D triangular mesh generator (we use only the coordinates and not their connectivity). We used 392,122 nodes, which is higher than the number of nodes used in the Cartesian embedded grid. The numerical results obtained were compared against those obtained by the Cartesian grid and they were in excellent agreement.
For the bifurcated geometry, we use both uniform Cartesian and irregular nodal distributions. In the case of Cartesian embedded grid, we obtain a grid independent numerical solution by using successively denser nodal distributions, with the denser one consisting of 678,967 nodes (for the fine grid it takes 0.68 s to create the nodal distribution). The total time for the simulation was set to Ttot = 30 s and the time step used is dt = 10 −4 s, which is smaller than the critical time that ensures stability for the explicit solver (the solution in each time step is computed in 0.52 s). The inflow velocity = 4 , with Um = 0.015 m/s results in Reynolds number Re = 212. Figure 23 shows the iso-contours for the stream function and vorticity field functions, at time t = 0.5 and t = 1 s and for steady state (steady state is reached when, for two successive time instances and for both vorticity and stream function field values, / < 10 ).  To further demonstrate the accuracy of the proposed scheme, we compare our numerical findings for both flow cases (bypass and bifurcation geometry) with those obtained using the finite element method. The flow domain is discretized using 392,122 and 678,967 nodes, respectively, and Taylor-Hood elements (as in the vortex shedding behind cylinders example). We use the IPCS flow solver to numerically solve the flow equations. We compare both numerical solutions, computing the NRMSE for the two velocity components (u x and u y ) on the vertices of the triangular elements (the same nodes are used in the meshless point collocation flow solver). Figure 24 shows the NRMSE for the velocity components (bypass and bifurcation flow examples) for the time interval between t = 0 and t = 1 s. To further demonstrate the accuracy of the proposed scheme, we compare our numerical findings for both flow cases (bypass and bifurcation geometry) with those obtained using the finite element method. The flow domain is discretized using 392,122 and 678,967 nodes, respectively, and Taylor-Hood elements (as in the vortex shedding behind cylinders example). We use the IPCS flow solver to numerically solve the flow equations. We compare both numerical solutions, computing the NRMSE for the two velocity components (ux and uy) on the vertices of the triangular elements (the same nodes are used in the meshless point collocation flow solver). Figure

Conclusions
In this contribution, we described a meshless point collocation algorithm for solving the stream function-vorticity formulation of N-S equations. We demonstrated our algorithm to work well in complex geometries like the ones shown in sections 4.1 (Vortex Shedding Behind Cylinders) and 4.2 (Flow in a Duct with Irregular Geometry). We demonstrate the accuracy of the proposed scheme by comparing our numerical results with those computed using the finite element method. An important advantage of our method is the ease and speed with which one can construct computational grids for flow domains with irregular shapes.
The meshless scheme based on DC PSE methods to compute spatial derivatives can be used in the case of Cartesian and Cartesian-embedded grids. The method works both efficiently and accurately for uniform Cartesian (embedded) grids and for irregular point clouds. We verified the accuracy of the proposed scheme through the following three benchmark problems: lid-driven cavity flow, backward-facing step and unbounded flow past a cylinder. The results were in good agreement with other published numerical studies.
Our proposed scheme offers several advantages over other commonly used methods: Computational efficiency, as demonstrated by time per iteration (Table 3) and total simulation time given for all examples in this paper. Our approach makes Direct Numerical Simulations (DNS) [54] possible (see lid driven cavity example with 1024 × 1024 grid resolution), even on personal computers. • Critical time step is easily computed with low computational cost.

•
Long critical time steps, comparable to those used in implicit schemes. The critical time step is related to the number of nodes located in the support domain and increases as the number of nodes increases (see Tables 1 and 2).

•
The imposition of Dirichlet boundary conditions is straightforward.

Conclusions
In this contribution, we described a meshless point collocation algorithm for solving the stream function-vorticity formulation of N-S equations. We demonstrated our algorithm to work well in complex geometries like the ones shown in Section 4.1 (Vortex Shedding Behind Cylinders) and 4.2 (Flow in a Duct with Irregular Geometry). We demonstrate the accuracy of the proposed scheme by comparing our numerical results with those computed using the finite element method. An important advantage of our method is the ease and speed with which one can construct computational grids for flow domains with irregular shapes.
The meshless scheme based on DC PSE methods to compute spatial derivatives can be used in the case of Cartesian and Cartesian-embedded grids. The method works both efficiently and accurately for uniform Cartesian (embedded) grids and for irregular point clouds. We verified the accuracy of the proposed scheme through the following three benchmark problems: lid-driven cavity flow, backward-facing step and unbounded flow past a cylinder. The results were in good agreement with other published numerical studies.
Our proposed scheme offers several advantages over other commonly used methods: Computational efficiency, as demonstrated by time per iteration (Table 3) and total simulation time given for all examples in this paper. Our approach makes Direct Numerical Simulations (DNS) [54] possible (see lid driven cavity example with 1024 × 1024 grid resolution), even on personal computers. • Critical time step is easily computed with low computational cost.

•
Long critical time steps, comparable to those used in implicit schemes. The critical time step is related to the number of nodes located in the support domain and increases as the number of nodes increases (see Tables 1 and 2).

•
The imposition of Dirichlet boundary conditions is straightforward.

•
Easy and straightforward way to impose vorticity boundary conditions using spatial derivatives computed using the DC PSE method (Equation (33))