Numerical simulation of 2D-vorticity dynamics using particle methods

Numerical simulation of 2D-vorticity dynamics using particle methods

Computers and Mathematics with Applications 69 (2015) 1484–1503 Contents lists available at ScienceDirect Computers and Mathematics with Application...

3MB Sizes 0 Downloads 17 Views

Computers and Mathematics with Applications 69 (2015) 1484–1503

Contents lists available at ScienceDirect

Computers and Mathematics with Applications journal homepage:

Numerical simulation of 2D-vorticity dynamics using particle methods E. Rossi a , A. Colagrossi b,∗ , G. Graziani c a

DM, University of Rome, Sapienza, Piazzale Aldo Moro 5, 00185, Roma, Italy


CNR-INSEAN, Marine Technology Research Institute, Rome, Italy


DIMA, University of Rome, Sapienza, via Eudossiana 18, 00184, Roma, Italy



Article history: Received 27 June 2014 Received in revised form 4 April 2015 Accepted 7 April 2015 Available online 9 May 2015 Keywords: Vorticity dynamics Vortex method Particle methods Meshless methods Smoothed Particle Hydrodynamics method Diffused Vortex Hydrodynamics

abstract In this paper two dimensional free vorticity dynamics is studied using two different particle methods. The first one is a Diffused Vortex Hydrodynamics (DVH) and the second one is a Smoothed Particle Hydrodynamics (SPH) method. These two methods present some similarities linked to their meshless nature but they are based on different numerical approaches. In this work advantages and drawbacks are highlighted by testing the particle methods on selected test-cases and by performing heuristic convergence measurements. The DVH method discussed in this work is characterized by the use of a regular distribution of points to perform the vorticity diffusion process. This redistribution avoids excessive clustering or rarefaction of the vortex particles providing robustness and high accuracy of the method. © 2015 Elsevier Ltd. All rights reserved.

1. Introduction This work concerns modelling viscous incompressible flow by using two particle methods: a Diffused Vortex Hydrodynamics (DVH) and Smoothed Particle Hydrodynamics (SPH). An extension of the previous work [1] is performed by analysing more deeply these two numerical schemes in terms of accuracy, CPU costs, rate of convergence, and range of applicability. Particle methods have been largely developed in the last years since they present interesting features which allow for the simulation of engineering problems involving complex interfaces dynamics such as vortex sheets and free-surface deformation (for a more general discussion on particle methods see e.g. [2]). The vortex method has several advantages. The first one is directly linked to the vorticity formulation of the incompressible Navier–Stokes equations for which the pressure field is no more a direct unknown of the problem. Consequently it is possible to discretize the fluid domain by means of vortex particles which are present in the rotational region only. Therefore the computational domain can be in generally much smaller in comparison with other methods. These advantages are no more true when considering more general flows, like for example compressible flow, for which a full discretization of the domain is needed to solve the evolution of the density field (see e.g. [3]). A second advantage is linked to the Lagrangian description. Indeed, the vortex particles can be moved in the domain in a Lagrangian way, and this handling of convection reduces the numerical dissipation which has always been inherent to mesh-based approaches for the nonlinear term of the Navier–Stokes equations. The vortex method is also well-known as a CPU cost effective method

Corresponding author. Tel.: +39 06 50 299 343; fax: +39 06 50 70 619. E-mail address: [email protected] (A. Colagrossi). 0898-1221/© 2015 Elsevier Ltd. All rights reserved.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


because of the less restrictive time step constraint (see [4]). Further, the boundary conditions at infinity are automatically satisfied and this is one of the advantages with respect to the other particles methods like SPH. This Lagrangian method has an essentially mesh-free nature even though ‘‘redistribution techniques’’ are needed in order to preserve the accuracy of the method (see e.g. [5]). Some of the above mentioned features attracted many researchers to follow this approach since Chorin applied it to slightly viscous flow [6]. The vortex method is classified by some aspects such as the element type for the discretization, the treatment of the diffusion term and the method of velocity field evaluation. The method proposed by Chorin in [6] was a so called Random Vortex Method (RVM) in which diffusion is performed giving a random displacement, with zero mean, to the vortex particles. More recent advances on RVM are shown in [7] where an efficient algorithm is developed for diffusion in the presence of complex two-dimensional solid boundaries. In particular the numerical method proposed in this work is based on the representation of point vortex particles where the velocity field is evaluated through the Biot–Savart law allowing for a total mesh-free representation. From this point of view the DVH is similar to other methods already described in the literature as the Viscous Vortex Particle Method (VVPM) presented in [8]. The main novelty in the DVH scheme regards the solution of the diffusion process and its use to regularize the particles spatial distribution. Indeed, the diffusion process is performed through superimposition of elementary ‘‘Lamb–Oseen’’ solutions (deterministic diffusion algorithm, Benson et al. [9]). In the adopted diffusion scheme each particle gives its diffusive contribution to the circulation of a new set of particles located on regular distribution of nodes. This algorithm avoids the excessive clustering or rarefaction of the vortex particles guaranteeing robustness and convergence property of the scheme. The procedure described above avoids using a remeshing procedure where a specific interpolation function and a frequency of remeshing in time need to be selected (see e.g. [8]). The Smoothed Particle Hydrodynamics method was introduced about three decades ago in astrophysics by Lucy [10] and Monaghan [11]. For more than ten years it was solely used for astronomical problems. Then, its applications spread to several physics and engineering problems, ranging from solid mechanics to multiphase flows. Monaghan [12] firstly applied it for simulating free-surface inviscid flows. With respect to Eulerian finite difference method, it was found to be more suitable to handle breaking wave and fragmentation processes, due to its Lagrangian nature, and furthermore, within this method, the free-surface boundary condition is treated in a straightforward fashion (see e.g. [13]) which is an advantages also with respect to other Particle Methods including the Vortex Method. Thanks to this valuable features, interest in SPH has been growing very fast in the last decade and several researchers have devoted their activity to improve this method. We underline that the present work is not aimed at comparing DVH and SPH in terms of abilities in solving specific fluid-dynamic problems, but more at showing the similarities of these two different particle methods. Indeed, the analysis performed shows similar behaviour on the convergence of the two schemes as theoretically predicted in [1]. The paper is organized as follows: Sections 2 and 3 are dedicated to a brief introduction of the DVH and SPH models. Section 4 deals with the simulation of a Lamb–Oseen vortex, using three different Reynolds numbers, and with the simulation of the merging process of a pair of co-rotating vortices. The abilities of DVH and SPH to recover accurate and convergent solutions are commented in detail. Although this work simply focuses on free vorticity dynamics, the last section is dedicated to the description of potential advantages in using the DVH also for simulating viscous problems around body of complex geometry. 2. Diffused Vortex Hydrodynamics Considering a flow domain Ω the Navier–Stokes equation for a 2D incompressible flow takes the following form in terms of vorticity, ω:

(∂t + u · ∇)ω(x, t ) = ν 1ω(x, t ),

∀x ∈ Ω


being u the velocity field and ν the fluid kinematic viscosity. Eq. (1) describes simultaneously both the advection and the diffusion of the vorticity field. According to the well known operator splitting scheme [6] these two steps are separately considered in the computational procedure. The approximate solution of the governing equation is obtained by the sequential solution of an Eulerian (inviscid) step and a purely diffusive one. In the Euler step an inviscid fluid is considered where the vorticity field is transported according to:

(∂t + u · ∇)ω(x, t ) = 0.


In the diffusive step the fluid is considered at rest and the pure diffusion of the vorticity is obtained by solving the heat equation:

∂t ω(x, t ) = ν ∆ω(x, t ).


After the splitting of the governing equation, it is clear that an appropriate procedure for both diffusion and advection of the vorticity must be devised. The natural choice for the numerical solution scheme appears now as viscous vortex method which is based on the original papers by Chorin [6,14], together with a deterministic diffusion algorithm, Benson et al. [9].


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

