An immersed boundary vector potential‐vorticity meshless solver of the incompressible Navier–Stokes equation

We present a strong form meshless solver for numerical solution of the nonstationary, incompressible, viscous Navier–Stokes equations in two (2D) and three dimensions (3D). We solve the flow equations in their stream function‐vorticity (in 2D) and vector potential‐vorticity (in 3D) formulation, by extending to 3D flows the boundary condition‐enforced immersed boundary method, originally introduced in the literature for 2D problems. We use a Cartesian grid, uniform or locally refined, to discretize the spatial domain. We apply an explicit time integration scheme to update the transient vorticity equations, and we solve the Poisson type equation for the stream function or vector potential field using the meshless point collocation method. Spatial derivatives of the unknown field functions are computed using the discretization‐corrected particle strength exchange method. We verify the accuracy of the proposed numerical scheme through commonly used benchmark and example problems. Excellent agreement with the data from the literature was achieved. The proposed method was shown to be very efficient, having relatively large critical time steps.


INTRODUCTION
The IB method was introduced by Peskin 1 for modeling blood flow through a beating heart. In the IB methods framework, flow equations are solved on a uniform Cartesian grid (Eulerian nodes), while the boundary of the solid is described/represented through a set of nodes (Lagrangian nodes) which, in general, do not coincide with the Cartesian grid points. Fluid and solid nodes communicate through volume forces, which are used to enforce the velocity boundary conditions and to ensure the conservation of linear momentum. Forces and torque are computed on both Eulerian and Lagrangian nodes, through interpolating-spreading operators, such as discrete Delta functions, 1 moving least squares (MLS), 2 and radial basis functions (RBFs). 3 A number of variants of Peskin's IB method have been developed and applied to various physical problems. A review on the different flavors of the IB method can be found in Reference 4. Therein, authors divide the existing IB methods into two groups. In the first group of methods (widely used in problems involving sharp boundaries and rigid objects), referred as "continuous forcing," 1,5,6 body forces are derived prior to the discretization step. The second group, termed as "discrete forcing" methods, [7][8][9][10][11][12][13][14][15] is based on a set of singular body forces, defined on the Eulerian fluid nodes, to enforce the desired boundary values.
Among the many flavors of IB methods, the boundary condition-enforced immersed boundary (BCE-IB) method 16 is used in the present study since it accurately imposes the prescribed boundary conditions (Dirichlet and Neumann) on the Lagrangian nodes. Different approaches developed in the framework of IB methods, such as the penalty force method 17 and the direct forcing method, appear to have certain drawbacks. The penalty force method requires the input of spring stiffness parameter, which is arbitrarily defined, case dependent and in many cases induces stability problems. 18 Moreover, the penalty force method and the direct forcing method do not exactly impose the no-slip boundary condition, which may result in unphysical flow penetration through rigid boundaries. 19 Furthermore, both penalty and direct forcing methods produce unphysical (temporally oscillating) pressure fields when simulating flows around moving/deforming bodies with complex geometrical shapes.
Attempts to eliminate the computational burden of body-fitted mesh generation have resulted in the development of meshless methods (MMs). In MMs, the spatial (flow) domain is represented by a set of nodes/points (randomly or uniformly distributed over the domain), without any predefined connectivity between them. Various strong and weak form MMs have been applied to numerically solve the Navier-Stokes (N-S) equations. One of the most popular MM is the meshless local Petrov-Galerkin (MLPG) method, which has been used to solve the N-S equations in primitive variables, 20 velocity-vorticity 21,22 and stream function-vorticity 23 formulation. The local boundary integral method (LBIE) has been successfully applied to solve incompressible flows, 24 while the multi-region LBIE has been combined with the RBFs interpolation method. 25 The indirect/integrated radial-basis-function network (IRBFN) method 26,27 has been used to solve transient partial differential equations (PDEs) governing fluid flow problems. Finally, the meshless point collocation method has been used in References 28-32 to numerically solve the N-S equations, along with the velocity-and pressure-correction method.
In this contribution, we develop a strong form meshless point collocation method along with the BCE-IB method 16 to numerically solve the nonstationary, incompressible N-S equations in complex geometries, in two and three dimensions. Our method relies on the ability of the MM to accurately and efficiently solve the N-S equations on (uniform or locally refined) Cartesian grids, and on the flexibility of the IB methods to treat complex geometries. The main advantage of the BCE-IB method is that it satisfies accurately both the governing equations and boundary conditions by correcting the velocity field to account for the IB. The velocity correction applies implicitly such that the velocity on the IB (Lagrangian points) interpolated from the corrected velocity values computed on the mesh nodes (Eulerian nodes) accurately satisfies the prescribed velocity boundary conditions. We numerically solve the flow equations in their stream function-vorticity formulation in 2D, and vector potential-vorticity formulation in 3D. The method applies to external and internal flows, as well as to moving objects (convex and non-convex). In strong form meshless point collocation methods, the spatial domain is represented by a set of nodes, without any predefined connectivity. Therefore, the MM has features similar to the finite differences (FD) method in terms of representation of the spatial domain, and at the same time has some unique and distinctive advantages over FD. FD methods developed for adaptive mesh refinement 33 are far more complex than MMs to implement. Additionally, when the Euler explicit time integration scheme is used, the inherent upwinding nature of MMs (neighboring nodes are used to compute spatial derivatives) allows for longer time steps. 34 Even for extremely fine grids and high Reynolds numbers, the critical time step is longer compared to that in FD methods, and it is comparable to the time-step used in implicit time integration schemes. 34 This particular feature of the proposed MM makes it attractive for flow problems with high Reynolds numbers, where Direct Numerical Simulations can be used. 34 We use the stream function-vorticity formulation (in 2D) and the vorticity-vector potential formulation (in 3D) to avoid the calculation of pressure. In the former case, the number of variables decreases from three (u, v, p) to two ( , ), while in the latter increases, from four (u, v, w, p) to six ( x , y , z , Ψ x , Ψ y , Ψ z ). Nevertheless, the reduced computational cost and the physics related to flows involving vorticity (vorticity captures more accurately the flow physics), highlight the need for the vorticity formulation of the governing equations (details in Reference 34). The numerical solution of N-S equations in primitive variables formulation requires solution in the entire physical space, which increases the computational cost. In the vorticity formulation, wall vorticity can be calculated with high spatial and temporal accuracy (this is not the case for fractional step methods where accuracy near the walls can be reduced). In many applications where wall vorticity is the flow parameter (it can capture the physics of viscous flows), vorticity based methods will result in a higher numerical accuracy for a given spatial and temporal resolution.
This article is structured as follows. In Section 2, we present the meshless point collocation BCE-IB method along with the governing equations for the stream function-vorticity and vorticity-vector potential formulation based on the meshless point collocation method, and we briefly describe the discretization-correction particle strength exchange (DC PSE) differentiation method. In Section 3, we verify the proposed scheme through comparison with benchmark problems, while in Section 4, numerical examples are presented to demonstrate the accuracy and the ease of use of the proposed scheme. Finally, in Section 5, we present our conclusions.

