The Lagrangian meshfree particle-based method has advantages in solving fluid dynamics problems with complex or time-evolving boundaries for a single phase or multiple phases. A pure Lagrangian meshfree particle method based on a generalized finite difference (GFD) scheme is proposed to simulate time-dependent weakly compressible viscous flow. The flow is described with Lagrangian particles, and the partial differential terms in the Navier-Stokes equations are represented as the solution of a symmetric system of linear equations through a GFD scheme. In solving the particle-based symmetric equations, the numerical method only needs the kernel function itself instead of using its gradient, i.e., the approach is a kernel gradient free (KGF) method, which avoids using artificial parameters in solving for the viscous term and reduces the limitations of using the kernel function. Moreover, the order of Taylor series expansion can be easily improved in the meshless algorithm. In this paper, the particle method is validated with several test cases, and the convergence, accuracy, and different kernel functions are evaluated.2 of 21 rising bubbles [3] and coalescing [4]. The ability of this method to model non-Newtonian fluid and large scale diffuse fluids has been demonstrated in some recent works [5,6] by introducing different symmetric models. Moreover, because the computations are based on the support domain, which is much smaller than the complete computational region, the ill conditioned system problem is rarely encountered. Among all Lagrangian meshfree methods, the smoothed particle hydrodynamics (SPH) method was one of the earliest methods developed and has been widely applied in different fields. The SPH method was first pioneered independently by Lucy [7] and Gingold and Monaghan [8] to solve astrophysical problems in 1977. Details of the SPH method as a computational fluid dynamics method can be found in recent reviews [9][10][11][12] and Liu and Liu's book [13]. Some successful applications of this method include coastal engineering, nuclear engineering, ocean engineering, and bioengineering. However, the accuracy of the conventional SPH method is unsatisfactory, and it is not easy to achieve an accurate high-order SPH approach.As a meshfree Lagrangian method, the particle distribution generally tends to be irregular in the computations, which leads to inconsistency and low accuracy [14,15]. For that reason, in some cases, only the first-order term of the fluid dynamics equations, the Navier-Stokes equations, is solved, and the viscous term, which contains the second-order differential, is obtained through the artificial viscosity with artificial parameters in the SPH method. This issue can also occur in the incompressible SPH (ISPH) method [16]. To improve the consistency and accuracy of these methods, different modifications have been developed. After using Taylor series expansion to normalize the kernel function, the corrective smoothed particle method (CSPM) [17,18] and the modified smoothed particle meth...