2.1. Particle discretization of the vorticity field In the framework of vortex method the vorticity field ω should be considered a measure, i.e. using the Lebesgue integral

ω(f ) =

f (x)ω{dx}


where f is a regular test function. In this case discretizing the vorticity field with N point vortices, each with circulation Γi and position xi , means to use the following approximation for ω N  

ω(f ) =

f (x) Γi δxi {dx},



where the vorticity field is discretized. Regarding δxi {dx}, this is the two dimensional Dirac delta measure

f (x)δx0 {dx} = f (x0 ).


It is possible to approximate the Dirac delta measure using a positive smooth function δϵ (x − x0 ) which converges in the weak sense to the Dirac delta measure, i.e.

 lim


f (x)δϵ (x − x0 ) dx = f (x0 ).


An example of such a function is:

  |x − x0 |2 δϵ (x − x0 ) = exp − 4π ϵ 2 4ϵ 2 1


and finally the vorticity field can be numerically expressed as:

ω(x, t ) =

N 

Γj (t )δε (x − xj ).


j =1

The evolution of the vorticity field is now computed through the advection and the diffusion of these point vortices. The adoption of the vorticity as the dependent variable avoids the discretization of the irrotational regions while allows for reserving the computational resources to the vortical zones. The adoption of an integral representation of the velocity field yields the exact enforcement of the far field boundary conditions. One of the drawbacks of the purely Lagrangian vortex method is the excessive clustering or rarefaction of the vortex particles (see e.g. [15]). The usual way to solve this problem is to introduce a redistribution procedure. In our case the adopted diffusion procedure avoids this requirement because a regular distribution of points is used, only for this step, and each particle gives its diffusive contribution to the new particles located on those points. At the end of the diffusive step the former set of vortex particles is replaced by the new vortices defined on the computational stencil. 2.2. Advection The inviscid evolution of the vortex particles is described by the Euler equation which states the material behaviour of the transported quantity. Therefore, each vortex carries around its circulation without diffusion. In order to evaluate the advection step, the velocity field must be computed to obtain the new position of the vortex particles. The velocity field can be calculated using:

    ω(y, t ) K(x, y)dy + ∇ Φ (x, t ) u(x, t ) = Ω

(x − y) 1   K(x, y) := ∇ ⊥ G(x, y) = − e3 × 2π |x − y|2


where ∇ ⊥ indicates the skew gradient operator, G is the free space Green function of the Laplace equation, in a 2D framework: G(x, y) = −1/(2π ) ln |x − y| and e3 the normal unit vector orthogonal to the fluid domain plain. The second component, ∇ Φ , is a correction to the velocity caused by boundary surfaces in the flow. Computation of this term typically requires the use of boundary integral methods, however, in this work only unbounded problem is considered and the potential flow term is always null.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Substituting (9) in (10) we obtain the expression for the velocity field for every single point vortex: u(xi , t ) =

N 

Γj (t )Kε (xi , xj ),


j =1

where Kε is a regularized version of the kernel K (see e.g. [16]). The operation count for the computation of the velocity of each particle is evidently given by N 2 if Eq. (11) is directly used. In order to achieve a faster solution the flow domain is recursively divided and a multipole expansion for the kernel K is adopted with a computational cost of order N log N which is comparable with the performance of other efficient numerical techniques (for details see [17]). Once the velocity field has been evaluated through (11) each vortex is displaced by the fluid velocity according to: d dt

xi (t ) = u(xi , t ).


The numerical time integration is performed using a fourth order Runge–Kutta algorithm. 2.3. Diffusion Viscous effects are accounted for in the diffusive step where the fluid velocity is neglected and the diffusion of the field vortices is evaluated by solving Eq. (3). The exact solution of the heat equation with initial datum a single point vortex at time t0 with initial circulation Γ0 and position x0 is |x−x0 | Γ0 − e 4ν(t −t0 ) . 4π ν(t − t0 ) 2

ω(x, t ) =


Therefore, even starting with all the circulation concentrated at time t0 on a single point vortex, for t0 + 1t we have a distribution of vorticity on all R2 , leading to a situation that is not possible to deal with numerically. To solve this problem we define a finite support for the Lamb–Oseen vortex rewriting Eq. (13) as: |x−x0 |2 Γ0 e− 4ν 1t , 4π ν 1t

 

ω(x, t0 + 1t ) =

ω(x, t0 + 1t ) = 0

if |x − x0 | ≤ Rd



where Rd is the radius of the truncated support. Using this approximation we obviously introduce an error ξ in the conservation of the circulation:

   Γ (t0 + 1t ) − Γ0   ξ =   Γ



where Γ (t0 + 1t ) is the circulation contained at time t0 + 1t in a circle of radius Rd . A simple integration of (13) gives for Γ (t ) the following expression:

Γ (t ) = Γ0 1 − exp −



4ν(t − t0 )


and by substituting this relation into (15) we obtain:

  R2 ξ = exp − d . 4ν 1t


Expression (17) for the relative circulation error will be used in the next section to determine the diffusive time step. Following the deterministic diffusion algorithm of Benson et al. [9], we can now use Eq. (14) to distribute the circulation of each vortex over a regular lattice characterized by a uniform spacing 1x. The vorticity field after diffusion is then lumped into point vortices located at the nodes of the lattice. The use of a Cartesian lattice with a uniform spacing is not real constraint for the proposed DVH model, indeed the diffusion procedure described above can be generalized on non uniform point distribution. However, in this work, since we want to focus the attention on the convergence of the DVH solutions, we preferred using of a simple uniform spacing point distribution. To preserve the conservation of the circulation of each vortex we must undergo a renormalization procedure of the diffused vorticity. Call xij the position of a generic node on the lattice, the contribution to the circulation of the node xij of the vortex located in xk with circulation Γk , in the time interval 1t using Eq. (13) and Stokes’ theorem, is: (k)



|xij −xk |2 Γk e− 4ν 1t 1x2 . 4π ν 1t



E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Thus to normalize this circulation in xij the following formula can be used: (k)′


Γk =  (k) Γij(k) Γrl



where the sum over r and l is made over set of regular points with distance from xk less than or equal to Rd . It is easy to see that a sum over i and j of the normalized circulations gives exactly the circulation Γk . After all the vortices have spread their diffusive contributions to the circulations of the points, to reduce the number of vortices generated in this diffusive process the nodes with circulation lower than a fixed cut-off Γcut −off can be eventually removed. 2.4. Choice of the time steps In the following we consider dimensionless variables, where U and L are the reference velocity and length of the problem. With this choice we also define the Reynolds number as Re = UL/ν . The spatial resolution linked to the diffusive Cartesian lattice is indicated in the following by the ratio L/1x. Because advection and diffusion of the point vortices are treated separately, we must deal with two different time steps: one for advection 1ta and one for diffusion 1td . The diffusion time step 1td depends on the Re number. First we consider that during a time step 1td a single vortex diffuses within a finite distance Rd . This means that a new Gaussian vorticity distribution is generated by this single vortex and it is further discretized on the Cartesian lattice during the diffusion step. The number of points inside the radius Rd for a lattice characterized by a uniform spacing 1x is: Nnode

  2  Rd , = π 1x


being ⌊x⌋ the largest integer not greater than x. It is important to underline that, Nnode (or equivalently the ratio Rd /1x) together with 1x are the key parameters of the DVH model, because they fix the discretization level used in this particle method. A typical number for Nnode in 2D is 51 which corresponds to a ratio Rd /1x equal to 4. Using expression (17) of the previous subsection the following relations hold:

 ξ = exp −


4 ν 1td

⇒ 1td


Rd /1x

 =



