In fluid simulations involving free-surface flows, Lagrangian finite element methods represent an effective numerical tool since fluid boundaries are tracked automatically by the position of the mesh nodes, allowing for a convenient treatment of the fluid domain characterized by fast evolutions. Among the different numerical methods recently proposed in the literature, the Particle Finite Element Method (PFEM) is a mesh-based Lagrangian approach for fluid modelling, particularly suited for problems with rapid changes in the domain topology. Starting from a set of fluid particles, the PFEM exploits the Delaunay tessellation and the alpha-shape method to create a computational mesh whenever the current one becomes overly distorted due to fluid motion [1]. An explicit time integration scheme is generally preferred in fast dynamics analyses. However, it is only conditionally stable and it requires small time-step sizes, which can be estimated based on the smallest element dimension in the computational mesh. Consequently, overly distorted elements can drastically reduce the stable time increment size. Moreover, the efficiency of the runtime remeshing plays a key role in Lagrangian explicit analyses. In 2D, the Delaunay tessellation guarantees optimal geometrical properties of the created triangles, such as the minimization of the maximum radius of an element circumcircle and the maximization of the minimum angle among all the elements. In contrast, the 3D Delaunay tessellation is less robust and it allows for the generation of badly shaped tetrahedra (slivers). These elements have dramatic impact in explicit approaches, leading to very small time-step sizes. Several mesh optimization algorithms have been introduced in the literature for the treatment of slivers (e.g. see [2]). Even though these techniques have been successfully applied in several engineering problems, they still have limited applicability due to the high computational cost when runtime remeshing is required or in the case of fixed nodes on the boundaries that cannot be relocated in the mesh smoothing operation. In the present work, the introduction of Virtual Elements (VEs) into the PFEM is proposed. The idea is to agglomerate [3] two badly shaped triangular (2D) or tetrahedral (3D) elements giving rise to a bigger VE element with 4 edges (2D) or six faces (3D) exploiting the possibility offered by VEs to have arbitrary shapes and number of edges (2D) or faces (3D).