A Review of Vortex Methods and Their Applications: From Creation to Recent Advances

: This review paper presents an overview of Vortex Methods for ﬂow simulation and their different sub-approaches, from their creation to the present. Particle methods distinguish themselves by their intuitive and natural description of the ﬂuid ﬂow as well as their low numerical dissipation and their stability. Vortex methods belong to Lagrangian approaches and allow us to solve the incompressible Navier-Stokes equations in their velocity-vorticity formulation. In the last three decades, the wide range of research works performed on these methods allowed us to highlight their robustness and accuracy while providing efﬁcient computational algorithms and a solid mathematical framework. On the other hand, many efforts have been devoted to overcoming their main intrinsic difﬁculties, mostly relying on the treatment of the boundary conditions and the distortion of particle distribution. The present review aims to describe the Vortex methods by following their chronological evolution and provides for each step of their development the mathematical framework, the strengths and limits as well as references to applications and numerical simulations. The paper ends with a presentation of some challenging and very recent works based on Vortex methods and successfully applied to problems such as hydrodynamics


Introduction
Vortex methods (VM) belong to Lagrangian methods, also called particle methods, used to solve continuous systems of the form: where Q denotes the unsteady solution vector evolving in space (x) and time (t), composed of the physical quantities Q i , and where F is the source term depending on space and time and on the solution Q and its partial derivatives. In a general point of view, Lagrangian methods differ from Eulerian ones by the way they discretize the problem. In an Eulerian approach, the computational domain is described as a mesh and the solution vector Q is approximated on each node according to a numerical scheme involving the solution at the neighboring mesh nodes. In the Lagrangian framework, the domain discretization is not realized on a mesh but on a set of particles that follow the dynamic of the system and are displaced with respect to the flow velocity u. More precisely, in the Lagrangian approach, the solution Q is discretized on numerical particles and, given a location x in the domain, each quantity Q among the vector solution Q is thus approximated at this location x by: where δ denotes the Dirac distribution. To write this expression in a regularized and discrete form, the Lagrangian methods define a volume v p for the particles as well as a smoothing function W, of compact support and radial symmetry of radius ε, which tends to the Dirac distribution when ε tends to 0. Therefore, the discrete approximation of the quantity Q at location x can be written as: In other words, the quantity Q at location x is evaluated as the sum of the quantity Q carried by the particles p located in the compact support of the smoothing radial function W of radius ε, centered at location x. The Lagrangian discretization of the continuous system (1) is therefore obtained by solving the following system of ordinary differential equations (ODEs) in order to update the particle positions x p , their volume v p and each quantity Q p they carry according to the source term F: The key common thread of all particle methods used for continuum flow simulations therefore relies on the approximation in the Lagrangian form of the governing continuous Equation (1) where the derivative operators involved in F are replaced by equivalent integral operators that are then discretized on the particle locations. Among the most popular particle methods used to solve continuous equations, one can cite the Vortex Methods (VM) and the Smooth Particle Hydrodynamics (SPH). They differ from each other by the governing Equation (1) they aim to solve and by the way they handle the resolution of the third ODE of system (4). SPH methods [1][2][3][4][5] solve the weakly compressible Navier-Stokes equations in velocity-pressure formulation and the particles carry the pressure and the flow density variables. They need a coupling with an equation of state to close the governing system. Regarding Vortex Methods, they are used to solve the incompressible Navier-Stokes equations in their velocity-vorticity formulation and the particles only carry the vorticity field.
If one starts from a present observation, one admits that in the computational fluid dynamic community, Vortex Methods do not benefit from the same recognition or the same expanded scope than Eulerian methods (spectral methods, finite elements, finite volumes, finite differences, etc.). However, they early raised the interest of the CFD research community due to their capability to precisely and naturally simulate flows dominated by the advection phenomenon. They have been the subject of plentiful of research works, which allowed these last recent years to lead to challenging applications involving important issues in fluid dynamics, in particular in hydrodynamics [6][7][8][9], in flow mediated transport [10,11], in turbulent wake dynamics of rotorcrafts [12,13], aircrafts [13,14] and vertical axis wind turbines [15,16], in aeroacoustics [17], in free jet flows [18] or in active [19] and passive [20] flow control. One of the first thorough reviews of Vortex methods has been proposed in 1991 by Sethian and Gustafson [21]. In 2000, Cottet and Koumoutsakos [22] proposed a further detailed book, gathering all the theory and numerical practice related to these methods. In 2011, Yokota and Obi completed this overview with a review paper dedicated to the simulation of turbulent flows with Vortex methods [23].
The present paper, starting from preliminary survey carried out in [24], proposes a synthetic review of Vortex methods from their genesis in the 1930s until their very recent advances. In the context of incompressible flow simulations, one of the first goals of this work is to provide an insight of a particle approach based on a solid mathematical framework with convergence analysis and conservation features, which is able to provide accurate simulations and which could be considered as a complement or an alternative to purely Eulerian methods. Concerning the Vortex methods themselves, the literature survey of this review aims to give an overall landscape (although not exhaustive due to the finite format of a review paper) of the existing Vortex Methods, depicting the different sub-families of VMs with their own advantages, drawbacks and applications. It also aims to highlight the intrinsic advantages of Vortex Methods that explain their intensive and continuous development from the 1970s, and attempts at showing how they progressively evolved in time to meet the requirements involved by the CFD community, namely the capability of handling complex physics while ensuring convergence to solution, accuracy, efficiency, and performance. This review is addressed both to scientists specialized in numerical analysis and discretization schemes, and those dedicated to incompressible fluid dynamics using as well direct numerical simulations or large eddy simulations.
The present review paper follows the chronological advances of Vortex Methods and is structured this way: First, Section 2 gives the global description and formulation of Vortex Methods in the initial context they were used for, that is to say the resolution of the inviscid Euler equations. In Section 3, one will explain how they expand to solve the Navier-Stokes equations, including the diffusion effects. Section 4 will address the crucial issue of particle distortion in Lagrangian methods and their need to perform a spatial adaptation of the particle locations in order to ensure the numerical convergence; this section introduces the so-called hybrid Vortex Methods, which include in their discretization schemes the presence of an underlying grid in order to circumvent the particle distortion phenomenon. In Section 5, we will focus on the treatment of boundary conditions in VMs, and more precisely on the enforcement of no-slip boundary conditions at the surface of a solid and rigid immersed obstacle. Section 6 will describe the progresses made in VMs in terms of algorithmic issues, allowing them to handle complex or computationally costly problems. Finally, Section 7 will propose a finite list of relevant and challenging simulations performed with Vortex Methods.