4 ln(1/ξ ) (L/1x)2



which allows for the evaluation of the diffusive time step from:

• the adopted spatial resolution L/1x • the ratio Rd /1x and • the Reynolds number, once the error ξ has been fixed (in particular we have taken ξ = 10−5 ). Conversely, the advection step 1ta can be chosen by considering the flow velocity U and the discretization 1x used during the diffusion step. To avoid that the vortex particle distribution becomes too irregular, a first evaluation of the advection time step can be chosen as:

1ta = Co

1x U

⇒ 1ta


= Co




where Co is a Courant number parameter of order of one. It is important to underline that for vortex method schemes based on operator splitting, as the present DVH, the stability constraint on Co is quite less restrictive with respect to other numerical methods (see [4]). Once the advection time step is evaluated it needs to be rearranged in order to synchronize diffusion and advection, therefore, 1td and 1ta need to be in integer ratio, i.e.

 N1t =

1td 1ta

⇒ 1tanew = 1td /N1t .


The above expression is used only if N1t ≥ 2. In the case N1t ≤ 1, the minimum between 1ta and 1td is chosen for both the two time steps (this is typical condition when simulating low Reynolds numbers flow where 1td becomes similar or even smaller than 1ta ). Thanks to the double constraint on the two time-steps, rearranging the above Eqs. (21)–(23) we can link the spatial resolution to the Reynolds number through: L


 =

1 Co N1t

(Rd /1x)2 4 ln(1/ξ )



E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


When the ratio Rd /1x is fixed, the spatial resolution increases linearly with the Reynolds number. Using N1t greater than one, allows for the reduction of the spatial resolution needed for a given Reynolds number. We now consider the so called Reynolds-cell number defined as Re1x =

U 1x


= Re





To capture the smallest vorticity scale for a generic flow problem, the Reynolds-cell should be of order of one [18]. As underlined in [19], the limit on the Reynolds-cell number for the mesh-based solvers has to be read only in terms of accuracy of the solution and not in terms of stability of the numerical scheme. Recently, a discussion regarding the link between the spatial resolution, Mach and Reynolds numbers in particle method context has been provided by [20]. Substituting Eq. (25) in Eq. (24) we get: Re1x = N1t Co

4 ln(1/ξ )

(Rd /1x)2



Once the ratio Rd /1x is set to a typical value of 4, ξ = 10−5 , and Co ∼ 1, the right factor is almost equal to 3. Therefore (26) implies that the value of N1t needs to be kept as small as possible in order to maintain Re1x in the right range for studying properly 2D vortex dynamics. 2.5. Evaluation of pressure field and kinetic energy by the DVH method As already commented in the introduction, in the vorticity formulation of the incompressible Navier–Stokes equations the pressure field is no more a direct unknown of the problem. Therefore, in order to evaluate the pressure field from the DVH solution, a Poisson problem for the pressure needs to be solved. In the context of bi-dimensional free vorticity dynamics it is possible to write a solution for a fluid contained in a ball BR of radius R, in the following way (for details see [21]):

  1 2  |u| − |u∞ |2 P := p − p∞ +   2        ∂(u(y) · n) − ω(y) u(y) · K(x, y) dy + P (x) = − + ν ∇ω(y) · τ G(x, y) dl  ∂t ∂ BR   BR     − P (y) τ(y) · K(x, y) dl,


∂ BR

where n and τ are the normal and tangent vectors on ∂ BR respectively. Taking the limit R → ∞ the integrals on ∂ BR go to zero and using a discrete vorticity field with N vortices, each with circulation Γi and position xi , the field P can be numerically evaluated as: P (x) = −

N 

Γi u(xi ) · Kϵ (x, xi ).



Regarding the DVH evaluation of the kinetic energy for the fluid contained in a ball BR of radius R, this global quantity can be calculated using the relation: Ek (R) =

Γ2 1 log R + 4π 2

ψωdl + O BR

  1



being ψ the stream function and Γ the total circulation in BR . For the test cases studied in this work, choosing a radius R large enough, the third term on the right-hand side of (29) can be neglected with respect to the other two terms. Formula (29) allows using only the vorticity within BR avoiding the calculation of the velocity field in the whole ball BR . The second term in Eq. (29) is known in the literature as the ‘‘excess energy’’ (see e.g. [21]). 3. The smoothed particle hydrodynamics model In the present section a SPH scheme is briefly described. The fluid is assumed to be barotropic and weakly-compressible and the reference equations are the Navier–Stokes equations. As discussed in [22,23] different state equations, p = f (ρ), can be used in the SPH scheme to model weakly-compressible fluids. Here, a simple linear state equation is used to match the pressure and density field: p = c02 ρ − ρ0


where, c0 , ρ0 and p0 are the sound speed, the fluid at rest density and pressure, respectively. The speed of sound c0 is set in order to guarantee density variation smaller than 0.01ρ0 . This is ensured through the following inequality:

c0 ≥ 10 max max t

p/ρ, max |u| .



E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Given a set of particles each characterized by its own mass mi , the particle densities can be expressed through the distribution:

ρi =

mj δε (|xi − xj |)



where δε = δε (|xj − xi |) is a positive, smooth approximation of the δ function as already introduced for the DVH (see Section 2.1). In this work a C2 Wendland kernel (see e.g. [24]) has been used. δε has a compact support of radius 2ε , where in the SPH literature ε is referred to as the smoothing length. By time differentiation of Eq. (32) we get the continuity equation in the SPH formalism:

ρ˙ i =

mj (ui − uj ) · ∇i δε (|xi − xj |).



Symbol ∇i indicates the differentiation with respect to the position of the ith particle, For the sake of simplicity in what follows ∇i δε (|xi − xj |) will be shortened to ∇i δε . From the numerical point of view, evaluation of the density through time integration of (33) is preferable with respect to the direct use of (32) and becomes crucial when dealing with free-surface flows (see e.g. [13]). For an isentropic fluid the sum of the time variation of kinetic energy Ek and internal energy Ei of the particles system reads as:

E˙k + E˙i =

˙i + mi ui · u






ρ˙ i = 0.


By substituting Eq. (33) in (34) and by using the kernel property

∇i δε = −∇j δε the acceleration of the ith particle is given by:

˙i = − u

 mj





 ∇i δε . 2




The acceleration given in Eq. (35) is due to the pressure forces, while the viscous one can be modelled through the viscous formula of Monaghan and Gingold [25] that preserves both linear and angular momenta. Using the above equations for ρ˙ i and u˙i the SPH scheme reads:

ρ˙ (t ) = −  m (u − u ) · ∇ δ j j i i ε i    j        pj µ  pi   u˙ i (t ) = − + 2 ∇i δε + 2 mj πij ∇i δε mj 2 ρi ρj ρ0 j j      pi  µ πij   ˙ (uj − ui ) · ∇i δε m − e ( t ) = −  j i   ρi2 ρ02 2  j  x˙ i (t ) = ui (t )


where ei is the specific internal energy of the ith particle. Symbol µ denotes the fluid dynamic viscosity while πij is the Monaghan and Gingold formula to model viscous forces:

πij = 2 (n + 2)

(uj − ui ) · (xj − xi ) ∥xj − xi ∥2


where n is the spatial dimension of the problem at hand. Thanks to the use of Eq. (34) the exact momentum and energy conservations in (36) are guaranteed regardless of the state equation adopted. The fluid particles are initially placed using the algorithm described in [26]. Thanks to this procedure, at the initial instant all particles have approximately the same volume, namely V0 , which is equal to the fluid domain volume divided by the 1/n number of fluid particles. Consistently, the particle mean spacing is denoted by 1x = V0 . The average number of particles in the kernel support is set by choosing the ratio ε/1x. In the present work ε/1x is set equal to 2 which in two dimensions corresponds to about 50 interacting particles. For ε → 0 and 1x/ε → 0 the system (36) converges to the Navier–Stokes equations (see e.g. [27,28]). Along with the volume distribution, the initial pressure and velocity fields are prescribed as well. The initial density distribution ρi (t0 ) is evaluated by means of the state equation and the particle masses are computed through the equation mi = V0 ρi (t0 ). The mass of the ith particle remains constant during the time evolution ensuring the total mass conservation of the particles system.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