THE PROPOSED ALGORITHM
We solve the flow equations in their stream function-vorticity (in 2D) and vector potential-vorticity (in 3D) formulation, by extending to 3D flows the BCE-IB method, originally introduced in the literature for 2D problems.

BCE-IB method for the stream function-vorticity N-S equations
Consider a rectangular computational domain (named Eulerian domain) Ω, which contains an IB (the Lagrangian domain), described in the form of a closed curve Ω L , as shown in Figure 1. The IB introduces an additional term in the flow equations. This extra term ( * ) in the stream function equation (Equation 2) is considered as a vorticity correction. The governing equations, in their stream function-vorticity formulation, that describe the fluid motion with the use of IB method, are written as: F I G U R E 1 A two-dimensional spatial domain and the immersed boundary (black dots) subject to the following boundary and initial conditions: with u = (u, v), , and being the velocity, stream function, and vorticity field, respectively, while and are the fluid density and dynamic viscosity. U B is the prescribed velocity on the Lagrangian points. The extra term * in Equation (2) is the source term related to the IB method, while x = (x, y) and X(s) (the immersed object can be viewed as a parametric curve X(s), with s being the parameter) are the Eulerian and Lagrangian points, respectively. The velocity U(X(s), t) at the Lagrangian points is interpolated from the velocity u(x, t) at the fluid (Eulerian) points as where (x − X(s)) is the Dirac delta function, representing the interaction between the fluid and the IB.
In our analysis, we use the Euler explicit time integration scheme to update the vorticity field (Equation 1). By applying the Euler explicit scheme to Equation (1), we obtain the predicted vorticitỹ(the variables with a tilde placed above the symbol are used to correct variables at the next time step) from the stream function n and vorticity n values at the previous time step (we use the notation n+1 = (t n+1 ) and n+1 = (t n+1 ), where t n+1 = t n + t), through the equatioñ Following, the predicted stream functioñfield is computed by solving the Poisson equation 2̃= −̃ (8) and the predicted velocityũ = (ũ,ṽ) field as:ũ =ỹ , We correct the velocity field, such that the prescribed velocity boundary conditions on the IB are satisfied. The corrected/updated velocity field u n+1 is written as where u * = ( u * , v * ) is the velocity correction so that the corrected velocity u n+1 satisfies the prescribed velocity U B on the Lagrangian points. Accordingly, the corrected/updated vorticity value n+1 is obtained as: with * being the vorticity correction. Since the updated vorticity field n+1 is defined as the vorticity correction * is computed as: Once the velocity correction is calculated, the updated velocity u n+1 and vorticity n+1 fields are calculated from Equations (11) and (12), respectively. The velocity u * is computed such that when the updated velocity u n+1 (x, t) (computed on the Eulerian nodes) is interpolated (U(X(s))) on the Lagrangian points, the updated velocity u n+1 (x, t) equals the prescribed velocity U B . Therefore, the velocity correction is distributed from the virtual boundary flux U through the delta function interpolation The IB is represented by a set of Lagrangian points X i = (X i , Y i ) , i = 1, 2, … , M and the flow domain by the Eulerian nodes x j = ( x j , y j ) , j = 1, 2, … , N in a Cartesian grid (uniform of locally refined) with grid spacing h x = h y = h. In our analysis, the Eulerian grid in the vicinity of the immersed object has a uniform Cartesian grid structure with grid spacing h. Furthermore, the Dirac function (x − X(s)) is approximated by a continuous kernel distribution D ( given as with the kernel (r) proposed by Lai and Peskin 17 Equation (15) is written in a discrete form where u * (x i , t) is the velocity correction at the ith Eulerian point, and Δs j is the jth boundary segment length. The unknowns in Equation (18) are the fictitious boundary fluxes U(X j , t). To satisfy the velocity boundary conditions on the Lagrangian nodes, the velocity U(X j , t) at the Lagrangian points, projected from the corrected fluid velocity field u n+1 (x, t) through the smooth delta function ij , must be equal to the prescribed boundary velocity U B on the Lagrangian points In Equation (19), u n+1 (x i , t) is the corrected velocity at the ith Eulerian point, expressed as the sum of the tentative velocityũ and velocity correction u * , as shown in Equation (11). Substituting Equations (11) and (18) into Equation (19), we get The equation system for the virtual boundary flux is now formed, and can be rewritten in matrix form as with R, F, and B defined as In the velocity correction procedure described above, the velocity correction u * is computed implicitly such that the velocity U(X(s), t) at the boundary point interpolated from the corrected velocity field u n+1 (x, t) equals the specified boundary velocity. The vorticity correction (Equation 14) is obtained by using the meshless DC PSE scheme (see Section 2.3) to compute the spatial derivatives. The DC PSE scheme is one of the most accurate methods for the computation of spatial derivatives. 35,36

BCE-IB method for the vector potential-vorticity N-S equations
The vector potential-vorticity formulation is the extension in three dimensions of the stream function-vorticity formulation. We consider a cubical (or rectangular) computational domain (Eulerian domain) Ω, which comprises an immersed body described in the form of a surface Ω L , as shown in Figure 2. The IB is facilitated by the introduction of an additional term * in the vector potential equation (Equation 26).
F I G U R E 2 A three-dimensional spatial domain and the immersed boundary (black dots) The flow equations are given as: where = ( x , y , z ) is the vorticity field, and = (Ψ x , Ψ y , Ψ z ) is the vector potential, defined such that the velocity u is written as: We update the vorticity field using the explicit Euler time integration scheme. By applying the explicit Euler scheme to Equation (25), we obtain the predicted vorticity components̃i, i = x, y, z through the velocity components u n i and vorticity n i values at the previous time step, through the equatioñ Following, we compute the predicted vector potential̃= (Ψ x ,Ψ y ,Ψ z ) field value by numerically solving the three Poisson equations of the vector potential components From the predicted vector-potential̃we compute the predicted velocityũ = (ũ x ,ũ y ,ũ z ) using Equation (27). We correct the predicted velocity fieldũ such that the prescribed velocity boundary conditions on the IB are satisfied. The corrected/updated velocity field u n+1 components are written as where u * = ( u * x , u * y , u * z ) are the components of the velocity correction. The corrected/updated vorticity value n+1 is given as where * = ( * x , * y , * z ) is the vorticity correction. The updated vorticity field n+1 = ( n+1 x , n+1 y , n+1 z ) is given as̃+ * = × u (32) and the vorticity correction * is given as * = × u * . The velocity correction is computed, the updated velocity u n+1 and vorticity n+1 fields are calculated from Equations (30) and (31), respectively. The IB is represented by a set of Lagrangian points X j (j = 1, 2, … , M), and the fluid by a set of Eulerian points x i (i = 1, 2, … , N) in a fixed uniform Cartesian grid, with mesh spacing h x = h y = h z = h. The Dirac function (x − X(s)) is approximated by a continuous kernel distribution D(x i − X j ) defined as The equations formulated for two-dimensional problems apply also for BE-IB method in three dimensions Equations (15)- (22) are also used in three-dimensions, by using h 3 in place of h 2 , and Δs j being the jth surface (instead of segment) element.

Spatial derivatives computation using DC PSE method
The DC PSE method computes spatial derivatives using a set of nodes that are arbitrarily distributed over the spatial domain. DC PSE originated as a Lagrangian particle-based numerical method, 36 and formulated as a correction of the particle strength exchange (PSE) operators. 37,38 PSE has been used for the evaluation of spatial derivatives of any degree of a field function (sufficiently smooth) discretized over scattered collocation points. To solve PDEs using the DC PSE MM, Bourantas et al. 35 reformulated the Lagrangian based DC PSE method to work in the Eulerian framework. The PSE operators are derived in two steps. First, by constructing an integral operator that satisfies continuous moment conditions, and that approximates a spatial derivative to a defined order of accuracy in the continuous domain. Second, by discretizing the integral operator over the particle positions using the mid-point quadrature method. A drawback from this two-step procedure is the introduction of an overlap condition for the PSE operators. The overlap condition requires that for the operator to be consistent, the inter-particle spacing h and the width of the operator kernel have to satisfy h ≤ a r , 0 < a < 1, r ≥ 1, where r depends on the accuracy of the operator. This results in a large number of particles being required for small kernel sizes . For the DC PSE method, the overlap condition can be relaxed by directly satisfying discrete moment conditions over the collocation points, requiring only that h = a independent on the order of convergence r. The discrete moment conditions for DC PSE can be thought of as a direct analog to the continuous moment conditions for PSE. The DC PSE operators overcome the limitation in PSE, where the error from operator discretization dominates the overall order of accuracy of the operator as prescribed by the moment conditions.
Here we briefly review the DC PSE operators for strong form formulations (excluding collocation point volumes) in 2D and how to construct them. For simplicity of example, we present the DC PSE formulation in 2D with standard kernel functions and form. The formulation of the DC PSE operators in n dimensions with arbitrary kernel functions is straightforward and the reader is directed to Reference 36 for an in-depth analysis. We begin by considering a differential operator, of arbitrary order, for a sufficiently smooth field where m and n are integers that determine the order of the differential operator. We define the DC PSE operator for approximating the spatial derivative D m,n f (x p ) as: where (x) is a spatially dependent scaling or resolution function, (x, (x)) a kernel function, and N(x p ) is the set of points in the support of the kernel function. The sign in Equation (35) is positive for m + n odd, and negative if even. The form of the operator in Equation (35), including the sign change, is chosen to match with Reference 36 for similarity to the original PSE operators. 38 However, the DC PSE formulation outlined here can also be applied in general to operators in the form . We construct the DC PSE operators such that we decrease the average spacing between nodes, h(x p ) → 0, the operator converges to the spatial derivative D m,n f (x p ) with an asymptotic rate r: where it is convenient to explicitly define the component-wise average neighbor spacing as h( where N is the number of nodes in the support of x p . Therefore, we need to find a kernel function (x) and a scaling relation (x p ) as to satisfy Equation (36). To achieve this, we replace the terms f (x q ) in Equation (35) with their Taylor series expansions around the point x p . This substitution gives: ) .
It is convenient to rewrite Equation (37) in the form: where ) .
We call Z i,j (x p ) the discrete moments. Now if we restrict the scaling parameter (x p ) to converge at the same rate as the average spacing between points h(x p ), that is we find that the discrete moments is (1), through normalization of the function argument.
These qualities are the motivation for the form of the normalized kernel function. Note that Equation (40) is a constraint than the original overlap condition h ≤ a r , 0 < a < 1, r ≥ 1 for the PSE method. Given Equation (40), the convergence behavior of Equation (38) is determined by the coefficients of the terms (x p ) i+j−m−n D i,j . We make the DC PSE operator satisfy Equation (38) by ensuring the coefficients of the (m, n) term is 1, and all other coefficients of ( This results in the following set of conditions for the discrete moments, where min is 1 if i + j is odd and 0 if even. This is due to the zeroth moment Z 0,0 canceling out for odd i + j. The choice of the factor (x p ) −m−n in Equation (38) simply acts as to simplify the expression of the moment conditions. For the kernel function (x) to satisfy the l = 6 conditions given in Equation (41) for arbitrary neighborhood node distributions, the operator must have l degrees of freedom. This leads to the requirement that the support N(x) of the kernel function has to include at least l nodes. In this article, as in Reference 36, we use kernel functions for the operator of the form This is a monomial basis multiplied by an exponential window function, where r c sets the kernel support and the a i,j are scalars to be determined to satisfy the moment conditions in Equation (41). The cut-off radius r c should be set to include at least l collocation nodes in the support N(x). To construct the operator Q m,n f (x p ) at node x p , the coefficients are found by solving a linear system of equations from Equations (41) and (42). With our choice of kernel function we have, where p(x) = {p 1 (x), p 2 (x), … , p l (x)} and a(x) are vectors of the terms in the monomial basis and of their coefficients in Equation (42), respectively. The normalization factor (x p ) −2 is absent, since it cancels out through the coefficients a(x p ). Using this formulation, the operator system becomes easy to obtain. For example, if we set r = 2 and approximate the first spatial derivative in the x direction, D 1,0 , we have l = 6 moment conditions and the monomial basis is p(x) = {1, x, y, yx, x 2 , y 2 }. The linear system for the kernel coefficients then is: where The scalar k > l is the number of nodes in the support of the operator, l is the number of moment conditions to be satisfied, and V(x p ) is the Vandermonde matrix constructed from the monomial basis p(x p ). E(x p ) is a diagonal matrix containing the square roots of the values of the exponential window function at the neighboring nodes in the operator support.
, the set of vectors pointing from all neighboring nodes in the support of x p to x p . So then explicitly Once the matrix A(x p ) is formed at each node x p , the linear systems can be solved and the coefficients a(x p ) used to construct the DC PSE operators at each node as in Equation (43) are computed. It is worth noting that the matrix A(x p ) depends on the number of moment conditions l and the local distribution of nodes in N(x). Therefore, if the system in Equation (44) is solved using a decomposition (such as LU), or inversion, of A(x p ) this solution can be reused for multiple right-hand sides, that is, for different derivative operators. The matrix A, similar to the moment matrix in moving least squares, provides information about the spatial distribution of the collocation nodes around the center point x p . The invertibility of A depends entirely on the invertibility of the Vandermode matrix V(x), due to E(x) being a diagonal positive-definite matrix. Finally, E affects the condition number of the resulting linear system.

Poisson's equation solution using direct and iterative solvers
To compute the stream function (in 2D) and vector potential (3D) field values, we numerically solve elliptic (Poisson) type PDEs. Iterative and direct solvers for elliptic type PDEs are well established. Direct solvers are robust and easy to use, but they can be computationally costly and memory intensive when the system of equations to be solved is large (a few million degrees of freedom). On the other hand, a decisive advantage of iterative solvers is their low memory usage, significantly less than a direct solver. Unfortunately, iterative solvers are less robust than direct solvers, and may fail to converge to an accurate numerical solution. Nevertheless, it is computationally more efficient to use iterative solvers to solve the Poisson equation when a large number of nodes are used. Many iterative solvers have been developed for Poisson type equations, 39 such as Gauss-Seidel or successive over relaxation (SOR). These solvers are accurate but computationally demanding, since they are characterized by slow convergence ( ( N 2 ) iterations are needed for convergence-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.
In the present study, we use either a direct or iterative solver depending on the number of degrees of freedom used. The discrete Laplace operator in the Poisson equation is accurately discretized using the DC PSE method, which performs accurately on uniform and locally refined Cartesian grids. The grid is fixed in time so a LU factorization is precomputed, and the computational cost is significantly reduced. For up to 1 million nodes, we use a direct solver to numerically solve the Poisson equation, since the factorization is fast and requires low memory. For larger number of nodes, we use the biconjugate gradient stabilized (BiCGSTAB) iterative solvers, since the linear systems are always positive definite, diagonally dominant and symmetric.

Critical time step for stable explicit time integration
The Euler explicit time integration scheme is conditionally stable. Therefore, the time step used should not be longer than the critical time step that ensures solution stability. To compute the critical time step, we rewrite Equation (1) as and Equation (25) as Equations (50)-(53) can be written in matrix form as Equation (54) is valid for both 2D and 3D flows. n+1 is a column vector of dimension N × 1 (N is the number of nodes), which in 2D corresponds to the vorticity field , and in 3D to the vorticity components ( x , y , z ) in the notation used, we omit the bold in vector variables. Matrix A in 2D is defined as and in 3D with i, j = x, y, z and matrices K, K i , and L, L VP i representing linear operators approximated using DC PSE as The explicit integration scheme defined by Equation (54) is stable if 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 Matrix A or A i is determined by the Laplacian operator L or L VP i , consisting of second order derivatives, and by the advection operator K or K i , which includes only first order derivatives. The relative weight of the two operators on the structure of matrix A or A i is dictated by discretization and velocity values (velocity is related to stream function and vector potential field values); a more refined discretization leads to a higher weight of operator L or L VP i (diffusion) and lower influence of operator K or K i (convection), while higher Reynolds number and higher velocity values result in higher influence of operator K or K i . When the nodal discretization is not sufficiently refined, the higher weight of the operator K or K i can lead to eigenvalues with a positive real part; in such case explicit time integration is not possible and discretization has to be refined. When the problem discretization includes a large number of nodes, matrix A becomes very large, and calculating the eigenvalues is not practical. However, the matrix A is diagonally dominant, with eigenvalues distributed close to the real axis. Therefore, we can conservatively estimate the critical time step using Equations (62) or (63) from the upper bound of the maximum eigenvalue of matrix A or A i computed using the Gershgorin circle theorem. 41 Given the composition of matrix A and A i , the upper bound of the maximum eigenvalue can be computed without assembling the matrix: Equation (64) shows that the eigenvalues of matrix A or A i depend on the spatial derivatives of the field variables computed using the DC PSE method. We have shown in Reference 34 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 3D geometries) in the support domain of each node, compared with the critical time step computed using a small number of nodes (around 30). 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 inherent upwind nature of MMs. Meshless methods are local support methods that use neighboring nodes to compute spatial derivatives. Furthermore, since matrices A and A i involve 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 CFL condition . 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.

Surface area of Lagrangian points in BCE-IB method
In 3D flow simulations, the surface area Δs i that is assigned to each Lagrangian node, is not as straightforward to define and compute as in two dimensions. The number and distribution of the Lagrangian points over the IB, affect the computation of ΔS i , and the accuracy of the IB method.
In 2D, the distribution and number of Lagrangian points are typically described through the ratio of ds h , with ds being the spacing of the Lagrangian points and h the Eulerian grid spacing. High ratio of ds h will cause unphysical fluid leakage when solving N-S equations (the spreading from Lagrangian to Eulerian points and vice-versa is not accurate). Low ratio will increase the computational cost (in this case rows in matrix R in Equation (22) will be linearly dependent and the matrix will become ill-conditioned and difficult to invert). There are many studies which provide recommendations regarding the ratio ds h ; Uhlmann 42 suggests that Lagrangian points spacing ds should be equal to the Eulerian grids resolution h (ds∕h = 1)), Su et al. 43 and Kang and Hassan 44 used ds h = 1.5 (Lagrangian points spacing is larger than the Eulerian grids resolution), while Silva et al. 45 reported on ds h ≤ 0.9 (Lagrangian points spacing is less than the Eulerian grids resolution). In 3D, the surface area Δs i appears as a key feature of the proposed scheme. In this study, we demonstrate an easy to apply method to compute Δs i . The area of each surface element Δs i should be close to ah 2 , with 0.1 ≤ a ≤ 4 (just like having in 2D 0.3 ≤ ds h ≤ 0.9, with h being the spacing of the Eulerian grid. In this study, we uniformly distribute the Lagrangian nodes, in order to assign the same value Δs to all Lagrangian nodes. To obtain a distribution of Lagrangian nodes with uniform surface element area, we generate a surface mesh (the nodes are the Lagrangian nodes used in the IB method) with element edges having lengths with small standard deviation from the mean value (practically we try to have element edges of equal length). The area element Δs for each Lagrangian node is the computed mean value of the facets area.

Solution procedure
The algorithmic procedure of the proposed IB method is summarized below. To update the solution from time instance t n to t n+1 : 1. Compute the tentative vorticitỹand stream functioñfield values by solving Equations (7) and (8). We use u n and n as initial flow values. 2. Compute the tentative velocity u by substituting̃in Equations (9) and (10). 3. Compute matrix R using Equation (22). 4. Compute virtual boundary fluxes U i at the Lagrangian points and then the velocity correction u * using Equation (15). 5. Compute the corrected/updated velocity u n+1 at the Eulerian points using Equation (11). 6. Compute the vorticity correction * using Equation (14). 7. Correct the vorticity at Eulerian points using Equation (12). 8. Use the corrected vorticity and velocity as the initial conditions, and repeat steps 1-7 for the computation of the next time increment. The process continues until a converged solution is achieved (steady case) or the given time is reached (unsteady case).
For the IB vector potential-vorticity formulation method the algorithm can be written as:  (Equation 14). 8. Compute the updated/corrected vorticity at Eulerian points using n+1 =̃+ * . 9. Use the corrected vorticity and velocity as the initial conditions, and repeat steps 1 to 8 for the computation of the next time increment. The process continues until a converged solution is achieved (steady case) or the given time is reached (unsteady case).

VERIFICATION OF THE PROPOSED ALGORITHM
In this section, we demonstrate the accuracy of the proposed method on both external (including moving objects) and internal flow problems. Our numerical results are compared against analytical solutions where available, and numerical results obtained using the well-established FE method. In all flow cases considered, we use the Euler explicit time integration scheme and report on the critical time step computed using the Gershgorin theorem. 41 We consider three benchmark problems: the unbounded flow past a cylinder; the in-line oscillating cylinder in fluid at rest, where analytical and experimental data are available; and the 3D lid-driven cavity flow with fixed sphere, where our numerical findings are compared against the numerical results using the ANSYS CFX commercial flow solver.

Unbounded flow past a cylinder
In this example, we examine the unbounded external flow past a cylinder with circular cross-section. 46 All lengths are non-dimensionalized by the cylinder diameter. The flow domain is large compared to the dimensions of the cylinder (Figure 3). The cylinder, with radius R c = 0.5, is located at point O(0, 0) of a square domain with dimensions −10 ≤ x ≤ 30 and −20 ≤ y ≤ 20. The flow is not perturbed near the inlet, which is located far from the cylinder. Therefore, the inflow velocity u inlet behaves like the potential flow u potential given as: or in terms of the stream function  The boundary conditions for the stream function and vorticity are shown in Figure 3. The nodal distribution that represents the flow domain is used to update the stream function and vorticity field values. We use a Cartesian grid that is locally refined in the vicinity of the cylinder as shown in Figure 4A. The adaptive mesh used in this example consists of five subdivisions Ω 1 -Ω 5 (subregions), each one having half the grid spacing of the previous one, and a locally refined layer in the vicinity of the IB.
We use successively denser grids to obtain a grid independent solution. The grids are described in Table 1. The initial resolution is denoted by h and the final one, after six refinements, by h f . The cylinder is represented by equally spaced Lagrangian nodes, such that the nodal spacing h f around the cylinder is h f ≈ Δs = 2 R c 180 . The spatial derivatives (up to second order) were computed using DC PSE method.
First, we consider the case with Reynolds number Re = 40 where the flow reaches a steady state solution, and the wake of the cylinder is symmetric and steady. We compute the critical time step (that ensures stability) of the explicit solver on the fly during the flow simulation using the Gershgorin circle theorem. The critical time step Δt critical and the time step Δt are listed in Table 1. We use a smaller time step Δt than the critical timer step Δt critical in our simulation to ensure the stability of the numerical solution. For Re = 40, the flow regime is characterized by a number of non-dimensional parameters (see Figure 5), including the length L sep of the recirculation zone, the distance (a) from the cylinder to the center of the wake vortex, the distance (b) between the centers of the wake vortices and the separation angle measured from the x-axis (see Figure 5).
To verify the accuracy of the proposed scheme, we compare the drag coefficient computed using the proposed method with results published in the literature. [46][47][48][49] The drag coefficient (C D ) of the immersed body is given by: Our numerical results (Table 2) at steady state are in close agreement with numerical results in the literature. [46][47][48][49] Figure 6 shows the drag coefficient C D for the different grid configurations. It can be seen that numerical solution obtained using 206,877 nodes ensures convergence of the drag coefficient. Figure 7 displays the instantaneous streamlines and the corresponding vorticity contours around the cylinder. We further demonstrate the accuracy of the proposed scheme for higher values of the Reynolds number where the flow structures become more complicated. 50 In the simulations, at Re = 550 we use a cylinder with radius R c = 1 (instead of R c = 0.5 for Re = 40) and the finest grid (856,965 Eulerian nodes and 3420 Lagrangian nodes) to study the convergence and accuracy of the numerical solution. We set the time step Δt = 2.5 × 10 −5 to be shorter than the critical time step Δt critical = 9 × 10 −5 to ensure stability. Figure 8 displays the streamline history of the flow field for Re = 550 at different time instances. A small secondary flow region appears at t = 3, along with a diving streamline of the flow. 50,51 In Figure 9, we demonstrate the development of the flow using vorticity contour plots. In the vorticity plots, a secondary vortex forms and becomes visible at t = 1. The evolution of the secondary vortex, which initially is confined, is mainly affected by the dynamics of the primary vortex. As the secondary vortex starts to grow, it penetrates the primary vortex. The primary vortex moves away from the immersed body, starts to diffuse and blocks the growth of the secondary vortex. Figures 8 and 9 were produced using a grid of 488,315 Eulerian nodes and 2900 equally spaced Lagrangian nodes on the cylinder boundary (the difference in the velocity field between the two finest grids was almost indistinguishable). Figure 10 shows the evolution of the drag coefficient versus time.
An important aspect of the explicit solver is to compute a critical time step that ensures stability. As stated in Reference 34, the DC PSE MM allows large time steps for explicit time integration schemes. In our simulations, we compute the critical time step using the Gershgorin circle theorem. We monitor the critical time step during the simulation to ensure that the time step used is always lower than the critical time step (see Figure 10B). We compute the critical time step

In-line oscillating cylinder in fluid at rest
In this example, we verify the accuracy of the proposed scheme on a moving boundary flow problem. We consider the inline oscillating circular flow problem, which has been previously investigated experimentally 52 and numerically. 15,[53][54][55] The problem geometry and boundary conditions are shown in Figure 11. The size of the computational domain is 14D × 14D, with the cylinder of diameter D located initially at the center. We apply no-slip boundary conditions on the boundaries of the flow domain. The trajectory of the cylinder center in the x-direction (oscillation direction) is given by the harmonic oscillation .
The oscillating amplitude is A = to the left, two thin boundary layers develop on the upper and lower sides ( Figure 12B), which later separate into two counter-rotating vortices with equal magnitude. Vortex generation stops as the body reaches its extreme left location ( Figure 12A). The cylinder moves backwards, and a mirroring process takes place on its right side. During this part of the cycle, the vortex pair generated earlier is split and pushed away by the cylinder (Figure 12C). We evaluate the accuracy of the proposed scheme, by quantitatively comparing our numerical findings with the experimental results reported by Dutsch et al. 52 Figure 13 displays profiles of the velocity components (u, v) at four locations (x = −0.6D, 0D, 0.6D, 1.2D) and three different phase angles (2 ft = 180 • , 210 • , 330 • ). There is satisfactory agreement of our numerical results with the experimental data. 52

3D lid-driven cavity flow with fixed sphere
We consider the 3D lid-driven cavity flow problem with a fixed sphere located inside a unit cubical domain. A sphere, with radius R c = 0.2, is immersed within the cube, located at its center. We numerically solve the flow equations using the vector potential-vorticity formulation. [56][57][58][59][60] The flow field is initially at rest. We apply no-slip boundary conditions (u = v = w = 0) on all boundaries, except on the top wall (z = 1), which moves with constant velocity (U = (1, 0, 0)). The boundary conditions for the vector potential field = ( while vorticity boundary conditions are computed using Equation (5) We discretize the flow domain using successively finer uniform Cartesian grids with spacing 41 3 = 68,921 nodes, 61 3 = 226,981 nodes, and 81 3 = 531,441 nodes, to achieve a grid independent solution. The results on the two finest grids are almost indistinguishable, indicating grid independence. Therefore, in our simulations we use a grid spacing of h = 1 60 , resulting in 226,981 nodes. The surface of the sphere was discretized using 1025 points. The time step used in the simulation Δt = 10 −3 , was shorter than the critical time step Δt critical = 10 −2 , computed using the Gershgorin circle theorem, for Re = 200. The initial values for all the flow variables were set to zero.
To demonstrate the accuracy of the flow solver, we compare our steady state results with those reported by Shu et al. 61 The variation of velocity component u along the z-coordinate, and the variation of velocity component v along the x-coordinate, are used to examine the accuracy of the proposed scheme. Tables 3 and 4 show the numerical results for the u-and v-velocity profiles along the vertical and horizontal lines, respectively, through the center of the cavity (x = 0.5, y = 0.5, z = 0.5) for Re = 100 and Re = 200. The numerical results obtained using the proposed scheme are in very good agreement with the results reported in Shu et al. 61 We show the velocity profiles of u component along the vertical centerline and v component along the horizontal centerline on the plane z = 0.5 for Re = 100 and 200 ( Figure 14). It can be found that all the velocity profiles generated using the proposed method agree very well with the results obtained using the established approaches reported in Shu et al. 61 Figure 15A displays the u velocity component profile along the vertical line along the z-axis through the center of the cavity for Re = 100. The results are in very good agreement with the numerical results computed using ANSYS CFX, highlighting the accuracy of the proposed IB method. Software implementation that maximizes the computation speed of the proposed method is beyond the scope of this study. In our simulations, we used in-house MATLAB version R2019a code executed on Intel i7 quad core processor with 16 GB RAM. For this code, the computation speed of our method was comparable with that of ANSYS CFX steady state solver with the root mean square error (RMSE) set to 10 −8 . It may be anticipated that the computation speed of our method would appreciably increase if software implementation is done using a compiled programming language such as C, C++, or Fortran. Figure 15B displays the streamlines in the presence of the sphere for Reynolds number Re = 100. The RMSE between the numerical results computed using the proposed scheme and ANSYS CFX are less than 1% for all three velocity components.
We examine the accuracy of the proposed scheme using a locally refined Eulerian grid. We discretize the flow domain using 548,485 nodes, locally refined in the vicinity of the immersed sphere ( Figure 16). The nodal distribution is dense to ensure a grid independent numerical solution. To obtain the Eulerian point cloud, we first generate a Cartesian grid of 51 × 51 × 51 = 132,651 nodes, with uniform resolution of h = 0.02 m. Second, we locally refine the Cartesian grid in two rectangular patches (0.2, 0.2, 0.2) ≤ (x, y, z) ≤ (0.8, 0.8, 0.8) with resolution h 1 = h 2 = 0.01. Finally, we locally refine the Eulerian grid points that form a thin layer around the IB ( Figure 16B). The finest grid resolution of the Eulerian nodes is h 3 = h 8 = 0.0025. The numerical results obtained using the locally refined Eulerian grid are in excellent agreement with those computed using the uniform Cartesian grid.

NUMERICAL EXAMPLES
In this section, we highlight the accuracy of the proposed method on benchmark flow problems in two and three dimensions.

Two cylinders moving toward each other
In this example, we examine the accuracy and the applicability of the proposed method on a case that involves flow past two moving solid objects. This flow example was first investigated by Russell and Wang,62 and later on by Xu and Wang. 63 The geometry (initial location of the two cylinders) and the boundary conditions (the far-field boundaries are rigid walls) are shown in Figure 17A. To eliminate the effect of the impulsive start, both cylinders oscillate about their initial position for two time periods, and then start moving toward each other (with velocity u s ) with a speed that results in Re = 40.
We conduct a mesh convergence study to obtain a grid independent numerical solution. We use a Cartesian grid with uniform resolution h, which is then twice locally refined ( Figure 17B) with half resolution, initially at −5 ≤ x ≤ 20 and −5 ≤ y ≤ 5, and following at −3 ≤ x ≤ 18 and −3 ≤ y ≤ 3. The coarse grid consists of 10,659 points. The initial resolution is h = 0.5 and after two refinements h r = 0.125. The moderate mesh has 43,509 points (the initial resolution is h = 0.25 and after two refinements h r = 0.0625), and the denser one has 175,785 points (the initial resolution is h = 0.125 and after two refinements h r = 0.03125). It was found that solutions on the final two nodal configurations are grid converged. Therefore, in our simulations we use 175,785 points to represent the flow domain. The critical time step (computed using Gershgorin theorem) for this nodal distribution and for Re = 40 is Δt critical = 5 × 10 −2 . Both cylinders move with predefined motion. The updated position of the lower and upper cylinder are described by the trajectories of their centers (x l , y l ) and (x u , y u ) through y l = 0 and x l = , and y u = 1.5 and x u = Both cylinders oscillate about their initial positions, and then start to move toward each other at time t = 16. Figure 18 displays the time evolution of drag and lift coefficients for the upper cylinder. The numerical results obtained by the proposed scheme are in good agreement with those reported by Xu and Wang. 63 The drag slightly increases when the two cylinders approach each other and decreases when they are close to each other. Then, drag rises back and begins to approximate the profile of a single cylinder, as pointed out in Reference 63. On the other hand, the two cylinders are repulsive when approaching each other and then become attractive after they have passed each other. Figure 19 displays the vorticity and stream function contours at several time instances.

Flow in rectangular duct with a sphere and non-convex object
As a final example, we consider the flow in a square duct with a sphere located downstream close to the inlet of the duct. The duct has dimensions 0 ≤ x ≤ 2.5 m and 0 ≤ y, z ≤ 0.4 m, the sphere is located downstream at O(0.5, 0.2, 0.2) and has radius R s = 0.1 m. We discretize the flow domain using 727,288 nodes, locally refined in the vicinity of the immersed sphere ( Figure 20). The nodal distribution is dense to ensure a grid independent numerical solution. To obtain the Eulerian conditions are applied, while a parabolic stream-wise velocity profile is set at the inlet At the outlet, we apply a fully developed flow condition . The boundary conditions for the vector Boundary conditions for the vorticity field can be computed using = × u| Ω . Following Wong and Reizes, 58  vector u inlet . The velocity field u may now be written as The fluid kinematic viscosity was set to v f = 10 −3 m 2 /s, yielding a fixed Reynolds number Re = 100, based on the sphere diameter D = 0.2 m and the mean flow velocity U mean = 4 9 U max at the inlet. Figure 22 shows the streamlines in the flow domain at various time instances (at t = 1 s and t = 3 s). In the absence of an exact analytical solution, the results obtained with the proposed IB method are compared against those computed using the incremental pressure correction scheme (IPCS) FE solver Oasis 64 at different time instances (from t = 0.1 to t = 3 with time interval of Δt = 0.1 s). The body-fitted mesh consists of 308,670 nodes and 1,774,875 linear triangular elements (P 1 ∕P 1 ). We computed the normalized root mean square error (NRMSE) for all three velocity components and found that the highest error over all time instances is NRMSE = 1.2 × 10 −2 (1.2%), highlighting the accuracy of our MM. To demonstrate the applicability of the proposed method to 3D flow cases that involve non-convex immersed objects, we replaced a sphere with Möbius torus when modeling a flow in rectangular duct ( Figure 23A). The flow domain and the boundary conditions are the same as in the flow behind the sphere example ( Figure 20). The spatial domain is represented by 642,509 nodes, locally refined in the vicinity of the IB ( Figure 23B). The immersed object (Möbius torus) is discretized using 9980 nodes. The computed streamlines in the flow domain at various time instances (at t = 1 s and t = 3 s) are shown in Figure 23C,D.

CONCLUSIONS
In this work, we developed a new meshless IB method by combining the BCE-IB and the DC PSE methods to numerically solve the nonstationary, incompressible N-S equations in their stream function-vorticity formulation. We have extended the proposed method to three dimensions using the vector potential-vorticity formulation. We computed the spatial derivatives of the unknown field functions (velocity and vector-potential) using the DC PSE method, which is an accurate interpolating meshless method MM. The method works efficiently and accurately for uniform Cartesian grids, and for Cartesian grids locally refined in the vicinity of the immersed objects. We used the Euler explicit time integration scheme and computed the critical time step that ensures solution stability using Gershgorin's circle theorem. The proposed scheme using DC PSE for computing spatial derivatives results in a relatively large time step. We demonstrated the accuracy of the proposed scheme in 2D (stream function-vorticity formulation) and 3D (vector potential-vorticity formulation) using three benchmark problems: unbounded flow past a cylinder for Re = 40 (where an analytical solution is available) and Re = 550 (where transient phenomena start to appear); in-line oscillating cylinder in fluid at rest (where experimental data are available); and 3D lid-driven flow in a cubical domain with a fixed sphere located in the middle of the cube. We highlighted the applicability of the proposed scheme in more computationally demanding flow cases such as the motion of two cylinders moving toward each other, and 3D flow in a rectangular duct with spherical and non-convex (Möbius torus) obstacles. Examples of application of the proposed scheme to 2D flow problems that involve non-convex objects can be found in Reference 65.
The proposed method appears as both computationally efficient, since large time steps are used, and accurate. To the authors knowledge this is the first time meshless point collocation method is combined with an IB method to solve the N-S equations in their vector potential-vorticity formulation.