General Description and Formulation
The first calculations based on vortex methods were performed, by hand, independently by Prager in 1928 [25] and Rosenhead in 1931 [26]. For that purpose, they discretized the flow into elements, namely the particles, characterized by their position and the local circulation they carry. One recall that the circulation Γ corresponds to the macroscopic rotation of the flow. It is a scalar field defined as the integral of the vorticity field ω on patch of vorticity A, bounded by a closed path: where the vorticity ω measures the local rotation of the fluid flow and is obtained by taking the curl of the velocity field: In the case of a two dimensional flow, the vorticity is perpendicular to the flow plane and corresponds to a scalar field. We note it by ω in that case. In the Cartesian coordinates system, the vorticity respectively expresses in 2D and 3D as follows: Based on the pioneer works [25,26], the method described in the sequel corresponds to the general definition of vortex methods, first introduced in the 1970s for the twodimensional formulation [27] and in the late 1980s for the three-dimensional formulation [28,29].
The design of vortex methods for inviscid incompressible flows is based on vorticityvelocity formulation of the Euler equations. The latter represent hyperbolic equations of conservation of mass (continuity), and balance of momentum and energy. They correspond to particular Navier-Stokes equations for inviscid (zero viscosity) flows, without thermal conductivity. In vorticity-velocity formulation and in the case of incompressible flows, they can be written as follows: ∂ω ∂t + (u · ∇)ω − (ω · ∇)u = 0 (8) div u = 0 (9) ω(·, 0) = ω 0 (10) |u| → u ∞ (11) where (u · ∇)ω denotes the advection term, (ω · ∇)u = [∇u][ω] corresponds to the stretching term (which vanishes in two-dimensions) and where Equation (9) corresponds to the incompressibility condition, resulting from the continuity equation. In fact, the Euler equations were originally formulated in the Lagrangian form which involves Dω/dt, the material derivative in time of the vorticity ω with respect to the advective field u. The discretization of these equations with a vortex method, that is to say with a particle method, is therefore natural and straightforward. It goes as follows: Let a regular approximation ω h of ω, and let x h p the trajectories of the particles transported by the flow, the approximated vorticity ω h may be numerically expressed at time t by: where the numerical particles locations x h p (t) are obtained by solving the following ODEs: and where α h p corresponds to the local circulation around x h p , defined as the product of the vorticity by the volume of the particle α h p = v p ω h p . According to the Kelvin's theorem, which states that in an inviscid flow the circulation around a closed curve remains constant in time when moving with the fluid, the local circulation α h p satisfies the following ODEs in 2D and 3D respectively: From the Kelvin's and Helmholtz theorems, the vorticity is advected at the local fluid velocity u h in the potential flow (one refers the reader to section 1.3 of Cottet and Koumoutsakos' book [22] for a detailed presentation of the first and second Helmholtz theorems that drive the motions of fluid particles carrying vorticity in inviscid flows and their further extension by Kelvin to include the effects of viscosity). The velocity can be recovered at all times from the vorticity through the following Poisson equation, based on the incompressibility condition: which can be solved through an integral representation: where G is the Green's function for the Laplacian operator, ∇ 2 G(x, x ) = δ(x − x ) and where denotes the convolution product. Then, from Equation (16) and denoting K = ∇G, one finally obtains the velocity field through the following Biot-Savart equation: with in 2D : (18) in 3D : To summarize, in vortex method, each vortex element is made to move at a single velocity according to the Biot-Savart law and Equations (12)- (14) and (17) model the transport of the vorticity in inviscid incompressible flow. However, the Dirac distribution involved in the vorticity definition (12) makes unaffordable the numerical discretization and approximation of the vorticity field in practice. Faced with this problem arose the idea of the vortex blob method, originally introduced by Chorin and Bernard [27]. In this numerical approach, the vorticity of each particle is distributed on a blob, that is to say on a disk of finite radius ε. This vorticity blob-distribution is defined by a smooth cutoff function ζ.
If the blob radius ε is small enough, one denotes ζ ε the mollified function defined in 2D and 3D by: Among the possible cutoff functions, the Gaussian-based ones are often favored because of their smoothness property, their radial symmetry and their fast decay. In vortex blob method, the vorticity ω h ε of the mollified particles is then given by the following expression (one refers to Figure 1 for a graphical support): with dα h Concerning the resolution of three dimensional inviscid flows using vortex method, most works mainly rely on the 3D formulation exposed here, which is a direct extension of the 2D method, but it has to be noticed that the very first proposition of 3D vortex method was based on a different approach, called the vortex filament method. The latter was introduced by Leonard [30][31][32] and leans on the theorems of Helmholtz and Kelvin, stating that tubes of vorticity retain their identity and move as material elements in an inviscid flow. In this approach, the quantity carried by the Lagrangian elements is a scalar, as in two dimensions, and the elements are now deformable lines (or filaments) instead of discrete points. Also, as for two-dimensional vortex schemes, the circulation of the particles remains constant in time. In practice, vortex filament methods were used to perform one of the first 3D simulations in CFD by Leonard [31] and were also used more recently by Angelidis and Neyret [33] with smoke simulations in graphic works. However, despite their very natural approach from a physical point of view, they suffer from several weaknesses: first of all, they are most of the time quite difficult to initialize for a given flow. Only specific problems like vorticity rings and jets are spared from this trouble. Also, some filament "surgeries" might become necessary in order to avoid tiny loops or to reconnect neighboring filaments, which make the method complex to use. Eventually, there is no straightforward way to include the physical diffusion effects in the model with filament method. In Vortex blob methods, the vorticity ω h at location x is approximated as the sum of the local circulations α h p carried by the neighboring particles of volume v p and located in the support of the ζ ε mollified function.

Convergence Issues
What distinguishes vortex methods from grid-based methods is the fact that they are built in such a way that many inviscid flow invariants are conserved [22]. This makes them particularly interesting since they ensure a correct qualitative answer, even in the case of underresolved simulations. Besides this important feature, vortex methods indisputably belong to the category of plain numerical methods by means of their complete numerical analysis and their theoretical proofs of convergence.
The works of Hald in 1979 [34] and Beale and Majda in 1982 [35] provide the complete proof of convergence as part of the vortex blob method in two dimensions. Their analysis included the influence of the numerical parameters involved in the model, namely the blob radius ε (that is to say the cutoff width) and the particle spacing, denoted h. It turned out that the vortex blob trajectories tend to the exact one when the number of eddies increases, more precisely when h = O(ε). This result allows us to reach high orders of convergence in two dimensions. An extension of this convergence result to the three dimension case is also given in [35] as well as in the work of Cottet [29], giving birth to new high orders three dimensional vortex methods. The convergence condition h = O(ε) is also clearly highlighted numerically by the simulations performed by Knio and Ghoniem [36] concerning the propagation and stability of three-dimensional vortex rings. The temporal discretization is added for the first time in the convergence analysis of the 2D and 3D Euler equations by Anderson and Greengard in 1985 [28]. The most important result which came out of this study is the fact that these methods do not suffer from any Courant-Friedricks-Levy (CFL) condition, which constraints the time step to be related to the minimal grid size in Eulerian methods. Indeed, in vortex methods the stability constraint on the time step only depends on the velocity gradient: where the constant C, also called LCFL number (for Lagrangian CFL), must be lower or equal to 1. This condition provides a much less restrictive constraint for the time step compared to purely grid-based schemes. After this general presentation, one can distinguish several advantages using the discrete vortex method. First, its Lagrangian approach with particle discretization makes it intrinsically close to the flow physics it aims to model and provides a direct scheme to solve the transport of particles and the quantities they carry, based on ODEs systems, without numerical dissipation. Secondly, the vorticity field is a compact field compared to the velocity. It is therefore much more efficient to directly solve and simulate the vorticity field. Then, if the fluid density is constant, the Euler equations for incompressible flows does not contain pressure term in vorticity formulation. Furthermore, in the discrete vortex method the boundary condition at infinity can be automatically satisfied contrary to most grid-based approaches where external flows must be solved within a finite simulation domain. Finally, the stability condition (25) provided by the vortex method scheme allows larger time steps compared to grid-based methods.

Viscous Vortex Methods
The first numerical schemes including viscous effects in vortex methods were designed with the desire to conserve the Lagrangian framework. They were consequently purely based on particles. Before starting with the description of the very first viscous vortex method, namely the Random Walk method developed by Chorin in 1973 [37], it is necessary to introduce the concept of viscous splitting of the Navier-Stokes equations. For sake of clarity, in the rest of this paper the approximated quantities ω h ε , α h ε and u h ε will be simply referred to as ω, α and u.

The Viscous Splitting and Random Walk Method
The original concept of viscous splitting was first introduced by Prandtl, in 1904 [38]. However, its design within the framework of vortex methods was originally proposed in the early 1970s by Chorin [37] in order to develop his Random Walk Method suited to model the diffusive phenomenon. The viscous diffusion equation is: where ν denotes the viscosity of the fluid. This equation is coupled with the Biot-Savart Equation (24) in order to recover the velocity field from the vorticity. The viscous splitting algorithm consists in performing substeps in which the advective and the diffusive effects are handled separately and successively. Thus, in the vortex method frame, the two steps of the viscous splitting algorithm for a two-dimensional flow are the following: (1) advection: (2) diffusion: where in the first step, the particles locations (x p ) are updated using the local flow velocity obtained with the Biot-Savart law and in the second one the vorticity (ω p ) is modified at these new positions by the diffusion. The convergence of the viscous splitting method was proved by Beale and Majda in 1981 [39].

Random Walk Method (RW)
The method of Random Walk, created by Chorin [37], consists in solving Equation (30) through the probabilistic interpretation of the Green's function, solution of the diffusion equation: which, in d-dimensions is known to be explicitly equal to: As mentioned previously, the RW method is inherently related to the viscous splitting approach since the vorticity field obtained after the particle push (Equations (27) and (28)) is assigned as the initial condition of the diffusion Equation (30). The diffusion effects are then modeled imposing a Brownian motion to the flow particles whose positions are updated over a time step ∆t according to the following scheme: where g is a vector of independent random variables obtained through a Gaussian probability distribution with zero mean and a variance σ = 2ν∆t. Moving the particles in such a way creates a concentration field that converges to the above Green's function solution of the diffusion equation, as the number of particles tends to infinity. Marchioro and Pulvirenti [40], Goodman [41] and Long [42] successively proved the convergence of the Random Walk method. Thank to its simplicity, its property to conserve the total circulation and also its capability to easily handle flows around complicated boundaries (see Section 5.1), this method has been extensively used. In terms of applications, we can cite the works of Chorin [37,43], Dommelen [44] and Stansby [45,46] who applied the Random Walk method for the simulations of flows around cylinders. We can also refer to Wang [47] who studied in his thesis several active flow control techniques to prevent the dynamic stall of airfoils and to Lewis [48] who used the Random Walk method to simulate flows over airfoil cascades. However, this viscous vortex method is a stochastic method: it therefore presents several drawbacks. First, it does not exactly conserve the mean position of the vorticity in free space. Next, the computed solutions are noisy due to the statistical errors. Indeed, increasing the viscosity ν makes the variance of the random variables larger and consequently deteriorates the solution. Finally, the stochastic feature of the Random Walk method limits its convergence capability. The error is indeed proportional to 1/ √ N as explained by Milinazzo and Saffman [49], where N is the number of particles. Therefore, accurate simulations need a very large number of fluid elements, which increases the computational cost and explains their rare application to three-dimensional simulations. One can however cite the use of RW method in the simulation of an impulsively started flow over a 3D cube at Re L = 100 performed by Gharakhani and Ghoniem [50] (with L the cube length) As alternatives to RW, multiple deterministic schemes abounded. The most emblematic ones are described in the following.

Particle Strength Exchange Method (PSE)
Among all the viscous deterministic vortex methods, the PSE, introduced by Degond and Mas-Gallic in 1989 [51] from the previous connected works of Cottet and Mas-Gallic [52,53], is considered as one of the most popular. The mathematical description of the PSE method is here given in a two-dimensional framework but the main ideas of this approach remain the same in three dimensions:

1.
Approximate the diffusion operator ∆ω by an integral operator: where η ε is a mollified diffusion kernel of order r, satisfying precise moment properties [22].
The particle approximation of the diffusion equation can be obtained from the numerical integration of (35). The resulting vortex scheme is given by: where the particles trajectories x h p are computed from (27) (based on Biot-Savart law), and ω h p (t) satisfies the following ODEs: where the right hand side is obtained from (35) through a discretization of the integral operator taking the particle positions as quadrature points.
Equation (37) can be interpreted as changes in the strength of particles due to the neighboring particles, which explains the name of "Particle Strength Exchange". Beyond its deterministic framework and its integral formulation, the accuracy and reliability of the PSE method is also based on its capability to conserve the natural invariants of the flow. Moreover, the method is not based on the viscous splitting algorithm, which therefore discards the intrinsic error of the time-steping process. The method exhibits an error bound of O(ε 2 ) for simple Gaussian smoothing. However, the quadrature rule used to discretize the integral makes the PSE method strongly dependent on the particles distribution, which has to respect the convergence condition h = O(ε) as well as possible. To circumvent this drawback, particle spatial adaptions often need to be added to the PSE method. This point will be addressed in details in Section 4.
A wide range of applications have been handled so far based on this deterministic approach which may be considered as one of the most efficient and accurate viscous Lagrangian vortex methods. Among them, one must cite the pioneer and reference works of Koumoutsakos and Leonard [54] who applied the PSE method to 2D viscous flows around translating and rotating cylinders for Reynolds numbers ranging from 40 to 9500. One can also cite the works of Ploumhans et al. [55] for the simulation of transient 3D flow past a sphere or those of Yokota et al. [56] for the study of 3D isotropic turbulence in a periodic domain. All these numerical studies highlight the necessity of particle spatial adaption in order to ensure the convergence of the vortex PSE method.

Core-Spreading Method (CSM)
This method is not based on integral operators approaching the Laplacian (like in PSE) but directly considers the exact solution of the diffusion Equation (26), given by the Green's function solution: In order to discretize the diffusion equation with a vortex scheme, if one considers the numerical scheme provided by the vortex blob method described in the first section of this paper, one has: where the smooth cutoff function ζ ε distributes the vorticity on a blob of radius ε.
The core-spreading method, first introduced by Kuwahara and Takami in 1973 [57] and later reformulated in 1980 by Leonard [31] has its origin in the observation that, if one chooses ζ ε as the Gaussian distribution: then scheme (39) reduces to the solution of the heat Equation (38), provided the variance of the Gaussian blob distribution expands in time according to ε 2 = 2νt, which explains the denomination of "core-expansion" or "core-spreading" for the method. However, Greengard [58] highlighted that a the straightforward implementation of CSM lacked of consistency and could only give a converged approximation of the Navier-Stokes equations under the very restrictive limit of vanishing viscosity. Indeed, in a general point of view, the expansion in time of the vortex blobs radius is not in agreement with the convergence criterion stating that ε must be a small parameter. To circumvent this inconsistency, a spatial adaptation procedure was proposed by Rossi [59] in order to keep control on distance between particles. Based on this improvement (and provided a careful treatment of the local spatial adaptation further described in Section 4.1), the work of Barba et al. [60] demonstrated the means for using CSM in the case of a Lamb vortex dynamics study, and Yokota et al. [56] also used it successfully for their study of 3D isotropic turbulence. This approach conserves the deterministic treatment of the viscous diffusion equation and preserves the meshless framework of vortex methods.

Diffusion Velocity Method (DVM)
This method, first introduced by Fronteau and Combis in 1984 [61] and later used in the vortex methods framework by Ogami and Akamatsu in 1991 [62], relies on the idea that the diffusion of vorticity is handled within the advection process by absorbing the diffusion term into the advection one in the Navier-Stokes equations ∂ω ∂t + u · ∇ω = ν∆ω.
This process gives rise to two different velocities: the advection velocity u and the diffusion velocity u d . The latter can be considered as an artificial velocity and is defined, in twodimensions, as follows: From this definition, the diffusion problem is now transformed into an advection equation: Thus, the transport of vorticity (due to advection and diffusion effects) is modeled by the vortices' movement provoked by the advection velocities and the diffusion velocities. The DVM has been compared to PSE approach in the work of Beaudoin et al. [63] for the 2D simulations of pollutant transport in groundwater, featuring that for such dispersion problem the PSE tends to provide more accurate results while DVM seems to be more stable. However, the DVM requires an important particle overlapping, as precisely highlighted in [63] and in the earlier work of Clarke and Tutty [64] simulating the impulsively started two-dimensional flow past an arbitrary shaped body using DVM associated with a random walk (RW) procedure. Concerning the comparison of diverse vortex diffusion methods one finally refers to a very recent 2D study brought by Lucchesi et al. [65] confronting a direct differentiation approach (non-conservative), with three different PSE schemes and a DVM method (all conservative by construction). An extension to the three-dimensional formulation of the DVM was proposed by Rivoalen et al. [66], where a diffusion vorticity ω d (also called "viscous vorticity") has to be added to the problem, giving the new vorticity transport equation: which, in the Lagrangian framework, is solved according to the following system of ODEs: However, it is important to notice that ω d needs the approximation of the Laplacian operator, which in practice in [66] is performed using a discretization of the integral approximation (as done in PSE method).
These three deterministic viscous methods, namely the PSE, the CSM and the DVM, represent the grid-free viscous approaches that allowed conservative and converged numerical simulations for two and three-dimensional problems. It is interesting however to mention several other methods that were proposed to solve the viscous diffusion equation with vortex schemes, despite their failure in successfully applying it to the three-dimensional case. One can for instance cite the one introduced by Fishelov in 1990 [67], called the Convolution of Cut-off Function method (CCF). It allows us to obtain the Laplacian operator by using explicitly second derivatives of the cutoff function, i.e., ∆ω ≈ ω ∆ζ ε , in the vortex blob scheme: But this method encounters several weaknesses. It has been first observed to be more sensitive than PSE to particle field distortion, but overall, the formulation (46) lacks the correction term −ω p ∑ q v p ∆ζ ε (x h q − x h p ) and is therefore not conservative [22]. These reasons explain why this approach failed in giving rise to deeper research works.
Among the other vortex schemes proposed to solve the diffusion equation, one can also mention the Vortex Redistribution Method (VRM) of Shankar and van Dommelen (1996) [68]. This method was thought as an alternative to PSE with no need to perform particle rezoning. As PSE, the modeling of diffusion is built on the exchange of circulation between neighbor particles, but in VRM the scheme formulation in not based on integral operators. The exchange of circulation between particles in modeled through the computation of functions that will determine the amount of circulation each particle has to redistribute to its neighbors, located in the area limited by the typical diffusion distance. The definition of the redistribution fractions is designed such that local circulation as well as linear and angular momentum are conserved. However, their evaluation implies many parameters whose computation might be very costly (even more that the cost of velocity computation through Biot-Savart law [22]), which in practice makes this approach difficult to use.
For an exhaustive description about grid-free viscous vortex methods, one can refer the reader to the detailed literature survey given in Barba's thesis [69] and to Cottet and Koumoutsakos [22].
The important number of diverse viscous vortex schemes illustrates the difficulty of modeling viscous incompressible flows with a grid free vortex method. As emphasized in the present section, the main limitation of all the proposed schemes does not come from their capability of modeling the viscous effects themselves, but from the necessity to ensure particle overlapping throughout the whole calculation. Studies on spatial adaptation procedures were therefore developed in parallel in order to address this key point. We describe them in the following section.

Spatial Adaptation and Hybrid Vortex Methods
In pure Lagrangian method, the particle distribution may undergo distortion, that manifests itself by the clustering or spreading of the flow elements in high strain regions. This accumulation or dispersion phenomenon is then translated in inaccurate computations because of a deteriorated communication between particles. As meaningful illustrations of the distortion effect one can cite the test case of the axisymmetric inviscid vortex patch, illustrated for instance by Rossinelli and Koumoutsakos [70]. In this example, one can observe the distortion of the initial and uniform particle distribution due to the important shear. One can also see in Henshaw et al. and in Krasny's studies [71,72] comparable behaviors concerning the appearance of non-physical small scales in vortex sheet calculations in the context of two dimensional inviscid flows. Therefore, and as mentioned in Section 2.2, the keypoint of vortex methods rely on the necessity of particle overlapping, that is to say, the need for the cutoff size ε to be at least of the order of h, the distance between particles.
Several solutions have been designed in order to overcome the distortion problem. One can separate these methods into two main classes.

The Meshless Rezoning Methods
The first class deals with methods involving a change in particle strength in order to correct the vorticity field such that the particle distribution better approximates the vorticity field. Among these methods one can first refer to the approach developed by Beale in 1988 [73] which consists in correcting the circulation at each time step through an iterative procedure. This technique has been successfully used by Chorin [74] to model incompressible viscous flow with a PSE scheme. However, the iterative process may be slow to converge in case of particle clustering or irregularities in the vorticity field. Both Beale and Chorin defined this circulation evaluation method as a problem of "scattered data interpolation".
In this spirit, Barba et al. [60] later proposed in 2005 a scattered data interpolation using Radial Basis Functions (RBF) as a rezoning approach in the context of vortex method. In this technique, as the vorticity ω is known only over a finite number of positions x p , the objective is to find an approximation to the assumed piecewise continuous function ω(x) as follows: where a p denote the coefficients that have to be determined by collocation using the known values of ω at x p , and where the basis functions φ p have compact support and are often chosen to be univariate and only dependent on the distance between two particles (radial functions). The solution of such interpolation requires to solve a N-linear equation (one refers to [69] for a detailed description of the algorithmic issues about such approach). The RBF interpolation is a grid-free formulation and the new particles may be located in any spatial arrangement.
Another type of meshless rezoning techniques are the Free-Lagrange methods. The latter were initially designed to model the diffusion phenomenon but they actually offer a good framework to handle distortion problems. They consist in creating a triangulation, whose nodes are the particles locations, which is then used to interpolate the circulation carried by the flow elements. For instance, the Delaunay triangulation, suggested by Russo and Strain in 1994 [75] and obtained connecting the particles when they belong to adjacent polygons in a Voronoi diagram, enables to successfully avoid the distortion of the particle distribution and provides a second-order accurate vortex scheme. Another example among these Free-Lagrange methods is the adaptive quadrature method introduced by Strain in 1997 [76]. This approach allows us to increase the order of accuracy of the usual triangulation methods by considering sub-elements that would combine several triangles together.
One can also cite the strategy introduced by Schrader et al. in 2010 [77] which can be considered as a generalization of the PSE method. It indeed proposes to correct the method given by Eldredge et al. in 2002 [78] consisting in approximating the derivatives on scattered particles locations. This corrected framework of the PSE scheme enables one to achieve the convergence without any overlapping condition and regardless of the resolution. However, the convergence of the latter was limited to a certain range of resolutions.
Finally, we can refer to the recent approach proposed by Rossi et al. in 2015 [79], called the Regular Point Distributions (RPD). This technique, introduced within the Diffused Vortex Hydrodynamics method (DVH) designed by the same authors, consists in defining an appropriated set of points on which are generated the new particles after each diffusion steps. This set of points is generated with a "package" algorithm which, in the presence of an immersed obstacle, creates points regularly around the body (of any geometry) in order to accurately prescribe the no-slip boundary condition at the surface (the enforcement of boundary conditions in vortex methods will be addressed in details in the next main section).

The Remeshing Methods
The second class of methods designed in order to remedy the distortion problem are based on a frequent physical re-location of the particles onto new positions where the overlapping condition is satisfied. These methods allow us to directly control the distance between the particles, thus ensuring the convergence of the solution. The following of this section is dedicated to the description of the re-location methods principally used by now, based on the so called "remeshing" technique.
From the pioneer works of Huberson and Jollès (1990) [80] and Koumoutsakos and Leonard (1995) [54], the remeshing approach (also called "redistributing" or "regridding") has been widely and is still successfully used in vortex methods. It consists in redistributing the particles on an underlying Cartesian grid. The description of the regridding process is here proposed for the one-dimensional case. Let us denote byx q the position of the particle q carrying the circulationΓ q (or vorticityω q ). These particles are remeshed at the grid points x k of a uniform Cartesian grid, with spacing h, using the interpolation kernel W of support [−M s , +M s ] like follows: In other words, this procedure means that every grid point k receives the sum of the weighted values of the circulation carried by the particles for which k is included in the interpolation support.
Let us recall that one of the strong advantages of the vortex methods lies in their capability to satisfy the conservation properties. In order to preserve this important feature, the interpolation kernel has to be constructed such that the successive moments are conserved. For instance, to satisfy the total circulation invariance (and more generally, the conservation of mass), W must ensure: The conservation of the first p moments is given by the following conditions: and W is therefore an interpolation kernel of order p (we refer to [22] for the proof of this assessment). For p = 0, 1 and 2, the conservation of the total circulation, the linear impulse and the angular impulse are cumulatively ensured.
In the following, we will describe the different families of remeshing kernels W by giving their construction as well as their properties.

Remeshing Kernels The Ordinary Interpolation Kernels
This group includes the kernels whose construction is only based on the satisfaction of conditions (49). These kernels, denoted Λ p where p refers to the maximum order of conserved moments, are piecewise polynomial functions of degree p. Their construction requires the calculation of p(p + 1) coefficients, obtained from the relations (49). For p = 1, one gets the linear interpolation kernel, satisfying mass conservation: As can be seen on Figure 2a where are also plotted Λ 2 , Λ 3 and Λ 4 , this type of kernels suffer from a lack of regularity. One can notice that Λ 2 is not even continuous, which is expected to lead to inaccurate solutions. In order to reduce the creation of oscillations due to the poor regularity of the kernels, Magni and Cottet [81,82] proposed to insert correction formulas in the expression of Λ p .

The B-Splines-Based Kernels
In 1946, Schoenberg [83] introduced smooth kernels, denoted M p , obtained using B-splines : For p = 1 one gets the top-hat function, M 1 . The p-order kernels M p thus correspond to the successive convolutions of M 1 : and each kernel M p is of class C p−2 ( Figure 2b). M 3 and M 5 particularly distinguish themselves for their non-oscillatory properties, as illustrated by Magni [84]. Nevertheless, the M p kernels only ensure the conservation of the two first moments, which strongly limits the accuracy of the remeshing procedure.

Extrapolation of B-Splines-Based Kernels
In order to achieve higher order remeshing schemes while preserving their smoothness properties, Monaghan [85] proposed in 1985 an extrapolation of the M p kernels. It consists in deriving linear combinations of the M p kernels and their derivatives, more precisely j , such that the successive continuous moments of W are canceled: With this approach, Monaghan obtained the so-called M 4 kernel defined by: This kernel is of class C 1 , it is piecewise cubic, and its support is of size 4 on [−2, 2]. Locally this kernel offers a third-order accurate scheme (although its global order is 2) and conserves the first three invariants when used in vortex methods (total circulation, linear impulse and angular impulse). In practice, the M 4 is a very good compromise between regularity, accuracy and computational cost, which explains that it is still widely used in vortex methods [6,8,86]. However, these extrapolated smooth kernels do not systematically satisfy the following interpolation property: For instance, the M 4 kernel does satisfy this condition, but it is not the case of the one derived from M 8 ( Figure 2c). Nevertheless, this property is essential since its dictates the algebraic conservation of the exact solution if the velocity is zero (that is to say if the particles do not move from their original location on the grid).

High Order Kernels
Much later, Bergdorf and Koumoutsakos [87] proposed in 2006 a higher order kernel, called M 6 and in 2014, Cottet et al. [88] used high order kernels in a systematic way. These high order kernels satisfy simultaneously the interpolation property, the smoothness and the moment conditions to any desired order. Their construction relies on the satisfaction of the five following properties: In [88], a new global notation is proposed for these high order remeshing kernels and stands as Λ p,r where p denotes the order to which moment conditions are satisfied and r refers to the regularity of the kernel. Figure 2d depicts some of these kernels. In particular, the kernels Λ 2,1 and Λ 4,2 respectively correspond to M 4 and M 6 .

The Directional Splitting
For multi-dimensions problems, the remeshing step is classically performed through a tensorial product of the chosen 1D kernel formula. For example, if the desired kernel has a 1D-support of S = 2M s points, then the number of operations needed to interpolate one particle on the grid is proportional to S 2 for a two-dimensional problem and to S 3 for a three-dimensional one. A new method was proposed in 2012 by Magni and Cottet [82] for the particle redistribution procedure. It consists in successively solving monodirectional sub-problems, including the advection and the remeshing step. In other words for a N-dimensional problem, one applies the following algorithm for each particle: for d ∈ {1, ..., N} : advection in direction d remeshing in direction d We refer to Figure 3 for a comparative pattern of tensorial product and directional splitting methods. As a consequence, if the chosen kernel contains S points in its 1Dsupport, the number of operations with the directional splitting method goes from O( The directional splitting consequently allows for a drastic reduction of the computational cost in terms of regridding operations. The full convergence proof of the 1D advection-remeshing scheme, considering a smooth non-constant velocity and remeshing kernels of type Λ p,r , was brought by Cottet et al. in 2014 [88], thus providing a robust and theoretical numerical framework to the directional splitting approach. The increase of the time-order of the directional splitting scheme (limited to 1 in the form described above) was proposed by Magni and Cottet [82] through the so called Strang splitting scheme, which was successfully used by Cottet et al. [88] in the context of threedimensional passive scalar advection problems. It consists in performing the directional advection/remeshing with half time steps, except for the last direction where one can gather the operations on [t n ; t n+1/2 ] and [t n+1/2 ; t n+1 ] together. The Strang splitting scheme is therefore of order 2 in time. The different substeps of this approach are given below for the 2D case (the same principle applies for the 3D case): With the Strang splitting method, a multidimensional remeshing thus costs O(3S) (respectively O(5S)) operations per particle in 2D (respectively in 3D).
The introduction of an underlying grid in order to ensure particle overlapping modifies the essence of the original vortex methods and rapidly gave birth to the so-called "hybrid" vortex methods. The last part of this section is dedicated to the description of the two main families of hybrid vortex methods, namely the Domain Decomposition method on one side and the Vortex-Particle-Mesh (VPM) methods on the other side, with the Vortex-In-Cell (VIC) and the semi-Lagrangian vortex methods among the latter.
remeshing by directional splitting remeshing by tensorial product Figure 3. Remeshing steps using a tensorial product (on top, depicted by plain arrows) and a directional splitting [82] (on bottom, depicted by dashed arrows). The red lines indicate the support of the remeshing kernel. In this example, the kernel has a 1D-support of 4 points.

Hybrid Vortex Methods
Hybrid vortex methods were developed in order to overcome the main weaknesses of the pure Lagrangian vortex methods which essentially lie in the treatment of the viscous effects (see Section 3) and, as we will see in the next section, in the near wall boundary conditions. In order to preserve the robustness of the global scheme, characterized by minimal numerical dissipation, these methods still involve a Lagrangian framework for the advection process but they also handle a fixed grid. The presence of this grid facilitates the prescription of the no-slip boundary conditions as well as the modeling of the diffusive term, using Eulerian schemes (like finite-differences, finite-volumes, finite-elements or spectral methods) while ensuring the particle overlapping condition. In the following we will consider a vortex method as hybrid whenever the resolution of the vorticity transport equation and the velocity Poisson equation will be handled both on the particles and on an underlying grid.
Before introducing some hybrid vortex methods, we need to define the common thread between all of them which stands in the way to exchange information between Eulerian and Lagrangian schemes. More precisely, hybrid methods require us, on one hand, to transfer the values of the flow fields a from the particles to the grid points (we will denote it the P2M process, for Particle to Mesh) and, on the other hand, to transfer the values of a from mesh points to the particle distribution (we will denote it the M2P process, for Mesh to Particle). It is important to note that the P2M and M2P processes do not imply a physical particle relocation, they materialize an exchange of information.
Let us first define the P2M scheme, also called the assignment scheme. It is defined by the following scheme: where V k corresponds to the volume surrounding the grid point x k (this volume equals to h 3 if the mesh is Cartesian with a uniform mesh size h) and the basis function φ k is defined by: and satisfies: If the considered grid is Cartesian, the basis function φ may be taken as one of the M p , M p or Λ r,p functions described in the previous subsection. In that case, the remeshing procedure presented previously may be considered as a particular case of the P2M.
On the other hand, the M2P scheme, also called the interpolation scheme, yields: These P2M and M2P formula lead to conservative schemes. They are depicted in Figure 4.

The Eulerian-Lagrangian Domain Decomposition Method
The first type of hybrid vortex approach presented in this section is the Eulerian-Lagrangian domain decomposition method, which was proposed by Cottet in 1991 [89] in the context of inviscid flows. The main idea relies on the fact that grid methods and Lagrangian methods are used in different parts of the computational domain. In the flow regions close to immersed obstacles, the use of grid-based schemes enables one to accurately impose the boundary conditions and on the other hand the grid-free vortex method, which is implemented elsewhere, ensures an accurate modeling of the advection-diffusion effects and also provides a simplified treatment of the domain boundary conditions. In the same spirit, Cottet et al. [90] considered an overlapping zone between the Eulerian (grid-based) domain Ω 1 and the Lagrangian (grid-free) domain Ω 2 ( Figure 5), implying a Schwarz method in the overall algorithm. It is important to note that the remeshing procedure performed in the Lagrangian domain Ω 2 is not as crucial in this context as in a pure Lagrangian framework (cf. Section 3). Indeed, the distortion effects mainly occur in the vicinity of the wall boundaries, where the highest velocity gradients are observed.

Ω1
Ω2 S2 S1 overlapping zone Lagrangian Ω 2 domains. In this case, inspired from Cottet et al. [90], Ω 1 is defined by a polar mesh and a uniform Cartesian grid lies in Ω 2 for particle remeshing.
Huberson and Voutsinas [91] also report simulations of flows past circular cylinders realized with this Eulerian/Lagrangian decomposition method, as well as the example of a wind turbine operating in a non-axial wind flow, where a finite-differences scheme is used in the near wake region. One can cite as well the works of Guermond et al. who considered matching domains (no overlapping) [92], which makes the definition of interface boundary conditions more difficult but enables a reduction of computational operations.

From Vortex-In-Cell (VIC) to Semi-Lagrangian Vortex Methods
The Vortex-In-Cell is a hybrid vortex method which was proposed in the pioneer studies of Christiansen in 1973 [93]. It relies on dealing with an underlying Cartesian grid and using either Eulerian schemes (on the grid) or Lagrangian vortex schemes (on particles), according to their respective strengths, to solve the different parts of the governing equations. This method is based on a splitting algorithm (inspired from the Chorin's viscous splitting algorithm described in Section 3.1), also called a fractional step algorithm or an operator splitting algorithm. Concretely, in a VIC method, the advection is performed by updating the positions and the vorticity of the Lagrangian particles through a set of ODEs. On the other hand, Eulerian schemes are handled on the underlying grid to efficiently compute the particle velocities and to solve the diffusion terms for which pure vortex method encounter difficulties, as widely discussed in Section 3. For three dimensional flows, the stretching term is often resolved on the mesh for the sake of simplicity.
The general Vortex-In-Cell fractional step algorithm can be described as follows: 1.

Poisson equation:
(a) Assign particle vorticity values to the grid using a P2M formula.
Compute the velocity field solving the following equation on the grid:

2.
Advection: (a) Interpolate the velocity field on the particles using a M2P formula.
Perform a Lagrangian advection of the particles and get their new positions and vorticity: Assign vorticity values and element volumes to the grid using a P2M formula.

3.
Stretching (in three-dimensions only): (a) Solve the stretching equation on the grid by differentiation of the velocity field: Note: One can also write the stretching equation in its conservative form, namely ∂ω/∂t = div(ω : u) ≡ (ω · ∇)u − udiv(ω)

4.
Diffusion: (a) Solve the diffusion equation on the grid and get the final grid values for the vorticity field: (b) Interpolate the particle values from the grid using a M2P formula.
As a typical example of study based on this type of VIC algorithm one can cite the work of Mortazavi and Giovannini [94]. In this study, dedicated to the simulation of 2D flow past thin and thick plates at Re H = 7000 (with H the thickness of the plate), the Poisson equation (step 1) is solved with a finite element method (on a mesh), the advection (step 2) is solved with a usual vortex blob approach (where the core diameter is taken equal to the mesh size) and the diffusion (step 4) is realized through a Random Walk method. Besides, let us precise that the algorithm structure written above corresponds to the original version of the VIC algorithm, following the idea of Christiansen. We note that this initial version actually suffers from the same difficulties as the one experienced by purely Lagrangian approaches because of the distortion effects provoked by the Lagrangian resolution of the advection equation (step 2). The emergence of the remeshing technique in the mid 1990s allowed us to overcome this problem and the VIC method, also known nowadays as the "Vortex Particle-Mesh method" (VPM) (or "Remeshed Vortex methods") then became particularly attractive and efficient to handle very different types of applications [95][96][97][98]. One can find a wide range of differences among all the existing Vortex Particle-Mesh methods, which can be explained by the very modular structure of the fractional step algorithm they are based on. Concerning the resolution of steps 3 and 4, a common trend is to use a finite-differences scheme (like in the four references given above [95][96][97][98]), while the Poisson equation, in step 1, if often handled with a Fast Fourier Transform method or a Fast Multiple Method (FMM) (these algorithm aspects will be explained more precisely in Section 6.1). The numerical order of a VIC method therefore depends on the order of each scheme adopted in each substep of the operator splitting algorithm. An important difference between the diverse existing VIC schemes concerns the remeshing procedure in step 2. For example, in Caprace et al. [99], the particles are remeshed on the underlying Cartesian grid only every 5 times steps, while in Mimeau et al. [100] the remeshing procedure relocate the particles on the grid points every time steps, implying the presence of a particle at each node of the Cartesian mesh, at any time. In this particular latter case, one talks about "Semi-Lagrangian Vortex methods".
Today, the Vortex Particle-Mesh methods may be considered as one of the most commonly used types of hybrid vortex methods thanks to their relevant advantages we can recapitulate in the following list:  (25)) for the whole stability, • Satisfaction of the particle overlapping condition all along the simulation through remeshing procedures, whose computational cost may be greatly decreased by the use of directional splitting schemes (see Section 4.2.2), • Limitation of the diffusive effects in the advection step thanks to the Lagrangian framework (provided the use of regular and high order remeshing kernels like M 4 or more generally Λ p,r ).
For all these reasons, the VPM methods (i.e., VIC, semi-Lagrangian) have been extensively used in DNS (Direct Numerical Simulations) to discretize the three-dimensional Navier-Stokes equations, especially in the case of highly transitional regimes with Reynolds numbers ranging from 1000 to 10,000. Among the large quantity of CFD studies realized with VPMs, one can mention for example the work of van Rees et al. [97] who carried out the 3D Taylor-Green vortex case at Re Γ = 1/ν = 1600 and successively compared their results to the one obtained with spectral methods, as well as the studies of Morgenthal and Walther [101] or Kudela and Kozlowski [102] where a VIC method is used to model incompressible flows past cylinders up to Re D = 3000 and Re D = 9500 respectively (with D the cylinder diameter). One can moreover refer to the work of Cottet and Poncet [103] which proposed one of the first 3D flow control studies using VIC, where an optimal rotation of a 3D cylinder is designed in order to reduce the drag force. In terms of flow control one also notices the investigations of Creusé et al., based on a VIC method, concerning the active control of a flow over a backward-facing step using pulsing jets [104]. We refer the reader to Section 7 for further applications handled with VIC and semi-Lagrangian vortex methods that clearly highlight the above listed advantages of these methods.
Another type of hybrid vortex method one can also cite to conclude this section is the one designed recently by Kornev in 2018 [105] where the mesh-free simulation is embedded into the grid based one. It consists in representing the large scale field on the grid, whereas the small scale field is evaluated through a pure Lagrangian vortex method, the separation of the large scale vortices from the small ones being based on a LES (Large Eddy Simulation) filtering procedure. By far, this method has only been tested on simple two dimensional canonical flows and its application to more complex physics remains to prove.

Treatment of the Boundary Conditions
The ability of correctly handle boundary conditions in a very important issue for any numerical method in CFD. As is well known, the boundaries and the flow viscosity are the only source of vorticity, which explains the efforts made in order to design a method able to accurately model the flow behavior in this region, essentially governed by the no-slip boundary condition. In the context of vortex methods, the major difficulty lies in the fact that the no-slip boundary condition only involves the velocity field in its expression whereas these methods are based on a vorticity formulation of the Navier-Stokes equations. In this section we will give a global description of the main approaches dealing with boundary conditions in vortex methods. We will therefore consider the velocity-vorticity formulation of the Navier-Stokes equations and we will restrict this survey to solid and rigid obstacles.
In the context of Lagrangian vortex methods, the objective is to express the boundary conditions in terms of the vorticity field, equivalent to the no-slip boundary condition u(x s ) = u s , where u s denotes the velocity of the solid body and x s any location on the solid surface. To write these boundary conditions one may choose a Dirichlet formulation (wall vorticity) or a Neumann formulation (wall-normal vorticity flux). The Neumann form is often preferred since it directly governs the local production of vorticity at a solid wall.
Indeed, the wall-normal vorticity flux, given by ∂ω ∂n , measures the vorticity that enters the flow from the boundaries. This flux is actually responsible for the total production of vorticity at the no-slip wall.
The approaches designed to model the boundary conditions can be divided into three main categories, exposed hereafter. The first family of methods belongs to pure Lagrangian methods; it is based on viscous splitting and mimics the vortex creation at the interface through the generation of vortex sheets. The two other families involve the presence of a mesh or a lattice and are therefore often used within Vortex-Particle-Mesh methods. The first one, known as the Boundary Elements Method (BEM) enforces the no-slip by defining vorticity flux boundary conditions on panels located at the wall, and the second one prescribes the boundary conditions by using Immersed Boundary Methods (IBM).

Boundary Treatment Based on Viscous Splitting and Vorticity Creation
The first successful approach designed to handle the boundary conditions in vortex methods for fluid flows is due to Chorin in 1973 [37], later improved by the same author in 1978 [106]. The method he introduced, called the vortex sheet/vortex blob method, uses the viscous splitting technique (described previously in Section 3.1) and relies on creating a vortex sheet at the solid boundary in order to vanish the slip velocity and thus to enforce the required no-slip boundary condition. Practically, the wall is discretized into segments. Each of these segments is associated to a local slip velocity where vortex elements are created at each time step to cancel it. We detail the different steps of the algorithm used in the improved method proposed in [106] to discretize the Navier-Stokes equations, with boundary conditions. More precisely, this approach is based on the coupling of the Prandtl boundary-layer equations near the surface with the Navier-Stokes equations everywhere else. The Prandtl equations are expressed in a coordinate system where x and y respectively denote the tangential and the normal coordinates, and where (u, v) corresponds to the velocity vector field in this system. The solid wall is assumed to be at y = 0 and the fluid domain is defined by the half-space y ≥ 0. In the boundary region, the Prandtl boundary-layer equations rely on the assumption that vorticity diffusion mainly occurs in the direction normal to the wall and high velocity gradients are essentially along the component tangential to the wall (see Figure 6); they therefore can be written as: They are completed by the continuity equation and the following boundary conditions: where Equation (66) imposes the no-slip boundary condition at y = 0 and Equation (67) expresses a far-field condition for the u component when y ≥ y ∞ . The value y ∞ is often taken of the order of the boundary layer thickness. Let us now consider a set of N vortex sheet elements S i of strengths α i . These elements should be seen as segments of a straight line, parallel to the x axis and of length ε, such that the two values of u at the extremities of S i differ from each other by α i . According to this definition, the vortex sheet approximation of ω is given by: where b ε is a one-dimensional cutoff function given by b ε = ε −1 b(x/ε) and where δ is the one-dimensional Dirac mass in the normal direction. The u component of S i , due to the presence of the other surrounding segments, is then obtained from a direct integration of Equation (65), using the right-hand side of the above ω approximation (Equation (68)): where H is the Heavyside function. The expression of v can then be deduced from u, using the continuity equation: Based on these definitions, one can give the substeps constituting the algorithm of the Chorin's vortex sheet/vortex blob method (each of these steps are depicted in Figure 7):

1.
Evaluation of the u and v components from formulations (69), (70) and advection of the vortex sheets with this velocity field.

2.
Computation of u on the wall (y = 0) and creation of new sheets at the boundary with strengths α i = −u(x i )dl in order to impose the no-slip boundary condition (these strength values directly come from the integration of Equation (65)). As explained by Chorin in [106], the velocity field is here extended across the wall by the antisymmetry u(x, −y) = −u(x, y). As a consequence ω(x, −y) = ω(x, y) and a vortex sheet is generated at y = 0.

3.
Diffusion of the vortex sheets in the y direction using a Random Walk approach (see Section 3.1). If a sheet arises in the far-field region (y ≥ y ∞ ), it is transformed into a vortex blob. If a sheet crosses the boundary, moving into the body, it is reflected on the other side. Step 1 Step 2 Step 3

Sheets convection Sheets creation
Sheets diffusion In the pioneer version of Chorin's method [37], the cancellation of velocity at the boundary was realized by generating vortex blobs instead of vortex sheets. But, the overlap of vortex blobs with the solid boundary led to inconsistency. Thanks to its more natural framework implying vortex sheets in the boundary region, this random vortex sheet/vortex blob method turned out to be widely used to model slightly viscous flow past obstacles. For example, one can refer again to the work of Smith and Stansby [107] who studied flows around cylinder, to Sethian and Ghoniem [108] who performed a validation study of this technique for flow past a backward-facing step at different regimes or to Cheer [109] who carried out simulations of flows past airfoils and blunt bodies. One can moreover cite the numerical convergence study of the vortex sheet/vortex blob method carried out by Mortazavi et al. [110] which investigates the influence of three discretization parameters (the blob strength Γ, the core radius h and the time step ∆t) on the accuracy. The role of the time step has specifically been outlined, with in particular the divergence of the numerical solution for too small values of ∆t, due to the Random Walk Method.
The three dimensional version of the vortex sheet/vortex blob method was introduced shortly after, in 1980, by Chorin [111]. In this version, the vortex sheets within the boundarylayer region are now defined as "tiles". These tiles are rectangles, parallel to the wall and their motion is dictated by the straightforward extension of Equation (65) in 3D. Like in two dimensions, the tiles undergo random walk. They are transformed into vortex blobs (also referred to as vortex segments in 3D, as opposed to vortex filaments) when they leave the boundary layer. This method was successfully used in several studies, for instance by Fishelov [112] or Gharakhani and Ghoniem [113], who respectively dealt with flow past infinite plate and flow inside a cylinder equipped with an eccentric inlet port and a driven piston.
However, the Chorin's approach suffers from several limits. The first one, already pointed out in the two-dimensional case by Anderson and Reider [114], refers to the interface zone (in which vortex sheet turn into vortex blob) where the continuity of the tangential velocity value is guaranteed but not the one of the normal velocity, which can lead to spurious oscillations close to the interface. Secondly, the divergence free condition for the vorticity field is only satisfied approximately, which prevents from getting accurate application of vortex methods to three-dimensional flows. Finally, the last and main limitation stands on the fact that the Prandtl approximation makes such approach only valid for attached flows.
As a variant of Chorin's vortex sheet/vortex blob method, one can cite the approach proposed by Walther and Larsen in 1997 [115] were the determination of the surface vorticity in step 2 is no more based on the Prandtl approximation but on the kinematic relation between the velocity and the vorticity. This method has been successfully applied in 2D to flow past a fixed and harmonically oscillating flat plate at high Reynolds number (Re c = 10 4 , with c the chord length of the plate).

Boundary Treatment Based on Vorticity Flux Conditions
The type of method presented in this section was designed later and re-use the framework of algorithms based on viscous splitting and vorticity creation. However, contrary to the Chorin's method, detailed previously, it handles vorticity flux boundary conditions to enforce the no-slip condition on the walls, instead of vortex sheets or vortex blobs. In a first part we will give the expression of the corresponding continuous problem and secondly we will focus on the numerical schemes and techniques used to approximate the integral solution of this continuous problem.

The Continuous Problem
The idea of a vorticity creation method based on a vorticity flux boundary condition was independently proposed in 1994 by Cottet [116] and Koumoutsakos et al. [117]. The so-called vorticity flux boundary condition is characterized by a relation between the normal derivative of the vorticity and the time derivative of the tangential component of the velocity. To derive the global Navier-Stokes problem, we follow the development described in [116] starting with the 2D Stokes problem. Let us consider two vorticity fields ω 1 and ω 2 in the same time interval [t n , t n+1 ]. Each time step is initialized by the exact initial vorticity field ω n and provides at the end the solution ω n+1 where ω 1 and ω 2 are the solutions of the 2 following sub-problems within the interval [t n , t n+1 ]: Sub-problem 1 : Sub-problem 2 : where s denotes the tangent vector to the surface S and u 1 corresponds to the velocity field associated with ω 1 . In other words, in sub-problem 1, the Stokes equations are solved with an homogeneous Neumann boundary condition, that is to say without creating any vortex sheet on the walls. Then, in sub-problem 2, vorticity is produced at the boundary thus enabling to cancel the slip resulting from sub-problem 1. At the end of the time step, ω n+1 is set to ω 2 (·, t n+1 ). We refer to [116] for the proof of equivalence with Chorin's algorithm and the proof of convergence using energy estimates. The generalization of the present algorithm to the Navier-Stokes problem is obtained straightforwardly by adding the advection term in Equations (71) and (72) and in replacing the vorticity flux boundary condition of (72) by:

The Integral Formulation
The issue now lies on the resolution of the continuous problems derived above in order to obtain an explicit expression of the approximated vorticity field such that the no-slip boundary condition is satisfied. The general existing algorithms are based on a viscous splitting of the Stokes or Navier-Stokes equations and may be described by the following steps:
The diffusion is solved in the domain with a PSE scheme for instance (the diffusion scheme needs to be based on a change of vorticity strength), 3.
The vorticity flux boundary condition is expressed in terms of an integral formulation derived from the continuous problem (cf. (72) or (73)), 4.
The PSE scheme in complemented by distributing this vorticity flux onto the existing vortex blobs so that vorticity enters the domain.
In the last step, particles just experience a change in their strength, they do not multiply themselves like in Chorin's vortex sheet method. In the following, we focus on steps 3 and 4.
Concerning step 3, the method designed to solve the continuous problem (72) was introduced and validated by Koumoutsakos et al. in [117] and [54]. It is denominated as boundary element method (BEM) or panels method. To present it, one re-writes the problem (72) by considering the vorticity diffusion equation supplemented with an homogeneous initial condition and a Neumann boundary condition: The integral solution of Equation (76) for the vorticity involves integrals over the surface S only. To discretize this integral, a boundary integral method is employed by defining S as a set of discrete panels, denoted with index i, length d and centered in x i . Then, the vorticity induced by a panel i can be given by the following integral: where, if one defines a coordinate system x = (x, y) with the x-axis parallel to the surface, the φ function yields: and where µ(x, t), the diffusion potential, can be approximated for points s along S by: with κ(s) denoting the surface curvature at point s and γ i = − (u · s) ∆t (s i ) corresponding to the vortex sheet strength derived from (76). We refer the reader to [117] for a full description of the integral evaluation of Equation (77).
Then one jumps to the step 4 of the algorithm, where the strength Γ j = ω(x j )v j of the particles are modified by distributing the vorticity flux: with N denoting the number of panels on the surface of the body. Finally, replacing µ i by approximation (79) in the above equation, the prescription of the no-slip boundary condition is realized through the following particle strength update: Concerning the 3D formulation, one distinguishes the flat surfaces (half-plane case) from the curved ones.
The first three-dimensional extension of the continuous formulation was given in the case of a flat plate. We note that in 3D, one must ensure that the vorticity is divergence free in order to satisfy the fundamental laws of the fluid mechanics. The direct 3D extension of problems (71) and (72) do guarantee this condition and we refer to Cottet and Koumoutsakos' book [22] for the proof. In order to solve the three-dimensional continuous problems, Casciola et al. [118] proposed a similar analysis as the one described in the above algorithm. First of all, this study supplies the expression of the integral equation for the vorticity creation. An approximation of such integral equation is then given and identified as the slip velocity at the surface of the flat body. Moreover, this work exposes an error estimate of the proposed approximation and also proves that the vorticity generated at the wall and then diffused in the flow is solenoidal.
The expression of three-dimensional continuous problems in the case of non-zero curvatures along the surface S has been first expressed for cylinders, where the curvature is constant and defined at κ = 1/R, where R is the cylinder radius. The natural coordinate system to consider herein is the cylindrical coordinate system where one denotes by (e r , e z , e θ ) the basis vectors, respectively corresponding to the unit normal vector (e r ) and the tangent vectors to the cylinder surface (e z , e θ ). According to Poncet and Cottet [119,120] it appears that the divergence free condition is ensured if the Neumann boundary condition ∂ω θ /∂n on the tangential component ω θ is modified into a Robin type boundary condition, that is to say by κω θ + ∂ω θ /∂n. The three-dimensional numerical scheme giving the solution of such formulation is derived in the same two references and the validations carried-out herein show the accuracy of the vorticity flux boundary conditions.
In conclusion, the handling of the no-slip boundary condition using vorticity flux formula stands as a suitable and successful method. It indeed allows us to treat in an accurate way the high velocity gradients at the solid-fluid interface while overcoming the main drawbacks of the prior approaches. It is important to highlight the fact that the BEM technique is generally related to the presence of a lattice, first because of the use of a PSE scheme (that needs particle overlapping to ensure convergence, as evoked in Section 3.2) and most of all, because its design is based on the definition of panels on which the flux integrals are discretized. BEM is therefore very often used within Vortex-In-Cell/Vortex-Particle-Mesh algorithms, but not systematically as illustrates the completely mesh-free proposition of Cooper and Barba [121] based on the use of the radial basis function (RBF) collocation method (see Section 4.1). The latter method evaluates the φ function of Koumoutsakos (Equation (78)) and solves the diffusion equation by using the radial basis functions and therefore proposes a treatment of the boundary conditions without any panel.

Improvements of Vorticity Flux Boundary Conditions
Based on the pioneer works of Koumoutsakos and Leonard presented above, Ploumhans and Winkelmans proposed substantial improvements in [122] (2D study) and in [55] (3D extension) in the early 2000s. These advances may be separated into two categories.
First, they proposed to improve the step 2 described previously in the viscous splitting algorithm, concerning the diffusion process. Indeed, the authors point out that if the standard PSE scheme is used, namely: then the particles j are completely and radially surrounded by other particles i, which is not true close to the wall. Thus, the use of standard PSE scheme induces spurious vorticity flux at the boundary (even if the total circulation is conserved). To overcome this problem, Ploumhans and Winkelmans [122] introduced a modified PSE scheme in which image particles are defined inside the body when the scheme is applied close to the wall. This process allows a zero vorticity flux during the PSE step, before computing it properly at step 3. The second type of improvements are based on the fact that BEM (or panels methods) rely on grids to perform the particles remeshing. In Koumoutsakos and Leonard [54], the BEM was proposed in the particular case of a circular cylinder and the particle redistribution was performed on a body-fitted lattice, that followed the cylinder boundary, making it quite restrictive in terms of applications. In their work, Ploumhans and Winkelmans [122] propose on one side a regular lattice able to redistribute the particles in the case of any general geometry for the immersed obstacle, and on the other side, they introduce a mapping on non-uniform lattices, thus allowing a coarser and controllable resolution for the far wake regions.
All these advances have been validated in the 2D case [122] at different flow regimes for flow past a cylinder, a square and a generic geometry (2D "Apollo" capsule) as well as in the 3D case [55] with the simulation of flow past a sphere at transitional Reynolds numbers (from Re D = 300 up to Re D = 1000).

Boundary Treatment Based on Immersed Boundary Methods
Driven by the necessity to overcome the difficulties related to the body geometry, the Vortex-Particle-Mesh methods turned, in the early 2000s, to another strategy able to enforce boundary conditions without being constrained neither by the type of underlying mesh nor the obstacle geometry. This strategy lies in their coupling with Immersed Boundary Methods (IMB).
One of the first types of Immersed Boundary Vortex methods to be introduced was an approach in which the boundary conditions are in the form of local forcing terms in the right-hand side of the momentum equations. In this approach, the flow equations are discretized in a unique way in the computational domain, independently of the nature of the cells (i.e., fluid, solid or mixed cells). One can particularly cite the approach proposed in 2003 by Cottet and Poncet [120] in the context of a VIC method, where the no-slip boundary condition is enforced from source points that are located on the boundary itself on all grid points that are at a distance less than the grid-size from the boundary. Few years later, Poncet [123] proposed a revisited immersed boundary technique allowing to satisfy the no-through-flow boundary condition on a body of arbitrary geometry, in an implicit formulation. More precisely, the approach consists in solving the linear system associated to the source-to-flow-through linear application; this system is of reasonable size since it is only constituted by the discretization points of the surface, not all the grid cells. In Poncet's paper, this vortex IBM is applied to flow past a complete aircraft geometry and to the vorticity evolution around the complex geometry of a hollow bridge structure.
From this first approach arose other immersed boundary vortex methods, presented hereafter.

The Brinkman Vortex Penalization Method
The Brinkman penalization technique is part of the continuous forcing immersed boundary methods. It was first proposed in a context completely detached from vortex methods by Caltagirone in 1994 [124], further developed by Angot et al. [125] and was expressed and used in a vorticity formulation by Kevlahan and Ghidaglia in 2001 [126]. The resulting Brinkman-Navier-Stokes vorticity equations are: where the forcing term ∇ × [λχ(u s − u)] is the penalization term, in which χ denotes the characteristic function of the solid bodies, u s corresponds to the velocity of the immersed obstacle and λ is the non-dimensionalized penalization factor. This penalization term actually comes from the Darcy equations that govern flows in porous media, and the factor λ is directly related, in the inverse proportion, to the permeability k of the considered media in the domain, i.e., λ ∼ 1/k. The penalization factor λ is set to zero in the fluid region and to infinity in the solid ones. With such approach, the force is activated only inside the solids, making the velocity vanishing in these regions, and the no-slip boundary conditions are automatically satisfied. The condition u = u s is thus prescribed on the boundary and rigid motion is ensured within the solid.
In the 2010s arose the first resolutions of the vorticity-penalization Equations (83) and (84) through VPM methods. In practice, the penalization equation ∂ω/∂t = ∇ × [λχ(u s − u)] is added to the fractional step algorithm and is solved on the grid, often using finite-differences and an implicit Euler scheme. In the last decade, this technique has been widely applied within VPMs to handle different types of problems involving boundaries. We can for instance cite the works of Gallizio [96] and Morency et al. [127] who coupled the penalization technique with a level-set method in order to handle moving bodies. Besides, in terms of fluid-structure interactions we can cite the works of Coquerelle and Cottet [128] who also coupled the penalization method with a level-set method and moreover with a collision model to simulate the interaction of an incompressible flow with rigid bodies. One also has to mention all the studies of Gazzola et al. focusing on obstacles' motion induced by vortical flow-structures [10] or dealing with fish schooling [6][7][8]86], where the penalization method is combined to a projection approach allowing to model the action of the fluid on the deforming fishes (more details are reported in Section 7 about these applications).
In terms of numerical advantages, the vortex penalization method naturally exhibits the ones of the immersed boundary methods, that is to say the absence of explicit prescription of the boundary conditions, the existence of a unique equation for the whole domain Ω and the possibility to handle complex geometries. Moreover, the Brinkman penalization method is particularly easy to implement and to parallelize. Finally, it allows the direct modeling of porous medium when giving intermediate values to the penalization factor, 0 < λ 10 3 . Figure 8 depicts and summarizes the modeling of flow in solid-fluid-porous media with the Brinkman penalization method and we refer the reader to Section 7 for an example of related application. However, the major drawback of this approach lies in its first-order accuracy. Besides, performing relevant simulations imposes a restrictive overall time step ∆t in the fractional step algorithm. To remedy these weaknesses, Hejlesen et al. [129] introduced an iterative penalization scheme allowing to accurately enforce no-slip boundary conditions through sub-iterations of the penalty method. The proposed scheme imposes the resolution of a Poisson equation at each iteration to re-evaluate the penalized velocity, but it is applied only in the vicinity of the body, which on one hand reduces the computational cost of each Poisson subproblem and on the other hand allows us to consider larger time-steps within the iterative penalization loop than the ones usually required in a standard Brinkman penalization. This iterative approach has proved to be more accurate for all the classical two-dimensional benchmarks tested by the authors (flow past an impulsively started cylinder at Re D = 9500 and flow past a flat plate at different regimes). Another version of iterative penalization method was also proposed by Gillis et al. in [130] where they formulate the penalization problem as a linear system that is solved iteratively using a recycled Krylov subspaces.

The Immersed Interface Vortex Method
Despite the improvements brought by the iterative Brinkman penalization methods described previously, the latter keep on smearing the interface of the immersed body and are still first order convergence methods. To overcome these weaknesses, a third type of immersed boundary technique was proposed in 2016 by Marichal et al. [131] and in 2019 by Gillis et al. [132] to be coupled to VPM: the immersed interface method (IIM). In a general point of view, the immersed interface method (originally proposed by Leveque in 1994 [133]) consists in treating the flow discontinuities at the body's surface through finitedifference stencils, where the jumps are precisely taken into account inside the stencils thanks to modified Taylor series. These jump-corrected Taylor series allow us to evaluate the derivatives of the flow fields, which remain valid across the interface. The prescription of the desired boundary conditions is then realized by injecting them directly in the jumpcorrected Taylor series discretizing the derivatives of the velocity field at each control point of the body's surface. In Gillis et al. [132], this approach is embedded in a two-dimensional VPM algorithm and tested in the case of flow past a 2D impulsively started cylinder at high Reynolds numbers (up to Re D = 40,000). On one side the results show a second order spatial accuracy, up to the boundary region, and on the other side they reveal to be in good agreement with previous VPM computations based on Brinkman penalization performed by Rossinelli et al. [134], while considering an eight-fold coarser grid size near the wall.

Algorithmic Issues
All the sections written so far enabled to outline the landscape of the different families of Vortex Method, with their respective characteristics in terms of numerical convergence, accuracy and reliability. In this part we will now focus on the algorithmic and computational techniques developed to efficiently handle them. This point constitutes a crucial issue for any numerical method which aspires to be considered as a usable and challenging one.

Fast Vortex Methods
Regardless of which type of Vortex Methods one considers, it is admitted that the intrinsically most expensive operation in any vortex algorithm relies on the resolution of the Poisson problem. Indeed, for N particles, the classical discretization of the Poisson Equation (15) through the Biot-Savart Law (17) implies to compute the interactions between each particle and all the others, leading to a O(N 2 ) computational cost. A wide range of studies have been dedicated to the task of reducing this cost. They are mainly split into two categories: the strategies dedicated to the Vortex-Particle (Lagrangian) framework and the one dedicated to the Vortex-Particle-Mesh context.
Among the algorithms able to solve the Biot-Savart law, the most efficient one is certainly the Fast Multipole Method (FMM). This technique, introduced by Greengard and Rokhlin in 1987 [135], exploits in the context of vortex methods the geometrical distributions of the vortices to evaluate the velocity. The key point is that the velocity field induced by a group of particles clustered around a center does not need to be computed directly from its individual members, but it can be approximated by a finite number of multipole expansions. This approach, based on tree data structures, thus allows us to achieve a O(N) computational cost. Parallel multicores FMM implementations have been widely developed and used in the context of Lagrangian Vortex methods, like in the works of Marzouk and Ghoniem [136], Yokota et al. [137] or Giannopoulou et al. [138], to deal with the velocity resolution. Besides its capability to accelerate the evaluation of particle interactions, FMM was also used by Hu et al. [139] to reduce the time consumption of the stretching term calculation in the vorticity equation. However, the use of FMM algorithms is not restricted to Lagrangian Vortex approaches and was also implemented within Vortex-Particle-Mesh methods, as it will be illustrated in the last Section 7.
As one evokes the Vortex-Particle-Mesh methods (VIC, semi-Lagrangian), they particularly and naturally focused on fast grid-based Poisson solvers in order to take advantage of the presence of an underlying mesh. The latter are called fast Poisson solvers, as opposed to the traditional "slow" solvers based for instance on Jacobi, Gauss-Seidel, Conjugate Gradients or Successive Over Relaxation (SOR) methods. These fast Poisson solvers can be defined as hierarchical methods and among them, one can distinguish the direct solvers from the iterative ones. The most popular direct solvers are the Cyclic Reduction method, consisting in discretizing the Poisson equation by means of finite-differences method over a rectangle by a successive splitting of the initial problem, and the Fast Fourier Transform method, which consists in solving the Poisson equation in the Fourier space using FFTs. These two direct grid-based Poisson solvers present a O(N logN) cost and have been widely used in VPM, for instance by Cottet and Poncet [120] for the Cyclic Reduction method (based on the Fishpack solver), or by Mimeau et al. [100] for the FFT approach (based on the FFTW fortran solver), where the computational cost of the Poisson problem was reduced to 5% within a whole time-step. On the other hand, one of the most famous iterative solvers relies on the Multigrid method, which recursively uses coarse discretizations of the Poisson problem to approximate its solution, following V-cycles. The complexity of Multigrid algorithms is O(N), which makes them very attractive in terms of computational cost. However, they may require more efforts in terms of implementation and use.

Performance-Computing Improvements Using GPUs
Because of their close link with computer graphics, the smooth particle hydrodynamics (SPH) methods were the first particle methods to have been implemented on graphic processors (GPUs) in the mid-2000s. Few years later, in 2008, arose the first implementation of Vortex methods in GPUs with the work of Rossinelli and Koumoutsakos [70]. In this pioneer contribution, the implemented 2D solver relies on a remeshed vortex method and runs exclusively on the GPU, including a periodic 2D multigrid Poisson solver developed for the GPU. This work, validated on 2D vortical flows was further successfully tested on 2D bluff body flows with vortex penalization techniques and a CUDA FFT library (cuFFT) by Rossinelli et al. [140]. It shown a good agreement and, above all, important speed-up with respect to identical simulations performed with prior solver based on CPU (Central Processing Unit) (up to 30 times speed-up of CPU-computations for flows past an impulsively started cylinder). In terms of GPU-simulations based on remeshed vortex methods, one can also cite the very recent work of Keck [141] who managed to solve a three-dimensional sediment flow on a 1537 × 512 × 512 Cartesian grid with a unique 32 gigabits Nvidia® Tesla V100 node, within only 6 h (this work will be presented into more details in Section 7).
Regarding mesh-free approaches, Yokota and co-authors also developed in the late 2000s a FMM-GPU algorithm (inspired from the one first proposed by Gumerov and Duraiswami in 2008) to realize Lagrangian vortex method simulations of 3D isotropic turbulence [137]. Few years later, they proposed in [142] a comparison between a spectral method and a FMM-based vortex method on 4096 GPUs for the simulation of 3D homogeneous-isotropic fluid turbulence with 69 billions of grid nodes or particles. The authors highlight the limits of spectral methods in simulating problems on such a great number of GPUs and show that, in such configuration, the parallel FFM-GPU solver allows us to achieve a much better computational performance and efficiency than spectral methods, due to the non-local FFT-based operations of the latter, and show the relevance of advanced Lagrangian vortex solvers for the simulation of isotropic turbulence.

Adapted Grids, Multiresolution Aspects and LES Models
The Lagrangian framework of Vortex-Particle methods intrinsically allows us to spatially adapt the particle resolution according to the flow physics (provided the overlapping condition is satisfied). One can refer to the recent work of Giannopoulou et al. [138], based on the DVH method (already evoked in Section 4.1), where a self-adaptive particle resolution (also called r-adaptivity) is carried out in order to progressively reduce the particle size when getting closer to the boundaries. The method is tested in the cases of two-dimensional flows past cylinders with circular, elliptical and triangular sections, and successively compared to a Finite-Volume method. However the results reveal a higher computational cost with the DVM approach, which can be in particular explained by the large number of vortices required at the boundary compared to the number of cells required in Finite-Volume methods as well as the smaller time steps required by the explicit DVH approach.
Concerning hybrid vortex methods (i.e., Eulerian/Lagrangian domain decomposition or VPM methods like VIC and semi-Lagrangian vortex methods), they have encountered, as shown throughout this paper, a growing notoriety since the late 1990s and rapidly got concerned by the cost of the grid involved in their schemes. Like in Eulerian methods, questions related to adapted grids, multiresolution simulations or subgrid scale models therefore naturally arose in order to give VPM methods the opportunity to handle complex physical problems and/or turbulent flows, which require a high level of mesh refinement and thus a high computational cost.

Adaptive Mesh Refinement (AMR)
This technique consists in defining blocks of uniform grid sizes in a dynamic way from a posteriori estimates. It has been first introduced within a 2D remeshed vortex method by Bergdorf et al. in 2005 [143]. In this work the two-dimensional grid adaptivity obeys a criterion based on the velocity derivatives in order to refine the underlying mesh in the flow regions of high gradients. The simulation results are compared to those of a r-adaptive vortex method with global mappings for the inviscid evolution of an elliptical vortex of compact support in an unbounded domain. It turns out that AMR-based method seems to be usable in a wider range of problems since the refined regions can be created on demand based on any criterion (provided a stepwise refinement transition from coarse to fine grid-sizes), contrary to the method based on adaptive global mappings that must adapt smoothly to the flow physics involved in the problem. Rasmussen et al. [144] also proposed a 2D multiresolution VIC algorithm with Brinkman penalization and validated it in the case of classical benchmarks involving bluff bodies; in this study the local refinement is not adaptive but fixed a priori. In 2010, El Ossmani and Poncet [145] introduced a 3D adaptive mesh refinement based on multilevel sub-blocks that adapt to the flow evolution through the periodic evaluation of the quantity α|ω| + β|∇ × ω|. They tested it to study 3D ellipsoid vortex dynamics and the wake development behind a sphere.
Based on these considerations, Bergdorf and Koumoutsakos [87] proposed in 2006 another multiresolution criteria for the grid adaptivity within remeshed vortex methods: the wavelets. The decomposition of a signal into wavelets (i.e., wave-like oscillations) is a mathematical tool complementary to Fourier transform. It is often used to extract information from data and to design compression/decompression algorithms. In the context of multiresolution analysis, the wavelet theory is applied to extract information from velocity and vorticity fields and to induce local refinements or compressions of the grid.
More precisely, the wavelet-based adaptivity proposed in [87] is performed using biorthogonal wavelets. The latter introduce a pair of analysis functions (involing the Fast Wavelet Transform) and a pair of synthesis functions (involving the inverse FWT). The velocity and vorticity field can therefore be decomposed in "scaling" and "detail" coefficients at different levels of resolution. Each point of the adaptive grid is thus represented by a scaling and a detail coefficient. The aim is to discard the detail coefficients which, compared to a fixed threshold, do not carry relevant information about the flow field. For each level of resolution, one only keeps the coefficients called active coefficients. To build up the wavelet-based AMR method, a block structure is used. Each block contains active coefficients situated at the same level of resolution and all blocks contain exactly the same number of active coefficients. These blocks are built so that they can be split into finer (refinement) or collapsed into coarser blocks (compression). In order to capture the emerging small scales, refinement of the grid is performed, written R(ω n , u n )| t r , based on a refinement threshold t r at the beginning of each time step through inverse FWT. For compression, noted C(ω n+1 , u n+1 )| t c , a threshold t c is used at the end of each time step to truncate terms in reconstruction, retaining active coefficients and discarding detail coefficients which do not carry significant information. The blocks where these detail coefficients are located then collapsed into coarser ones performing FWT. Figure 9 depicts this process. The concept of wavelet adaptivity has been later implemented by Bergdorf, Koumoutsakos and co-authors in a 2D remeshed vortex library, called MRAG-I2D [134], where the space-adapted grids are also coupled with Local Time-Stepping (LTS) integration schemes in order to speedup the computational efficiency. These schemes, exploiting the time step locality imposed by the different stability conditions, allow us to integrate the coarser elements with larger time steps than the finer ones and thus enable one to perform fewer integration steps. Local time-stepping schemes are shown to speedup the computation by one order of magnitude at each level of resolution from a certain number of levels. The non-uniformity of the grid imposes the use of a Fast Multipole Method (FMM) to solve the velocity Poisson equation. The MRAG-I2D library efficiently maps these different algorithmic and numerical aspects onto heterogeneous CPU-GPU platforms. It has been used for the simulation of bluff body flows with Reynolds numbers ranging until 40,000, flow-mediated interactions between two solids, analysis of optimal shapes in anguilliform swimming and reinforcement learning in the context of self-propelled bodies (these applications will be further detailed in Section 7).

Refinement step Compression step Vortex Particle
Mesh algorithm Figure 9. Sketch of refinement/computations/compression process within one time step of Vortex-Particle-Mesh method with wavelet-based Adaptive Mesh Refinement (AMR).

Subgrid Scale Vortex Methods
When the flow under study becomes turbulent, and therefore computationally expensive to simulate, the Large Eddy Simulations (LES) allow us to only capture and resolve the most energetic scales while modeling the effect of the Sub-Grid Scales (SGS) on the captured scales [146]. Before presenting some SGS vortex methods, let us first briefly introduce the general concept of LES methods. As evoked above, this approach consists in only resolving the smallest scales down to the cutoff length ∆ of a chosen filter F . Typically the spacial discretization is on the order of ∆. Thus the resolved (or filtered) quantityQ is defined as: and the remaining smallest scales must be modeled. The filtered Navier-Stokes equations, written here in their velocity-pressure formulation, therefore can be written: where τ ij = u i u j −ũ iũj , referred to as the subgrid scale stress tensor, models the smallest scales behavior in the flow. Among the different classes of subgrid scale models, the eddyviscosity models are the most commonly used as for instance the famous Smagorinsky model, given by: where C is a fixed constant, called the Smagorinsky coefficient, and |S| is the rate-ofstrain tensor. Vortex methods allow us to carry out under-resolved viscous simulations and large eddy dynamics. To illustrate this statement we will describe the guidelines of three different approaches.
The first one is the anisotropic subgrid scale model, based on Vortex method, proposed by Cottet in 1996 [147]. The 2D model derivation is based on the resolution of the regularized/smoothed Euler equations using vortex blob method (see Section 2.1): where the function ζ ε = ζ(x/ε) ε 2 results from a cutoff function ζ, that one assumes positive and decaying at infinity, of mean value 1 and with radial symmetry. The approximation of the exact vorticity field ω is then defined by the regularization ω ε = ω ζ ε . In the spirit of LES methods, at a finite core size ε it is appropriate to interpret ω ε as the filtered vorticity, and the core smoothing function ζ ε as the spatial filter, of size ∆ ≡ ε. Performing the convolution of (88) with ζ ε one gets: where E represents the truncation error resulting from the fact that the convolution was only applied on ω in the advection term div(u ε ω). Based on the expression of error E, one can then derive an evaluation of the enstrophy production, given by: where the subscript ε has been dropped for clarity. This expression indicates that the enstrophy will increase because of the exchange of vorticity between point x and y when they will satisfy: meaning that the truncation error E(x) in the mollified Euler equations intrinsically embed antidiffusive (or backscatter) mechanisms. The idea proposed in [147] is therefore to clip antidiffusivity while leaving contributions in the direction of dissipation untouched, leading to the following 2D anisotropic eddy viscosity model: where the index-denotes the negative part of the quantity between brackets. This diffusion scheme, very reminiscent of PSE methods, is conservative. As explained by Cottet in [147] and in all the following works of the same author related to this subject, the proposed eddy viscosity model turns out to be more accurate than usual Smagorinsky-type models, which would produce stronger diffusion for identical flow problems. A few times later, Mansfield et al. [148] proposed another development of a LES vortex scheme in order to model the effects of subfilter scale (SFS) velocity and vorticity fluctuations. To do so, the authors start from the 3D filtered vorticity transport equation: where( ·) denotes the low-pass filtered quantities and is the subgrid scale vorticity stress. ∂R ij /∂x j refers to the subfilter scale torque g. This SFS torque is modeled based on a Smagorinsky-like eddy diffusivity scheme: with ν sma = (C∆) 2 |S| the Smagorinsky eddy diffusivity. The authors then propose a particle discretization of this SFS model based on a vortex blob discretization and introduce a dynamic variation of the model coefficient C, directly evaluated in a Lagrangian way from the particles data. Based on this approach, the same authors performed dynamic LES of colliding three-dimensional vortex rings in [149]. The results coincide with experimental observations and highlight in particular the capability of the model to correctly simulate the generation of small-scale turbulent structures as well as the formation of small-scale ringlets propagating radially after collision. The same approach was also used by Pinon et al. [150] (with a fixed model-coefficient C) in order to perform LES of a 3D flow of a round-jet in a crossflow, providing results in good agreement with previous experimental and numerical studies.
With the increasing use of Vortex-In-Cell and semi-Lagrangian vortex methods, a novel way of performing LES within vortex methods appeared, where the subgrid scale vorticity stress in the rhs of the filtered vorticity transport equation is directly solved on the grid, using finite-differences schemes. This is what proposed Cocle et al. [151] who investigated several types of existing SGS models within a VIC method. Among the tested models, the usual Smagorinsky model is not able to be activated only during the complex phases of the flow, and turns out to be too diffusive, as expected. A good compromise highlighted in [151] is the use of the so called "regularized variational multiscale" model. In this approach, the artificial diffusion does not operate on the complete LES vorticity field, as in usual Smagorinsky models, but acts solely on the high frequency part of the LES vorticity field: with ω s denoting the "small-scale" vorticity field defined as ω s = ω −ω (whereω is the filtered vorticity obtained through a discrete filter of size ∆) and where the subgrid viscosity ν sma is evaluated through a Smagorinsky model (see definition (87)). As reported in [151], this multiscale subgrid model included in a VIC method allows us to perform large eddy simulations with good spectral behavior. It indeed preserves the low and medium wavenumbers while ensuring the SGS dissipation at the high frequencies. This SGS model has been further used within others VPM methods, to handle challenging applications in aerodynamics (some of them are presented in Section 7). All the algorithmic strategies exposed in this section, as well as the cited related works, are summarized in the diagram of Figure 10.

Applications and Issues
Previous section illustrates the capability of purely meshless vortex methods to excel in the calculation of isotropic turbulence. They also proved their ability in simulating simple geometrical bluff body flows. Besides these types of applications, they seem to show a certain difficulty in handling complex flow physics implying the coupling of the Euler or Navier-stokes equations with other models.
On the contrary, hybrid vortex methods like Eulerian/Lagrangian domain decomposition or VPM methods, continually succeeded during the last two decades in undertaking more and more challenging applications. This last section proposes a non-exhaustive presentation of recent and advanced simulations performed with hybrid vortex methods in order to illustrate their capacities and their versatility. Table 1, given at the end of this section, outlines the main physical and numerical characteristics related to each of them.

Hydrodynamics
The simulations exposed hereafter deal with Vortex-Particle-Mesh methods used to study fluid-structure interactions in hydrodynamics and more specifically in fish schooling. The involved Reynolds numbers (based on the fish length L) stand in the range 100 ≤ Re L ≤ 5000.
Propelling and energy harvesting of articulated bodies (referred with number (#1) in Table 1) by Bernier et al., 2019 [9]. This study examines the hydrodynamics resulting from passive and actuated bodies in fluids. More precisely it consists in handling the passive propulsion of a free swimming articulated fish within the wake of a leading bluff body. The free swimmer is modeled by a multi-body articulated system including complex kinematic constraints like rotational joints and linear rails for self-propulsion as well as damper-like elements for energy harvesting. The parametric simulations show an efficient combination between power production (for the damped passive multi-body swimmer) and flow control (for the leading body).
Optimization of fish schooling (#2) by Gazzola et al., 2012b [86] and van Rees et al. 2013 [7]. The simulations performed in these two studies both include a CMA-ES (Covariance Matrix Adaptation Evolutionary Strategy) optimization algorithm within their VPM method in order to study optimal swimming regimes for anguilliformed robots. For the modeling of the swimming bodies, both study involve independent local deformable grids on which are defined the fish shapes and carrying the values of deformation velocities and characteristic/level-set functions, then interpolated on the underlying Cartesian computational grid. In [86], the CMA-ES algorithm is used to identify the fish shape-pattern which maximizes the escape distance. The resulting best motion kinematics exhibit the C-shape pattern as the optimal one, as confirmed by in vivo observations. In [7], the CMA-ES procedure is applied to the identification of the fastest and most efficient shapes for self-propelled swimmers. The results provide achievable robotic-shapes that allow us to outperform real larval zebrafish geometries by 40% and 135% for speed and efficiency, respectively.
Reinforcement learning for self-propelled swimmers (#3) by Novati et al., 2017 [152]. Another AI model-free procedure used by Novati et al. within a VPM method to optimize fish swimming lies in reinforcement learning (RL) algorithms. Their simulations involve two fishes, swimming in a synchronised tandem configuration (leader-follower). The coupling of a VPM with a Q-learning (RL) algorithm allows us to automatically define an optimal policy allowing the follower to maximize a desired and specified long-term reward (like the reduction in energy spent).

Aerodynamics
Vortex methods are known to be intrinsically well suited for advection-dominated flows. Problems related to aerodynamics are therefore natural applications for particle approaches. This subsection presents some studies showing the capability of hybrid vortex methods to perform LES of turbulent aerodynamic flow, characterized by Reynolds numbers higher than 10 4 .
Dynamics of the wake of aircraft wings (#4) by Cocle et al., 2008 [14]. In this study the LES simulations, performed from a VPM method based on the SGS model (96), handle the 3D space-development of wake vortices behind an elliptically loaded wing in ground effect at Re c = 10,000. They clearly depict the complex physics occurring in the wake region between primary and secondary vortex and the apparition of small-scales vortices due to the break down of the primary vortex under the effects of the secondary vortex spiral motion.
Dynamics of rotorcraft wakes (#5) by Caprace et al., 2020 [12]. Large eddy simulations of a 3D flow past an advancing helicopter rotor are realized through a VPM approach based on SGS model (96) and coupled with an Immersed Lifting Line (ILL) method. The latter, derived from the algorithms commonly used in wind energy community, is fully compatible with the Vortex-Particle-Mesh framework. It indeed consists in computing the vorticity corresponding to the local forces developed by the modeled wing/rotor and then in introducing this vorticity in the flow using vortex particles, that is to say in a Lagrangian fashion. The simulations performed in this study allowed, in particular, to precisely describe the vortex dynamics involved in the generation of the two main wake vortices of the rotorcraft thanks to thorough 3-D visualizations. For further details about the ILL approach in a VPM framework, one also refers to [99], by the same authors.
Dynamics of rotorcraft and finite aircraft wakes (#6) by Stock et al., 2010 [13]. Another type of large eddy simulations of rotorcraft wakes was earlier proposed by Stock et. al using a hybrid Eulerian/Langrangian decomposition domain approach (see Section 4.3). In the near-body region the computation is realized through a fully-compressible, Eulerian, Navier-Stokes flow solver based on a body-fitted structured grid. In the wake region, a pure Lagrangian vortex method is used, with the following characteristics: the velocity evaluation is performed through a GPU-FMM solver; since the Reynolds numbers considered in this study are in the order of Re c ≈ 10 6 , no diffusion scheme is applied; a boundary element method (BEM) is used to set on each triangular panels the bound vortex sheet strength involved in velocity computation; an anisotropic SGS model based on Cottet's model [147] (3D extension of Equation (93)) is implemented for LES; finally, an overlap region and a coupling algorithm is carried out at the Eulerian/Lagrangian interface. LES of finite aircraft wing wakes are first performed, showing results that coincide rather well with experiments (showing however discrepancies in the region of the tip-vortex core), before presenting LES of 4-blades advancing rotor wakes for two angles of rotation of the whole device in forward flight.

Scalars Transport
Scalar advection problems involve the coupling of the Navier-Stokes equations with an advection-diffusion equation for the scalar transported at the flow velocity. In such problems, the production of small scales is dictated by the value of the Schmidt number, Sc = ν/κ, defined as the ratio between the viscosity of the fluid and the diffusivity of the scalar. If Sc > 1, the size of the smallest scalar fluctuations is smaller than the smallest scale of the turbulent flow. In such situation, it appears quite natural therefore to perform numerical simulations involving different grid resolutions for the scalar and the momentum. The two following studies showcase the capability of VPM methods, more precisely here semi-Lagrangian vortex methods (see Section 4.3), to efficiently handle such problems.
Multi-scale solver for scalar turbulent transport on CPU-GPU (#7) by Etancelin et al. 2014 [153]. The computing strategy of this study is based on a semi-Lagrangian vortex method and consists in solving the Navier-Stokes equations (momentum and Poisson equation) on multi-CPUs while the scalar transport is computed on several GPUs. Passive scalar in a turbulent plane jet of Sc = 10 and Re w = 1000 (with w the jet width) is simulated in a hybrid CPU/GPU computation with 128 3 and 512 3 grid points for the velocity and the scalar, respectively. In [154] a thorough HPC study is performed by the same author for turbulent jet problems (see Figure 6.13, p. 133) varying the respective grid sizes of the flow quantities and the scalar and varying the number of CPUs and GPUs in the parallel computation. It shows in particular that the computational time for the resolution of the scalar on a 1024 3 grid on 8 GPUs takes half the time needed for the resolution of the flow on a 256 3 grid on 32 CPUs.
Multi-scale solver for sediment transport on GPU (#8) by Keck 2019 [141]. This study examines the three-dimensional dynamics of sediment-laden fresh water above salt water. The modeling of such problem involves the coupling of the Navier-Stokes equations with two transport-diffusion equations concerning respectively the salinity S and the sediment-particles concentration C. The Schmidt numbers considered in this study for the concentration C are Sc (C) = 3.5, 7, 14, 28 and the ratio between the two scalars diffusivities is equal to 25. The problem is solved via a semi-Lagrangian vortex method, using GPU based computations. The obtained 3D sediment isocontours show, as expected from literature, salt fingers initially present on both sides of the sediment-laden water/salt water interface that then disappear under Rayleigh-Taylor instabilities for Sc (C) = 3.5 and 7. For the two other Schmidt numbers Sc (C) = 14 and 28, no more sediment fingers can be seen. They are replaced by downward moving sediment plumes of mushroom-shape, characteristic of Rayleigh-Taylor instabilities (see Figure 11 (a1 to a4)). These computations are performed on a lonely GPU core. An extension of this study to multi-scale hybrid CPU-GPU computations is envisioned.

Porous Flows
The last examples given in this non-exhaustive list of applications focus on the ability of VPM, and more precisely hereafter of semi-Lagrangian vortex methods, to handle porous flows. This can be done through the coupling of the governing incompressible flow equations with Darcy or Brinkman equations, which are solved on the grid.
Passive flow control using porous media (#9) by Mimeau et al., 2014 and2017 [20,155]. Based on the results obtained in the 2D case [155], the study [20] proposes to simulate 3D incompressible flows in solid-porous-fluid media through a semi-Lagrangian vortex method and then to perform passive control of bluff body flows by partially covering the body surface with a porous coating. The modeling of the physical problem is directly realized through the Brinkman-Navier-Stokes equation (Equation (83) The results show that the effects of the added porous layer on the wake stabilization and drag reduction noticeably depend on the permeability and the location/geometry of the covering layer. The parametric simulations allow us to identify an efficient "porous ringinlay" control device whose porous coating is only located in flow separation areas (see Figure 11 (b1 to b3)).  [156]. Contrary to all the previous studies reported in this review paper, this last illustration of semi-Lagrangian capabilities deals with Darcy-Brinkman-Stokes equations, that is to say with problems not dominated by advection. More specifically, this last study aims to simulate dissolution processes at the pore scale of actual rocks. To do so, the authors couple a Lagrangian formulation of the chemical reaction equations (handled through a semi-Lagrangian solver) with the so called superficial velocity formalism (handled with an Eulerian scheme). The authors address the crucial issue of how to estimate the numerical diffusion induced by the remeshing step in the semi-Lagrangian solver and how to switch from a low-order to a high-order kernel in order to correctly adapt the physical diffusion. The 3D DNS simulations highlight the effects of fluid inclusions in a non-percolating and percolating carbonate rock.

Conclusions
The intensive development of the vortex methods between 1970 and the mid 1990s was motivated by the will to design robust and accurate numerical schemes able to overcome the main weaknesses of traditional grid methods, especially in the case of dominant advection effects in the flows. The major advantages arising from these pure Lagrangian methods are the following: appealing physical attributes, automatic adaptivity of the discretization elements (particles), compactness of vorticity support, accurate handling of boundary conditions at infinity, negligible numerical dissipation, conservation of the flow invariants, relaxed stability condition on the time step and usefulness of vorticity in terms of CFD results interpretations. However, during this period, vortex methods experienced some difficulties. First they had to face an inaccurate modeling of the viscous effects due to the particle distortion phenomenon. They were also confronted with the expensive resolution of the Poisson equation allowing to recover the velocity field from the vorticity. Finally, they failed in providing both accurate, efficient and generalized treatment of boundary conditions on solid walls due to the velocity formulation of the no-slip boundary conditions. For all these reasons, and despite the efforts made to design solutions to overcome the difficulties, these solutions remained partial and the purely Lagrangian vortex methods did not manage to become a common tool for the simulation of incompressible real and complex flows. In order to bypass their intrinsic limitations, they started to move towards hybrid formulations in the early 1990s, resulting in the introduction of underlying grids with remeshing. The presence of Cartesian grids indeed enables the use of Fast Poisson solvers as well as grid-based schemes for the discretization of the diffusion term, and makes possible the handling of immersed boundary methods for the treatment of the solid bodies, while ensuring a permanent control of the particle overlapping through regular redistribution operations on the grid. In the meantime, the particle discretization of the flow and the Lagrangian advection scheme allow us to keep most of the advantages of the Lagrangian methods. The introduction of grid and remeshing process in vortex methods marks the birth of hybrid methods, with among them the Eulerian-Lagrangian decomposition domain methods and the Vortex-Particle-Mesh methods (like the Vortex-In-Cell and the semi-Lagrangian vortex methods). All these advances allowed by the hybrid approaches increased the interest and potential of VM in handling challenging CFD problems. Important efforts were therefore supplied to design efficient algorithms in the context of hybrid vortex methods and eventually gave birth to high performance scientific libraries, able to perform simulations of complex real flows arising in diverse application fields like hydrodynamics, turbulent wake dynamics, scalars transport or porous flows.

Future Research Perspectives
As described before, the academic research on vortex methods has brought multiple achievements on the accuracy, stability and applicability of this family of approaches. They have also been boosted with powerful high performance computing tools increasing their capability to handle complex problems. Several perspectives appear to fill the future of the research related to this family of methods. The most current issue would be to strengthen their abilities for realistic and industrial applications. Below, some research perspectives are listed.
• Turbulent and high Reynolds number flow simulations. The recent researches highlight the road map for this target with two possibilities. The first option consists in proposing modern turbulence models for the velocity-vorticity filtered Navier-Stokes equations (i.e., LES-like approaches), based for example on those already developed by Mansfield et al. [148] on one side or Cocle et al. [151] on the other side. The second view, that takes its idea from the inherent nature of vortex methods, is the multiresolution approach (based on the pioneer works of Cottet and Wray [147,157]) which consists in using two different levels of resolution for vorticity (fine level) and the velocity (coarse level). This multi-resolution would have a deep interest in problems where important Schmidt numbers are involved like in heat transfer (e.g., thermal conduction and convection [158]), passive-scalar transport [153] or sedimentation flows [141]. • Fluid Structure Interactions (FSI). The handling of no-slip boundary conditions has long been one of the main difficulties of vortex methods. The implementation of the Brinkman penalization method has significantly improved this drawback in the context of Vortex-Particle-Mesh methods. However, this approach is limited because of its low order and accuracy. The very recent introduction of immersed interface method (IIM) [131,132] and immersed lifting line approach (ILL) [12,99] introduced a real improvement and allows us to start a real progress in the challenging field of FSI. • Environmental and industrial particular flows. Vortex methods belong to the family of particle methods with direct treatment of particles. Therefore, their application to problems like sediment transport [141,159] or internal multi-phase flows in industry with particle mixing, dispersion, deposition, and particle-wall collision [160], seems natural and deserves a special place in the future perspectives. • Multi-phase flows. Besides the industrial fluid/particles two-phase flows evoked above, vortex methods also have a role to play in the handling of multi-phase flows in general, with for instance the development of models able to capture interfaces between two phases and to handle surface tension, variable fluid properties and high mass density-ratio flows. To the authors knowledge, only few works based on vortex methods have been dedicated to this topic [161]. • Finally, future vortex computations need to be continually strengthened with novel high performance computing parallelism techniques (hybrid multi-CPU/multi-GPU) and powerful algorithms (like the very recent 2D-3D Poisson solver released in 2021 by Caprace et al. [162]) in order to achieve more complex flow approximations. Funding: This research received no external funding.

Conflicts of Interest:
The authors declare no conflict of interest.

Abbreviations
Here is the list of the abbreviations used in this manuscript, in order of appearance: particle v p volume of particle p (denoted v p in 2D) x p location of particle p W ε denomination of the smoothing radial symmetric function in Lagrangian methods ε radius of the smoothing radial symmetric function Γ circulation × cross product ω vorticity field (denoted ω in 2D) ω 0 initial vorticity field (denoted ω 0 in 2D) (·) h quantity approximation u ∞ far-field velocity α p local circulation around x p (denoted α p in 2D) convolution product