4. Numerical validation 4.1. Lamb–Oseen vortex test case In this section an infinite fluid subjected to a radial force field: fr = −β 2 r rˆ is studied. β is a constant parameter, r is the radial coordinate and rˆ is the radial unit vector. A velocity field merely tangential to circumferences centred in r = 0 is considered. An incompressible isotropic solution is searched and under these hypotheses the Navier–Stokes equations reduce to:

 2 1 ∂p u   − β 2r , − = − r ρ0 ∂ r      1 ∂ ∂u u ∂u   =ν r − 2 . ∂t r ∂r ∂r r


where u is the tangential velocity while the radial one is identically null. In this framework the two equations are decoupled and it is possible to integrate u by the second equation and to evaluate the pressure field by the first one. Furthermore, we underline that the evaluation of u is not affected at all by the radial force field which is involved only in the solution of the pressure field only. Therefore, the DVH solutions are not affected by the parameter β , while this is not the case for the SPH method; we will come back on this point at the end of this subsection. The Lamb–Oseen vortex is characterized by an initial velocity field equal to: u0 (r ) = ω0


1 − e−(r /L)





where L is the length scale of the problem. The corresponding vorticity contribution is a Gaussian distribution:

ω(r , 0) = ω0 e−(r /L) . 2


Therefore ω0 is the vorticity in r = 0 at t = 0. The exact solution of Eq. (38) for the initial datum is:

   f (t ) :=



, g (r , t ) := f (t ) 2 4π t ω0 + Re L    L2  uexact (r , t ) = ω0 1 − e − g ( r ,t ) ,



(see Fig. 1(b)) where the Reynolds number is defined as Re = π ω0 L2 /ν . The vorticity is then:

ωexact (r , t ) = ω0 f (t ) e−g (r ,t ) ,


a Gaussian vorticity shape is therefore maintained for all the time evolutions (see Fig. 1(a)). Considering a null radial force field, β = 0, the pressure evaluated with the first equation of (38) is: pexact (r , t ) − p∞

ρ0 ω02 L2


f (t ) 8

[2e−g (r ,t ) − e−2g (r ,t ) − 1] 2 log(2) + + 2[Ei(1, 2g (r , t )) − Ei(1, g (r , t ))] g (r , t )

f (t ) log(2) 4



being Ei the exponential integral function, p∞ the pressure at infinity and the functions g (r , t ) and f (t ) defined in Eq. (41). As expected, in the core of the vortex the pressure drops and the negative pressure peak reduces as time increases (see Fig. 1(c)). Such a simple problem is used to measure heuristically the convergence of DVH and SPH varying the Reynolds number. Indeed, the exact solution can be used to evaluate the relative errors over vorticity and velocity during the simulations. In particular the relative error Ef for the function f (x, t ) is defined using a combination of L1 norm in space and L∞ in time:

 Ef = max

t ∈(0,tend ]

|fexact (x, t ) − f (x, t )|dx  , |f (x, t )|dx Ω exact


where Ω is the computational fluid domain. The latter is defined as a circular domain with radius R = 20L, in this way Ω is large enough to ensure that the vorticity field is almost negligible on ∂ Ω . The final time, tend , is set equal to Re/ω0 , in this time range the maximum vorticity in Ω is reduced to about 7% of the initial value ω0 .


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 1. Vorticity (top), Velocity (middle) and pressure (bottom) profiles for the Lamb–Oseen Vortex at t = 0 (solid lines) and t = Re/ω0 (dashed lines).

Fig. 2. Pressure profiles for the Lamb–Oseen Vortex at t = 0 using β = 0, 0.05ω0 (left) and β = 0.10ω0 (right). For β ̸= 0 the pressure pR := p(r = R) is set to zero while for β = 0 a positive constant p∞ is considered.

In the analytical solution Re simply changes the time scale, however, in the next sections it is shown that increasing the Reynolds number the advection step becomes more and more important introducing errors both in the DVH as well as in the SPH numerical scheme. Therefore, different error levels are expected upon changing the fluid viscosity. Since in the DVH the boundary condition at infinity is automatically satisfied, there is no need to confine the fluid domain Ω . The vorticity field is discretized only for r ≤ R, with R = 20L while beyond this distance the vorticity is neglected because of its small values. Conversely, this is not the case for the SPH method which requires a closed domain. This is realized by considering an annular region r ∈ (R, R + 2ε) where fictitious equispaced particles are positioned and where values for u and p are prescribed. This problem can be demanding for SPH method since it is mainly dominated by the viscous forces (see [29]). Since SPH is based on low order integral interpolation formula (see e.g. [30]) the error on second order operators is quite sensitive to the particle disorder. To limit this error, the number of particles inside the support of δε , needs to be large enough (see [30,24,31]). In order to avoid the excitation of the so called tensile instability [32,33] in the SPH method, a positive pressure field needs to be set in the initial condition. This can be done in two ways (i) using Eq. (43) and setting p∞ ≥ ρ0 ω02 L2 log(2)/4, (ii) using a radial force field with a large enough intensity, β , and setting pR := p(r = R) ≥ 0. The second strategy, with pR = 0, has been used in [34] and in [28] since in the problems studied the fluid domain was confined by a free surface. Fig. 2 shows the pressure profile for different β values. Increasing enough this parameter the parabolic profile dominates the pressure law (43). Even if the values of β and p∞ are not affecting the velocity solution of Eq. (38), in Sections 4.1.2 and 4.1.3 it is shown that the SPH numerical errors can be largely modified by these two parameters.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Fig. 3. Left: Relative maximum error over vorticity as a function of L/Rd for various ratios Rd /1x and ξ = 10−5 . Right: Relative maximum error over vorticity as a function of L/Rd for various ξ and fixed ratio Rd /1x.

Fig. 4. Maximum relative error over vorticity (4(a)) and over velocity (4(b)) as functions of L/Rd for various Reynolds numbers and fixed ratio Rd /1x.

4.1.1. Numerical results for the DVH In this section the DVH model is tested on the Lamb–Oseen vortex problem. The test has been made in two steps: a set of simulations without advection has been performed to better control the error coming from the diffusive step. This test corresponds to the limit Re → 0, then a new set of simulations has been made taking advection into account and considering increasing Re numbers. As discussed in Section 2.3, the diffusion of each point vortex is evaluated by discretizing Eq. (14) on a regular distribution of points. Thus to reduce the errors made in this process we can either increase the number of points contained in the support of radius Rd (i.e. increase the ratio Rd /1x) or decrease the fixed relative error ξ (see Eq. (15)) for the diffusion of each vortex. This second choice implies a reduction of the diffusive time step (see Eq. (21)), which means that the diffusion equation is resolved with higher accuracy in time, while increasing the ratio Rd /1x implies higher accuracy in space. We call Eω and Eu the relative error over vorticity and velocity, respectively, measured according to formula (44). Fig. 3(a) shows the relative errors Eω obtained for the Lamb–Oseen problem for the case Re → 0, changing the spatial resolution L/Rd . Four different ratios Rd /1x are considered: 3, 4, 5 and 6 and ξ is set equal to 10−5 . From this plot it is possible to see that there are essentially two regimes varying the spatial resolution L/Rd : (i) a first regime, for low spatial resolutions, where the errors reduce (i.e. the numerical solver is convergent) and (ii) a second regime where Eω saturates. The convergence steepness in the first regime depends on the ratio Rd /1x. In particular a convergence rate equal to 5 for Rd /1x = 3 up to 10 for Rd /1x = 6 is obtained. On the second regime the solution is no more affected by the number of vortices used to discretize the fluid domain; the only way to reduce the errors is to increase the ratio Rd /1x or to reduce ξ . This kind of convergence behaviour is typical of the meshless methods where the convergence requires the reduction of some parameters, in this case 1x/L, the ratio 1x/Rd and ξ (see e.g. [35]). In the next section we show that also the convergence of the SPH method depends on a similar requirement. For the highest ratios Rd /1x = 6 the saturation level of about 10−4 is reached very soon and the error EΩ is mainly dominated by the factor ξ . Fig. 3(b) shows Eω for three different values of ξ : 10−3 , 10−5 and 10−7 with fixed Rd /1x = 5. It is straightforward to note that the saturation level decreases with ξ so that we can expect a convergence to the exact solution if we simultaneously decrease 1x/Rd and ξ . Fig. 4(a) and (b) show the relative errors Eω and Eu for Re numbers: 10, 100, 1000 with fixed Rd /1x = 4 and ξ = 10−5 . For Eω the case without advection (Re ≪ 1) is also shown for comparison. As expected Eω increases with Re, due to the errors linked to the advection steps, in particular for Re equal 10 the errors are practically the same as the case Re ≪ 1


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 5. Maximum relative error over vorticity (5(a)) and over velocity (5(b)) as functions of L/Rd for various ratio Rd /1x and fixed Reynolds numbers.

(convergence rate equal about to 8) while for Re equal to 1000, the saturation of Eω is no more visible in the adopted spatial resolution range and the convergence rate reduces to about 4. Furthermore, for Re = 1000, the L/1x range covered by the first regime becomes quite wider with respect to the low Reynolds number cases, this means that a large number of vortices is now needed (about 105 ) to reach an error level Eω < 10−3 . Regarding the errors on the velocity field, Eu , these are higher than Eω , however Eu shows an almost constant convergence ratio of about 1.8 for all the three Re numbers, without any reduction due to saturation error level. Using Eq. (24) a link between the ratio N1t and Reynolds number can be derived:

 N1t =

 1 (Rd /1x)2 Re, 4 Co ln(1/ξ ) L/1x


therefore, upon fixing Rd /1x and the resolution L/1x, we must perform, for every diffusive step, N1t advection steps cumulating errors as N1t increases. Indeed, if N1t is large enough, as in the case Re = 1000, where N1t is equal to 30 for the lowest resolution, irregularities in the vortex particles distribution become quite visible. Fig. 5(a) and (b) show Eω and Eu for three different values of Rd /1x: 3, 4 and 5, fixed Re = 100 and ξ = 10−5 . As previously shown for the case Re ≪ 1, convergence with Rd /1x is obtained also here in the presence of advection steps. For the lowest ratio Rd /1x = 3 saturation is visible for both Eω and Eu . Summarizing, given the Reynolds number, accurate results require the right combination between spatial resolution L/1x and Rd /1x (linked to the number of nodes inside the diffusive support). A good combination can be obtained through Eq. (45) by reducing N1t as much as possible. 4.1.2. Numerical results for SPH using a radial force field In this section the SPH model is tested on the Lamb–Oseen vortex problem. To avoid tensile instability the strategy adopted in [34] and in [28] is used adding a radial force field (see Eq. (38)) and setting the parameter β equal to 0.1ω0 . The pressure, pR , at the domain boundary (i.e. r = R) equal to zero. We underline that, in the SPH model, the smoothed length 2ε is equivalent to the radius Rd used in the DVH. Therefore to check the accuracy of the numerical solution as well as the convergence of the SPH scheme, different tests changing the spatial resolution L/(2ε) and the ratio (2ε)/1x are performed. Fig. 6(a) shows the relative errors Eu (see Eq. (44)) on the velocity field for the smallest Reynolds number (Re = 10). For this case a relative error of order 10−2 is measured using 2ε/1x = 4 and 8 while, using only 9 interacting particle neighbours (i.e. 2ε/1x = 2), the error is 10−1 . Further, a small convergence rate is always obtained and a saturation level is quite visible for the two highest ratio 2ε/1x. The results are in line with the analysis performed in [30] indeed, as it is well known in the SPH literature, the method converges when both 1x as well as 1x/ε simultaneously go to zero (see e.g. [27]). However, the results show also that a reduction of both 1x and 1x/ε is not enough to get convergence to the analytical solution. This can be related to the shape of the kernel δε , as shown for the DVH; however this aspect has not been investigated in this work. With respect to DVH, considering the same spatial resolution L/1x and the same Nnode (i.e. Rd = 2ε ), the SPH errors on Eu are almost two order of magnitude larger. Increasing the Reynolds number, Eu increases, and for Re = 1000 (see Fig. 6(c)), this error is always greater than 10−2 for all the simulations performed, while with DVH this error still remains of order 10−4 for the highest resolution. In the SPH method vorticity is not a primary variable and to evaluate this field a further integral interpolation is needed through the formula:

ρω(xi ) = e3 ·

 j

(uj − ui ) × ∇ Wj (xi ) mj .


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Fig. 6. SPH simulations: error on the velocity field Eu for the Lamb–Oseen problem for three different Reynolds numbers and using three different ratios 2ε/1x = 2, 4, 8.

Fig. 7. SPH simulations: error on the vorticity field Eω for the Lamb–Oseen problem for three different Reynolds numbers and using three different ratios 2ε/1x = 2, 4, 8.

Fig. 8. SPH simulations: vorticity field at time t = ω0 /Re for the Lamb–Oseen problem for three different Reynolds numbers.

Therefore, larger errors on ω are expected with respect to that on velocity. This is confirmed by the results depicted in Fig. 7(a), (b) and (c) where Eω is one order of magnitude larger that Eu (see Fig. 6(a), (b) and (c)). In particular for Re = 100 and 1000 the error Eω is even always divergent when using the smallest ratio 2ε/1x. Fig. 8 shows the vorticity field at the end of the simulations for the spatial resolution L/(2ε) = 3 and for the ratio 2ε/1x = 8, three Reynolds numbers have been analysed: 10, 100, 1000. From these plots it is quite visible that increasing Re spurious numerical vorticity starts to form inducing errors in the numerical solution. The development of this ‘‘numerical turbulence’’ and its connection to extra-dissipation mechanism has been already commented in [24] in the context of viscous gravity waves. In left plot of Fig. 9(a) the time history of the SPH kinetic energy for Re = 1000 is compared with the analytical and DVH ones. At the end of the simulation the SPH model has dissipated about 10% more energy than that expected while the DVH method provides a plot that is practically superimposed to the analytical one.


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 9. Left: Time histories of the kinetic energy Ek for the Lamb–Oseen problem at Re = 1000, the SPH and DVH results are compared with the analytic solution. Right: Time histories of the SPH total circulation Γ for the Lamb–Oseen problem for three different Reynolds numbers.

Fig. 10. Lamb–Oseen problem at Re = 1000, effect of the parameter β and p∞ on the SPH solution. (a): SPH Vorticity field at time t = ω0 /Re using β = 0 and p∞ = ρ0 ω02 L2 log(2)/4. (b): Time histories of the kinetic energy Ek . (c): Time histories of the SPH total circulation Γ .

Right plot of Fig. 9(b) displays the time history of the total circulation error. For Re = 1000 the error on the conservation of the total circulation, Γ , starts to become non-negligible going beyond 0.5% (regarding the conservation of circulation in SPH method see also [34]). 4.1.3. Numerical results for SPH using a constant background pressure Conversely to the previous section, here the Lamb–Oseen problem is solved with the SPH method without using the radial force field, i.e. assuming β = 0. However, to avoid the tensile instability a constant background pressure p∞ = ρ0 ω02 L2 log(2)/4 is considered to guarantee positive values of the pressure field (43). In this section only the case Re = 1000 is considered being the case where the SPH numerical solutions, presented in the previous section, show the highest error levels. Fig. 10(a) shows the vorticity field at the end of the simulation at t ω0 /Re = 1. Using the parameters [β, p∞ ] = [0, ρ0 ω02 L2 log(2)/4], the vorticity field is much less noisy with respect to the case [β, pR ] = [0.1ω0 , 0] (see right plot of Fig. 8). This affects also the energy decay and the circulation conservation which are now much closer to the analytical solution as depicted in plots 10(b) and (c). Plots 11(a), (b) depict the relative errors Eu , Eω for the velocity and vorticity fields. For the two highest values of the ratio (2ε)/1x, 4 and 8, these errors decrease almost of one order of magnitude with respect to those shown in Section 4.1.2.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Fig. 11. Lamb–Oseen problem at Re = 1000. SPH Error on the velocity field Eu (a) and on vorticity field Eω (b) using three different ratios 2ε/1x = 2, 4, 8.

Fig. 12. Pressure profile for the Lamb–Oseen problem at time t ω0 /Re = 0.5 for Re = 1000. Top: SPH and DVH results are compared with the analytic solution. Bottom: Comparison of the SPH results using the ratios (2ε)/1x equal to 4 and 8. SPH results are obtained using β = 0 and p∞ = ρ0 ω02 L2 log(2)/4 (see Section 4.1).

The improvement of the SPH results is due to the lower intensity of the pressure field with respect to that obtained using the radial force field with β = 0.1ω0 (see Fig. 2). The over-pressure induced by the force field enforces the particles to reorganize their spatial positions, therefore, high frequency modes have been excited producing numerical turbulence and extra-dissipation. This mechanism is more deeply discussed in [33]. Finally, top plot of Fig. 12 shows the pressure profile at t ω0 /Re = 0.5 obtained through the analytical formula (43), the DVH and the SPH solvers. For the latter a spatial resolution L/(2ε) = 3 and the ratio 2ε/1x = 8 have been adopted. The SPH profile matches the analytical solution but close to the vortex core a non-physical pressure-drop is detected. When using the DVH solver the pressure profile has been calculated by formula (28). Adopting the ratios L/Rd = 3 and Rd /1x = 4 the DVH pressure profile is pretty well superposed to the analytical solution. On the bottom plot of Fig. 12 also the SPH pressure profile obtained with the ratio 2ε/1x = 4 has been reported. Using such a ratio, numerical high-frequency pressure noise appears in the SPH solution. 4.1.4. Numerical results for SPH considering fluid of different nature Despite of the drawbacks of the SPH method discussed in the previous section, this method has the advantage of a direct evaluation of the fluid deformation due to its pure Lagrangian nature. Fig. 13 shows the fluid deformation at the end of the simulation t ω0 = Re for the three Reynolds numbers considered. In all the three cases the same amount of vorticity is damped, however, from these snapshots, it is quite visible the increasing relevance of the advection transport with the Reynolds number. The direct control of the fluid deformation can be quite useful for many applications. For example, on the two particle sets highlighted in Fig. 13 it is easy to give different fluid properties (few changes are needed) with SPH; while this operation it is generally not straightforward using vortex methods. The results reported in Fig. 14 refer to a case where two different fluids are used. The Reynolds number of the first fluid is Re1 = π ω0 L2 /ν1 = 104 while the second fluid has a Reynolds number which changes in time, starting with the same


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 13. SPH simulations: fluid deformation for the Lamb–Oseen problem for three different Reynolds numbers. Particles coloured in blue were initially positioned on x < 0 semi-plane while red particle belonged to x > 0 at t = 0. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Fig. 14. SPH simulations: mixing of two fluids with different viscosity initialized with the Lamb–Oseen velocity–pressure fields. Left: Fluid Deformations: Fluid-1 particles are coloured in blue and are initially positioned on x < 0 semi-plane. Red particles belonged to x > 0 at t = 0 and are relative to the Fluid-2 which viscosity changes in time (see Eq. (47)). Right: Vorticity field at time t ω0 /Re2 = 20. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

Re2 (t = 0) = Re1 and decreasing linearly in time by the factor 1/1000 at a fixed time t0 , Re2 (t0 ) = π ω0 L2 /ν2 = 10:

Re(t ) =

 t  Re1 + (Re2 − Re1 )

t < t0

Re   2 t0 = 2 Re2 /ω0 ;

t ≥ t0 Re1 = 104 , Re2 = 10.



Because of this change in the viscosity of the second fluid the vorticity field becomes more complicated than that presented for the Lamb–Oseen problem. Indeed, vorticity is generated by the interaction between the fluids and also negative values develop (see right plot of Fig. 14). This is just an example to show the flexibility of SPH with respect to a Vortex Method. Even if the DVH allows to get more accurate results with lower CPU costs, it is important to underline the advantages of the SPH method like in this simple example just discussed, for which the DVH scheme requires complex changes in the model while for the SPH the modifications are trivial (for SPH simulations of multi-fluids dynamics see also [36]). 4.2. Merger of a pair of co-rotating vortices After the simple dynamics of an isolated vorticity patch discussed in the previous section, the more complex behaviour of a merger of a pair of co-rotating vorticity patches is considered here. At large Reynolds numbers, the dynamics of the merging process is not significantly influenced by the vorticity diffusion, the phenomenon is almost driven by the advection term while the diffusion process has a secondary role and could be relevant only for very low Reynolds numbers. In the first snapshot of Fig. 15 the initial condition of the problem is shown. The vorticity is equal to ω0 inside two circular regions of radius R and zero elsewhere. The total circulation of each patch at t = 0 is equal to Γ0 = π R2 ω0 , the Reynolds number for this problem is defined as Re = Γ0 /ν and it is set equal to 18 850, in order to get a low-viscosity evolution (for more details see [16]). The results presented in Fig. 15 show the vorticity evolution obtained by the DVH model using a spatial discretization R/Rd = 50 (Rd /1x = 4, ξ = 10−5 ) which corresponds to a number of vortex particles equal to about 251.000. At the end of the simulation, because of the diffusion process, the number of the vortex particles is about 2 millions. The two initial vortical regions are highly stretched while rotating around the origin of the axis. At time t ω0 = 30 the two

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Fig. 15. Merging of a pair of co-rotating vorticity patches: vorticity field evolution using DVH solver.

patches are almost merged while two thin vortical structures are shed. To capture the latter the spatial resolution needs to be properly selected. For such complex evolution an analytical solution is not possible to be obtained, however, in order to assess the validity of the DVH model and to monitor the accuracy of the solutions, some first integrals of motion involving global quantities can be computed:

• • • •

total circulation: Γ = Ω ω dx,  second vorticity moment: I = Ω r 2 ω dx,  enstrophy: S = 1/2 Ω ω2 dx,  excess energy E = 1/2 Ω ω ψ dx,

being ψ the stream function. For the present DVH model the conservation of the total circulation Γ is in practice always preserved (relative errors less than 10−6 since we use a non null Γcutoff , see Section 2.3). For the problem studied here (2D unbounded without free-stream velocity), the following relations between Γ , I, S and E hold (see e.g. [21]): dI dt


4 Γ0 Re

Γ (t ),

dE dt


2 Γ0 Re

S (t ).


Therefore, the accuracy of the numerical models (DVH and SPH) can be inferred by computing the relative errors:

    I(t ) − I − 4Γ0  t Γ (τ )dτ   0   0 Re   E = max    I t ∈(0,tend ]  I0    t 2Γ0     E = max  E (t ) − E0 + Re 0 S (τ )dτ .     E  t ∈(0,tend ] E0


Left plot of Fig. 16 shows the convergence trend of the errors EI and EE obtained by the DVH model. An almost second order convergence is achieved for both the errors as a function of the spatial discretization R/Rd . On the right plot of the same figure the SPH results for EE are shown as a function of the spatial discretization R/(2ε). Since the two ratios Rd /1x and 2ε/1x are both equal to 4, the two different resolutions, R/Rd for the DVH model, and R/(2ε) for the SPH, can be considered equivalent (same R/1x). As expected the SPH error level is higher than the DVH one, however the second order convergence is achieved also by this solver. Regarding EI , the SPH presents a larger error without any reduction for increasing resolution R/(2ε), as reported in the left plot of Fig. 17. This is mainly due to the development of spurious vorticity in the whole domain (see right plot of Fig. 17). Even if this vorticity amount is very limited the x2 factor inside the definition of I amplifies its effect. The development of this numerical vorticity by the SPH method is widely documented in the literature (see e.g. [24,37]) and can be controlled by increasing the ratio ε/1x.


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 16. Merger of a pair of co-rotating vorticity patches. Left: DVH errors EI and EE (see Eq. (49)) as a function of the spatial discretization R/Rd . Right: errors EE comparison between DVH and SPH results. The SPH spatial resolution is indicated by the ratio R/(2ε).

Fig. 17. Merger of a pair of co-rotating vorticity patches. Left: Errors EI comparison between DVH and SPH results varying the spatial resolutions R/Rd for the DVH and R/(2ε) for the SPH. Right: Vorticity field at t ω0 = 50 using the SPH solver.

We underline that, in order to get the velocity–pressure fields at t = 0 for the SPH, those fields have been evaluated using the DVH solver through formulas (11) and (28), once the initial condition on the vorticity field was prescribed. Similarly to the case discussed in Section 4.1.3 the pressure at infinity, p∞ , has been chosen in order to have the whole pressure field positive in the initial condition. A circular domain of radius 10R has been discretized for the SPH method using the same technique discussed in Section 4.1 for handling the outer boundary. The highest adopted spatial resolution was R/(2ε) = 25 corresponding to a number of particles of about 3 millions. When using the DVH model, simulations up to a spatial resolution R/(2ε) = 50 were possible, thanks to the discretization of the rotational portion of the fluid domain only, allowing also to simulate a last case with Reynolds number equal to 100 000 (see last plot in Fig. 18). Finally, Fig. 18 shows the vorticity field at time t ω0 = 36 obtained by the DVH method, using four different Reynolds number starting from Re = 188.5 up to Re = 100,000. Fig. 19 shows the vorticity field at time t ω0 = 100 for the two highest Reynolds numbers Re = 18 850 and Re = 100,000. At these lowest viscosity levels, vortex filaments develop and intersect each other during the time evolution producing a complex vorticity pattern. It is just an example to show the capability of the DVH in simulating problems for a wide range of the viscosity scales. For the most viscous case the vorticity has been diffused over a big portion of the computational domain already at t ω0 = 36 while for the highest Reynolds number the enstrophy is still close to 90% of its initial value. The two-branch filaments originated during the merger process become thinner and thinner as the Reynolds increases and the vorticity maintains values very close to the initial one inside these regions which become more persistent during the time evolution. 5. Computational resources All the simulations presented in this work have been performed on Desktop PC with 6 cores Xeon(R) CPU X5650, 2.3 GHz, using a simple OpenMP parallel programming for both SPH and DVH solvers. The computational cost of the implemented DVH method is equal to about 80 µs for vortex particle and for time iteration on a single core. For the SPH the cost is 20 µs for fluid particle and for time iteration, again on a single core. The max number of particles used for both methods is of order 106 and the memory allocated for both DVH and SPH did not exceed 1 GB.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


Fig. 18. Merger of a pair of co-rotating vorticity patches. Effect of the Reynolds number on the vorticity field at t ω0 = 36 using the DVH solver.

Fig. 19. Merger of a pair of co-rotating vorticity patches. Effect of the Reynolds number on the vorticity field at t ω0 = 100 using the DVH solver.

Table 1 Numerical parameters for DVH and SPH solvers to simulate the merger of a pair co-rotating vortices at Re = 18850.

R/1x Rd /1x, 2ε/1x N particles N iterations tend ω0 1t ω0 CPU time on six cores



100 4 500,000 3086 50 16.2 10−3 6h

100 4 3,145,000 13 333 50 3.75 10−3 40 h

Considering the merger of a pair of co-rotating vortices, Table 1 reports the main parameters of the simulations performed with both SPH and DVH. The lower CPU costs of the DVH solver are due to the smaller number of particles needed (only the rotational part of the fluid domain needs to be discretized) as well as to the minor number of time iterations (larger time steps 1t).


E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503

Fig. 20. Distribution of point around two objects obtained with the packing algorithm described in [26].

6. Extension of the DVH for solving viscous flow around a body In this section a brief note on how the DVH model can be extended for simulating vortex shedding around bodies of arbitrary geometries is given. A key point of the DVH scheme is that the viscous diffusion process is used to regularize the vortex particle distribution using an ad hoc set of points. In this work we use a simple Cartesian lattice but in general other distribution can be used. For example when considering the flow around a body with complex geometry, an adaptive stencil can be used to regularize the geometric distribution of points around the body. In particular we used the so-called ‘‘glass-configurations’’ which can be obtained through the packing algorithm described in [26]. This algorithm can be easily embedded in the DVH scheme since it is based on a simple particle dynamic system. An example of ‘‘glass-configurations’’ of points around objects is given in Fig. 20. The main characteristic of the ‘‘glass-configurations’’ is that the inter-distance between the points is close to the mean global value. Therefore, it is possible to associate almost a same elementary volume to all the points even if the tessellation of the space is not perfectly regular. This geometric property is used in the diffusion process to derive the vorticity field from the circulation of the vortex particle system. When a solid body is present during the diffusion process, the boundary condition on the vorticity ∂ω/∂ n = 0 needs to be enforced. This can be handled, if the considered vortex is close to the body (i.e. the distance of the vortex from the body is less than Rd ), by using a straight wall approximation together with a distribution of image vortices inside the body. If the body has geometrical singularities an adaptation of this method to account for the presence of edges must be implemented. The velocity of the ith particle can be approximated through: u(xi , t ) = u∞ + ∇ Φ (xi , t ) +

N 

Γj (t )Kε (xi , xj ),


j =1

being u∞ the free stream velocity and ∇ Φ the potential flow velocity which allows for the enforcement of the no-flow condition on the solid surfaces. Concerning the no-slip boundary condition the method proposed in Chorin [14] has been followed. The fluid adherence to the solid surfaces is enforced by generating at each time step a new set of vortices on the body contour whose strength is set to nullify the tangential velocity component on the body surface. 7. Conclusions The numerical simulation of two dimensional free vorticity dynamics is studied using two different particle methods: the Diffused Vortex Hydrodynamics (DVH) and the Smoothed Particle Hydrodynamics (SPH) method. Both the methods present advantages and drawbacks and this work is mainly dedicated to discuss such aspects. The DVH model is directly linked to the vorticity formulation of the Navier–Stokes equations, the pressure field is no more an unknown of the problem and the fluid domain is discretized by vortex particles which are present only in the rotational region. Conversely, the Smoothed Particle Hydrodynamics method solves the problem for the velocity–pressure unknowns and this implies the discretization of the whole fluid domain. However, thanks to its Lagrangian nature, SPH directly controls the fluid domain deformations. This feature allows for the possibility to treat free-surface and multi-phase flows (not discussed in this work) in a natural way with respect to other numerical methods. The DVH presented in this work is characterized by the use of a set of regular points to perform the vorticity diffusion process. This vortex redistribution process avoids excessive clustering or rarefaction of the vortex particles ensuring robustness, convergence and high accuracy of the method. This is a novelty with respect to many vortex methods presented in the literature. The last section was dedicated to the validation of the proposed numerical scheme and to the heuristic check of their convergence. For this reason a simple Lamb–Oseen vortex is firstly used to compare the numerical results against an analytical solution. Then the more complex problem of the merger of two co-rotating vortices at high Reynolds number is investigated.

E. Rossi et al. / Computers and Mathematics with Applications 69 (2015) 1484–1503


In all the proposed tests the present vortex method shows a convergence always close to second order. Regarding the SPH method the problems studied here are quite demanding since the second order operators are based on a low order integral interpolation formula which is quite sensitive to particle disorder. This leads to the generation of spurious vorticity within the fluid domain (see [24]) which implies higher levels of error with respect to DVH. For SPH an order of convergence between one and two is achieved. The SPH drawbacks are compensated by the direct evaluation of the fluid deformations due to its pure Lagrangian nature, which permits the application of the method to a larger set of fluid dynamic problems in a more easy way. Acknowledgements The research leading to these results has received funding by the Flagship Project RITMARE – The Italian Research for the Sea – coordinated by the Italian National Research Council and funded by the Italian Ministry of Education, University and Research within the National Research Program 2014–2015. The authors thank Benjamin Bouscasse from CNR-INSEAN for his suggestions and contributions on the research presented in this paper. References [1] [2] [3] [4] [5] [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37]

A. Colagrossi, G. Graziani, M. Pulvirenti, Particles for fluids: SPH versus vortex methods, Math. Mech. Complex Syst. 2 (1) (2014) 45–70. P. Koumoutsakos, Multiscale flow simulations using particles, Annu. Rev. Fluid Mech. 37 (2005) 457–487. J.D. Eldredge, T. Colonius, A. Leonard, A vortex particle method for two-dimensional compressible flow, J. Comput. Phys. 179 (2) (2002) 371–399. G.H. Cottet, P.D. Koumoutsakos, Vortex Method: Theory and Practice, Cambridge University Press, 2000. L.A. Barba, A. Leonard, C.B. Allen, Advances in viscous vortex methods-meshless spatial adaption based on radial basis function interpolation, Internat. J. Numer. Methods Fluids 47 (5) (2005) 387–421. A.J. Chorin, Numerical study of slightly viscous flow, J. Fluid Mech. 57 (1973) 785–796. P. Ramachandran, M. Ramakrishna, S.C. Rajan, Efficient random walks in the presence of complex two-dimensional geometries, Comput. Math. Appl. 53 (2007) 329–344. J.D. Eldredge, Numerical simulation of the fluid dynamics of 2D rigid body motion with the vortex particle method, J. Comput. Phys. 221 (2007) 626–648. M.G. Benson, P.G. Bellamy-Knights, J.H. Gerrard, I. Gladwell, A viscous splitting algorithm applied to low Reynolds number flows round a circular cylinder, J. Fluids Struct. 3 (1989) 439–479. L. Lucy, A numerical approach to testing the fission hypothesis, Astron. J. 82 (1977) 1013. J.J. Monaghan, Smoothed particle hydrodynamics, Annu. Rev. Astron. Astrophys. 30 (1992) 543. J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (2) (1994) 399–406. A. Colagrossi, M. Antuono, D. Le Touzé, Theoretical considerations on the free-surface role in the smoothed particle hydrodynamics model, Phys. Rev. E 79 (5) (2009) 056701. 1–13. A.J. Chorin, Vortex sheet approximation of boundary layers, J. Comput. Phys. 27 (1978) 428–442. J.D. Eldredge, A. Leonard, T. Colonius, A general deterministic treatment of derivatives in particle methods, J. Comput. Phys. 180 (2) (2002) 686–709. G. Graziani, M. Ranucci, R. Piva, From a boundary integral formulation to a vortex method for viscous flows, Comput. Mech. 15 (1995) 301–314. L. Van Dommelen, E.A. Rundesteiner, Fast, adaptive summation of point forces in the two dimensional Poisson equation, J. Comput. Phys. 83 (1989) 126–147. R.H. Kraichnan, D. Montgomery, Two-dimensional turbulence, Rep. Progr. Phys. 43 (1980) 547–619. H.D. Thompson, B.W. Webb, J.D. Hoffman, The cell Reynolds number myth, Internat. J. Numer. Methods Fluids 5 (1985) 305–310. D.J. Price, Resolving high Reynolds numbers in smoothed particle hydrodynamics simulations of subsonic turbulence, Mon. Not. R. Astron. Soc. 420 (2012) L33–L37. G. Riccardi, D. Durante, Elementi di Fluidodinamica, Springer-Verlag, 2006. D. Molteni, A. Colagrossi, G. Colicchio, On the use of an alternative water state equation in SPH, in: Proc. 2nd International SPHERIC Workshop, Madrid, Spain, 23–25 May 2007. M. Antuono, A. Colagrossi, S. Marrone, D. Molteni, Free-surface flows solved by means of SPH schemes with numerical diffusive terms, Comput. Phys. Comm. 181 (3) (2010) 532–549. A. Colagrossi, A. Souto-Iglesias, M. Antuono, S. Marrone, Smoothed-particle-hydrodynamics modeling of dissipation mechanisms in gravity waves, Phys. Rev. E 87 (2013) 023302. J.J. Monaghan, R.A. Gingold, Shock simulation by the particle method SPH, J. Comput. Phys. 52 (1983) 374–389. A. Colagrossi, B. Bouscasse, M. Antuono, S. Marrone, Particle packing algorithm for SPH schemes, Comput. Phys. Comm. 183 (8) (2012) 1641–1653. R. Di Lisio, M. Grenier, M. Pulvirenti, The convergence of the SPH method, Comput. Math. Appl. 35 (1–2) (1998) 95–102. A. Colagrossi, M. Antuono, A. Souto-Iglesias, D. Le Touzé, Theoretical analysis and numerical verification of the consistency of viscous smoothedparticle-hydrodynamics formulations in simulating free-surface flows, Phys. Rev. E 84 (2011) 026705. F. Macià, J.M. Sánchez, A. Souto-Iglesias, L.M. González, WCSPH viscosity diffusion processes in vortex flows, Internat. J. Numer. Methods Fluids 69 (2012) 509–533. N.J. Quinlan, M. Basa, M. Lastiwka, Truncation error in mesh-free particle methods, Int. J. Numer. Methods Eng. 66 (2006) 2064–2085. R. Fatehi, M.T. Manzari, Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives, Comput. Math. Appl. 61 (2) (2011) 482–498. J.W. Swegle, D.L. Hicks, S.W. Attaway, Smoothed particle hydrodynamics stability analysis, J. Comput. Phys. 116 (1995) 123–134. D. Le Touzé, A. Colagrossi, G. Colicchio, M. Greco, A critical investigation of smoothed particle hydrodynamics applied to problems with free-surfaces, Internat. J. Numer. Methods Fluids 73 (7) (2013) 660–691. M. Antuono, A. Colagrossi, D. Le Touzé, J.J. Monaghan, Conservation of circulation in SPH for 2D free-surface flows, Internat. J. Numer. Methods Fluids 72 (2013) 583–606. S. Mas-Gallic, P. Raviart, A particle method for first-order symmetric systems, Numer. Math. 51 (1987) 323–352. N. Tofighi, M. Yildiz, Numerical simulation of single droplet dynamics in three-phase flows using ISPH, Comput. Math. Appl. 66 (4) (2013) 525–536. M. Ellero, P. Español, N.A. Adams, Implicit atomistic viscosities in smoothed particle hydrodynamics, Phys. Rev. E 82 (4) (2010) 046702-6.