A Second Order Semi-Lagrangian Particle Finite Element Method for Fluid Flows

<!-- metadata commented in wiki content


==A Second Order Semi-Lagrangian Particle Finite Element Method for Fluid Flows==

Jonathan Colom-Cobb<sup>(a,*)</sup>, Julio Garcia-Espinosa<sup>(a,b)</sup> , Borja Servan-Camas<sup>(a)</sup> Nadukandi, P.<sup> (a,c)</sup>;

:(1) Centre Internacional de Mètodes Numèrics a l'Enginyeria (CIMNE)

Edifici C1, Campus Norte, UPC, Gran Capitán s/n, 08034 Barcelona, Spain

:(2) Universidad Politécnica de Cataluña (UPC)

C. Gran Capitan s/n, Campus Nord, 08034 Barcelona, Spain.

:(3) Repsol Technology Lab

Ctra. Extremadura km 18, 28935 Móstoles, Madrid.

(*) Corresponding author.
-->

===Abstract===

In this paper, a second order SL-PFEM scheme for solving the incompressible Navier-Stokes equations is presented. This scheme is based on the second order velocity Verlet algorithm, which uses an explicit integration for the particle’s trajectory and an implicit integration for the velocity. The algorithm is completed with a predictor-multicorrector scheme for the integration of the velocity correction using the Finite Element Method. A second order projector based on least squares is used to transfer the intrinsic variables information from the particles onto the background mesh, while a second order interpolation scheme is used to transfer the accelerations from the mesh to the particles. Convergence analyses are carried out to assess the second order convergence.

<span id='_GoBack'></span>'''Keywords: '''Semi-Lagrangian, Particles, PFEM, FEM, Second-order

===1. Introduction===

Solving the convective terms in transport equations using Lagrangian particles avoids the instabilities introduced by the standard approaches in the Eulerian framework. Instead, it leads to the problem of solving the particles trajectories, which offers minimum numerical erosion. However, Lagrangian and semi-Lagrangian methods are usually based on first order in time integration schemes, which limits its rate of convergence. As an exception, a recent work [<span id='cite-_Ref535245996'></span>[[#_Ref535245996|4]]] within the semi-Lagrangian Particle Finite Element Method (SL-PFEM) framework, has proposed a second order scheme. It is based on Strang´s symmetric operator splitting, a third order projector to transfer data from the particle to the mesh, and on estimating the particle´s trajectories using a linear approximation to the instantaneous streamline determined by the velocity field.

<span id='_Ref532985341'></span><span id='_Ref535245996'></span><span id='_Ref532897341'></span><span id='_Ref536171694'></span><span id='_Ref536171696'></span><span id='_Ref535310419'></span>Until [<span id='cite-_Ref535245996'></span>[[#_Ref535245996|4]]], the SL-PFEM framework has been developed based on explicit integrators that compute streamlines to approximate the particle´s trajectories. This approach has been successfully applied to many different problems [<span id='cite-1'></span>[[#1|1]], <span id='cite-2'></span>[[#2|2]], <span id='cite-3'></span>[[#3|3]], <span id='cite-4'></span>[[#4|4]]]. In fact, the more convection dominates the flow, the better the streamlines approximate the pathlines and the larger the time step that can be used. However, in many problems, convection does not dominate and streamlines computed explicitly at the beginning of the time step are not a good approximation of pathlines. In these cases, small time steps are required to obtain accurate results.  An example of this was presented in [<span id='cite-5'></span>[[#5|5]]], where a linear wave problem was used as an example. <span id='cite-_Ref532897315'></span>[[#_Ref532897315|Figure 1]] has been extracted from [<span id='cite-_Ref532897341'></span>[[#_Ref532897341|5]]], and it shows how pathlines can be approximated by explicit streamlines and by an explicit acceleration method. It can be observed that the explicitly streamline is not such a good approximation and therefore it requires to use very small time steps. On the other hand, including information of the acceleration, even in an explicit way, can considerably improve the approximation. This is the motivation of this work, where an algorithm based on the second order velocity Verlet integrator [<span id='cite-6'></span>[[#6|6]], <span id='cite-7'></span>[[#7|7]], <span id='cite-8'></span>[[#8|8]]] is proposed, since it has some desirable properties such as being second order accurate over time, time reversible, and symplectic when applied to Hamiltonian systems [<span id='cite-_Ref535310419'></span>[[#_Ref535310419|8]]].

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image1.png|600px]] </div>

<div id="_Ref532897315" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 1: Approximation of pathlines of fluid particles under an acceleration field induced by a linear wave. Top; using explicit streamlines. Down: computing trajectories with explicit acceleration. Figure extracted from [</span><span id='cite-_Ref532897341'></span>[[#_Ref532897341|<span style="text-align: center; font-size: 75%;">5</span>]]<span style="text-align: center; font-size: 75%;">]</span></div>

This paper is organized as follows: first, the problem of particles moving with a velocity field is introduced. Second, the velocity Verlet scheme is presented to integrate the particle´s equations of motion. Third, the incompressible Navier-Stokes equations are introduced in a semi-Lagrangian framework, where the divergence free condition is enforced using a predictor-multicorrector scheme inspired in the fractional step method. Fourth, a global least-square projector is analysed. And fifth, a number of convergence verification, validation and demonstration cases are carried out.

===2. Particles moving with a velocity field===

Along the paper, Eulerian variables (mesh values) are written with lower case letters, particle´s variables are written with upper case letters and vectors are written with bold letters. Let <math display="inline">\mathit{\boldsymbol{u}}(\mathit{\boldsymbol{x}},t)</math> be a velocity field, let <math display="inline">\left\{ \lambda \right\}</math> ''' '''be a set of particles each of them identified with a label <math display="inline">\lambda</math> . The position occupied by each particle at time <math display="inline">t</math> is denoted by <math display="inline">{\mathit{\boldsymbol{X}}}_{\lambda }\mathit{\boldsymbol{(}}t\mathit{\boldsymbol{)}}</math>, and the velocity of the particle by <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }\mathit{\boldsymbol{(}}t\mathit{\boldsymbol{)}}</math>. If the particle velocity is given by the velocity field <math display="inline">\mathit{\boldsymbol{u}}</math>, then''' ''' <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }\left( t\right) \mathit{\boldsymbol{=\, }}\mathit{\boldsymbol{u}}({\mathit{\boldsymbol{X}}}_{\lambda }\mathit{\boldsymbol{(}}t\mathit{\boldsymbol{)}},t)</math>, and the particles’s trajectories are the pathlines defined by the velocity field <math display="inline">\mathit{\boldsymbol{u}}</math>. The particle’s equations of motion are:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{d{\mathit{\boldsymbol{X}}}_{\lambda }\left( t\right) }{dt}={\mathit{\boldsymbol{U}}}_{\lambda }\left( t\right) =</math><math>\mathit{\boldsymbol{u}}({\mathit{\boldsymbol{X}}}_{\lambda }(t),t)</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref7798987'></span>(1)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{d{\mathit{\boldsymbol{U}}}_{\lambda }\left( t\right) }{dt}={\mathit{\boldsymbol{A}}}_{\lambda }(t)</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref7798989'></span>(2)
|}


where the particle’s acceleration <math display="inline">{\mathit{\boldsymbol{A}}}_{\lambda }\left( t\right)</math>  is the rate of change of <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }\mathit{\boldsymbol{(}}t\mathit{\boldsymbol{)}}</math>. An acceleration field <math display="inline">\mathit{\boldsymbol{a}}(\mathit{\boldsymbol{x}},t)</math> can be derived from the velocity field as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\mathit{\boldsymbol{a}}\left( \mathit{\boldsymbol{x}},t\right) ={d}_{t}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{x}},t\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(3)
|}


where <math display="inline">{d}_{t}</math> is the total derivative

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{d}_{t}\mathit{\boldsymbol{u}}={\partial }_{t}\mathit{\boldsymbol{u}}+</math><math>\mathit{\boldsymbol{u}}\nabla \mathit{\boldsymbol{u}}</math>'''.'''
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527448368'></span>(4)
|}


Given an acceleration field, <math display="inline">\mathit{\boldsymbol{a}}</math>,''' '''derived from a velocity field <math display="inline">\mathit{\boldsymbol{u}}</math> as in Eq. <span id='cite-_Ref527448368'></span>[[#_Ref527448368|(4)]], the particles’ velocities and positions can be obtained integrating the equations of motion with initial conditions:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }\left( 0\right) =\mathit{\boldsymbol{u}}({\mathit{\boldsymbol{X}}}_{\lambda }(0),0)</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(5)
|}


===3. Numerical integration of particles’ equations of motion===

Let <math display="inline">\mathit{\boldsymbol{a}}\left( \mathit{\boldsymbol{x}},t\right) =</math><math>{d}_{t}\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{x,}}t\right)</math>  be an acceleration field derived from an unknown velocity field <math display="inline">\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{x,}}t\right)</math> . Given a set of particles <math display="inline">\left\{ \lambda \right\}</math>  with initial position and velocity <math display="inline">{\mathit{\boldsymbol{X}}}_{\lambda }(0)</math> and <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }\left( 0\right) =</math><math>\mathit{\boldsymbol{u}}({\mathit{\boldsymbol{X}}}_{\lambda }(0),0)</math>. Then the position and velocity at time <math display="inline">t</math> is obtained integrating the equations of motion:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{X}}}_{\lambda }\left( {t}^{n+1}\, \right) ={\mathit{\boldsymbol{X}}}_{\lambda }\left( {t}^{n}\right) +</math><math>\int_{{t}^{n}}^{{t}^{n+1}}{\mathit{\boldsymbol{U}}}_{\lambda }\left( t\right) dt</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(6)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n+1}\, \right) ={\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n}\right) +</math><math>\int_{{t}^{n}}^{{t}^{n+1}}{\mathit{\boldsymbol{A}}}_{\lambda }\left( t\right) dt</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(7)
|}


In molecular  dynamics, a well-known numerical integration scheme for the equations of motion is the velocity Verlet integrator [<span id='cite-_Ref536171694'></span>[[#_Ref536171694|6]],<span id='cite-_Ref536171696'></span>[[#_Ref536171696|7]]], which reads:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{X}}}_{\lambda }\left( {t}^{n+1}\, \right) ={\mathit{\boldsymbol{X}}}_{\lambda }\left( {t}^{n}\right) +</math><math>\Delta t{\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n}\right) +\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{A}}}_{\lambda }\left( {t}^{n}\right) +</math><math>O(\Delta {t}^{3})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527385792'></span>(8)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n+1}\right) ={\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n}\right) +</math><math>\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{A}}}_{\lambda }\left( {t}^{n}\right) +\right. </math><math>\left. {\mathit{\boldsymbol{A}}}_{\lambda }\left( {t}^{n+1}\right) \right) +O(\Delta {t}^{3})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527385793'></span>(9)
|}


where eq. <span id='cite-_Ref527385792'></span>[[#_Ref527385792|(8)]] is explicit in time and provides the particle position. We can rewrite  eqs. <span id='cite-_Ref527385792'></span>[[#_Ref527385792|(8)]] and <span id='cite-_Ref527385793'></span>[[#_Ref527385793|(9)]] as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{X}}}_{\lambda }^{n}+</math><math>\Delta t{\mathit{\boldsymbol{U}}}_{\lambda }^{n}+\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(10)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{U}}}_{\lambda }^{n}+</math><math>\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{a}}}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})+\right. </math><math>\left. {\mathit{\boldsymbol{a}}}^{n+1}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1})\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(11)
|}


where <math display="inline">{\mathit{\boldsymbol{X}}}_{\lambda }^{n}={\mathit{\boldsymbol{X}}}_{\lambda }\left( {t}^{n}\right)</math> , <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }^{n}=</math><math>{\mathit{\boldsymbol{U}}}_{\lambda }\left( {t}^{n}\right)</math> , <math display="inline">{\mathit{\boldsymbol{a}}}^{n}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n}\right) =</math><math>\mathit{\boldsymbol{a}}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right)</math> .

===4. A Semi-Lagrangian approach for solving the incompressible Navier-Stokes equations===

====1.1 Lagrangian formulation of the incompressible Navier-Stokes equations====

The incompressible Navier Stokes equations can be written as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{d}_{t}\mathit{\boldsymbol{u}}\boldsymbol{=-}\mathit{\boldsymbol{\nabla }}\left( \frac{P}{\rho }\right) +</math><math>\nu \Delta \mathit{\boldsymbol{u+f}}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527449161'></span>(12)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\mathit{\boldsymbol{\nabla }}\cdot \mathit{\boldsymbol{u}}=0</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref535921834'></span>(13)
|}


Where  <math display="inline">P</math> is the fluid pressure, <math display="inline">\rho</math>  is the fluid density, <math display="inline">\nu</math>  is the kinematic viscosity and <math display="inline">\mathit{\boldsymbol{f}}</math>''' '''are mass forces. Eq. <span id='cite-_Ref527449161'></span>[[#_Ref527449161|(12)]] can be expressed as <math display="inline">{d}_{t}\mathit{\boldsymbol{u=}}\mathit{\boldsymbol{a}}</math>,''' '''where the acceleration field is given by:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\mathit{\boldsymbol{a}}\boldsymbol{=-}\mathit{\boldsymbol{\nabla }}\left( \frac{P}{\rho }\right) +</math><math>\nu \Delta \mathit{\boldsymbol{u+f}}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(14)
|}


Solving the Navier Stokes’ equation in a Lagrangian framework by integrating the fluid particles dynamics is not as simple, since the acceleration field depends on the fluid velocity and the pressure. In order to ensure the divergence free condition of the flow field, the pressure must fulfil an integral equation on the domain. This proves to be very complicated when imposing the continuity condition in a Lagrangian framework.

====1.2 Semi-Lagrangian integration of the incompressible Navier-Stokes equations====

Below, a Semi-Lagrangian (SL) approach is introduced to overcome the difficulties of estimating the pressure gradient and viscosity terms in a Lagrangian framework. To do so the acceleration field will be obtained in an Eulerian framework by projecting the particle´s velocity field onto the Eulerian space.

Let <math display="inline">\left\{ \lambda \right\}</math>  be a set of fluid particles. Then the equation of motions can be integrated using the velocity Verlet integrator:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{X}}}_{\lambda }^{n}+</math><math>\Delta t{\mathit{\boldsymbol{U}}}_{\lambda }^{n}+\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527451108'></span>(15)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}-{\mathit{\boldsymbol{U}}}_{\lambda }^{n}}{\Delta t}=</math><math>\frac{1}{2}\left( {\mathit{\boldsymbol{a}}}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})+\right. </math><math>\left. {\mathit{\boldsymbol{a}}}^{n+1}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1})\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527451071'></span>(16)
|}


As stated above, this scheme is second order accurate over time. Let´s assume that all variables are known at time <math display="inline">{t}^{n}</math>. Then the new position of the fluid particles at time <math display="inline">{t}^{n+1}</math> can be easily calculated using Eq.<span id='cite-_Ref527451108'></span>[[#_Ref527451108|(15)]]. Reordering Eq. <span id='cite-_Ref527451071'></span>[[#_Ref527451071|(16)]]:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }^{n+1,i+1}={\mathit{\boldsymbol{U}}}_{\lambda }^{n+1/2}+</math><math>\frac{\Delta t}{2}{\mathit{\boldsymbol{a}}}^{n+1}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref3545805'></span><span id='_Ref527454686'></span>(17)
|}


where <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }^{n+1/2}={\mathit{\boldsymbol{U}}}_{\lambda }^{n}+</math><math>\frac{\Delta t}{2}{\mathit{\boldsymbol{a}}}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})</math>. In order to solve eq. <span id='cite-_Ref527454686'></span>[[#_Ref527454686|(17)]], the term <math display="inline">{\mathit{\boldsymbol{a}}}^{n+1}</math> will be solved in an Eulerian framework. Then an equivalent Eulerian equation is imposed:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{u}}}^{n+1}\left( \mathit{\boldsymbol{x}}\right) ={\mathit{\boldsymbol{u}}}^{n+1/2}\left( \mathit{\boldsymbol{x}}\right) +</math><math>\frac{\Delta t}{2}{\mathit{\boldsymbol{a}}}^{n+1}(\mathit{\boldsymbol{x}})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527457971'></span>(18)
|}


Where <math display="inline">{\mathit{\boldsymbol{u}}}^{n+1/2}</math>  is obtained via projection:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{u}}}^{n+1/2}\left( \mathit{\boldsymbol{x}}\right) ={P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{n+1/2}\right\} \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(19)
|}


<math display="inline">{P}_{\left\{ \lambda \right\} }^{n+1}</math> is a projector operator that transfer the values of the variables from the set of particle with positions <math display="inline">\left\{ {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}\right\}</math>  onto the mesh.

We now introduce what we call the coherence condition for the projector:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">{\psi }^{n+1}\, (\mathit{\boldsymbol{x}})={P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ \psi \left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}\right) \right\} \right]</math> ,
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref7778240'></span>(20)
|}


where <math display="inline">\psi</math>  is an arbitrary variable defined on the Eulerian framework. This condition means that the projection of the nodal values interpolated at the particle´s return the original nodal values. If this condition is fulfilled, then from Eq. <span id='cite-_Ref527457971'></span>[[#_Ref527457971|(18)]] and <span id='cite-_Ref7778240'></span>[[#_Ref7778240|(20)]]:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{u}}}^{n+1}\left( \mathit{\boldsymbol{x}}\right) ={P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{n+1/2}\right\} \right] +</math><math>\frac{\Delta t}{2}{P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{a}}}^{n+1}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}\right) \right\} \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref7798820'></span>(21)
|}


which is equivalent to say that <math display="inline">\, {\mathit{\boldsymbol{u}}}^{n+1}\left( \mathit{\boldsymbol{x}}\right) =</math><math>{P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}\right\} \right]</math> . There is another outcome from the fulfilment of the coherence condition. <span id='cite-_Ref7777391'></span>[[#_Ref7777391|Figure 2]] shows the conceptual implementation of a generic SL-PFEM. In general two iterative loops are required to avoid any sort of splitting error. However, if the coherence condition is fulfilled, there is no need of iterating at the outer loop. This is because the correction of the particles´ velocities, given by  <math display="inline">\frac{\Delta t}{2}{\mathit{\boldsymbol{a}}}^{n+1}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1})</math> , will be projected back to the mesh as it is, with no further need of correction.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image2.jpeg|600px]] </div>

<div id="_Ref7777391" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 2: Conceptual description of the SL-PFEM</span></div>

Eq. <span id='cite-_Ref527457971'></span>[[#_Ref527457971|(18)]] can be solved iteratively by further splitting into:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\hat{\mathit{\boldsymbol{u}}}}^{n+1,i+1}={\mathit{\boldsymbol{u}}}^{n+1/2}+\frac{\Delta t}{2}\left( -\right. </math><math>\left. \nabla \left( \frac{{P}^{n+1,i}}{\rho }\right) +\nu \Delta {\hat{\mathit{\boldsymbol{u}}}}^{n+1,i+1}+\right. </math><math>\left. {\mathit{\boldsymbol{f}}}^{n+1}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527460900'></span>(22)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{u}}}^{n+1,i+1}-{\hat{\mathit{\boldsymbol{u}}}}^{n+1,i+1}=\frac{\Delta t}{2}\left( -\right. </math><math>\left. \nabla \left( \frac{{P}^{n+1,i+1}}{\rho }\right) +\nabla \left( \frac{{P}^{n+1,i}}{\rho }\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527458693'></span>(23)
|}


Where <math display="inline">i</math> is the iteration index. Now in order to obtain the pressure <math display="inline">{P}^{n+1,i+1}</math>, the continuity condition is imposed. Introducing the continuity equation ( <math display="inline">\nabla \cdot {\mathit{\boldsymbol{u}}}^{n+1,i+1}=</math><math>0</math>) into Eq.<span id='cite-_Ref527458693'></span>[[#_Ref527458693|(23)]] and reordering

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\Delta \left( \frac{{P}^{n+1,i+1}}{\rho }\right) =2\frac{\nabla \cdot {\hat{\mathit{\boldsymbol{u}}}}^{n+1,i+1}}{\Delta t}+</math><math>\Delta \left( \frac{{P}^{n+1,i}}{\rho }\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref527460073'></span>(24)
|}


Once the iterative process between Eqs. <span id='cite-_Ref527460900'></span>[[#_Ref527460900|(22)]]-<span id='cite-_Ref527458693'></span>[[#_Ref527458693|(23)]] has converged, from eq. <span id='cite-_Ref527457971'></span>[[#_Ref527457971|(18)]] we obtain:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{a}}}^{n+1}\left( \mathit{\boldsymbol{x}}\right) =2\frac{{\mathit{\boldsymbol{u}}}^{n+1}\left( \mathit{\boldsymbol{x}}\right) -{\mathit{\boldsymbol{u}}}^{n+1/2}\left( \mathit{\boldsymbol{x}}\right) }{\Delta t}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(25)
|}


====1.3 The SL-PFEM method for solving the incompressible Navier-Stokes equations====

In the previous section, a projector operator was introduced to project the particle´s variables onto the Eulerian space. Then, a strategy to solve the acceleration due to the pressure gradient and viscous terms can be formulated in an Eulerian framework.

The finite element method (FEM) is introduced to solve Eq. <span id='cite-_Ref527460900'></span>[[#_Ref527460900|(22)]] and <span id='cite-_Ref527460073'></span>[[#_Ref527460073|(24)]]. First, we need to discretize it in space using a finite element mesh. Then, any intrinsic dependent variable <math display="inline">\psi</math>  can be expressed as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\psi }_{h}(\mathit{\boldsymbol{x}})=\sum _{b}^{}{N}^{b}\left( \mathit{\boldsymbol{x}}\right) {\psi }_{b}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(26)
|}


where <math display="inline">\lbrace b\rbrace</math>  is the set of mesh nodes, <math display="inline">{N}^{b}\left( \mathit{\boldsymbol{x}}\right)</math>  are the classic Finite Element (FE) linear shape functions, and <math display="inline">{\psi }_{b}</math> are the nodal values. And let <math display="inline">{\boldsymbol{\psi }}_{h}</math> be the vector of nodal values <math display="inline">{\psi }_{b}</math>.

Once the spatial discretization is set with a FE mesh, the projector operator must transfer values from particles onto the background mesh following Eq. <span id='cite-_Ref7798820'></span>[[#_Ref7798820|(21)]]:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\boldsymbol{u}}_{h}^{n+1}{=P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}\right\} \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(27)
|}


Where <math display="inline">{P}_{h,\left\{ \lambda \right\} }</math> must fulfil the condition given in Eq. <span id='cite-_Ref7778240'></span>[[#_Ref7778240|(20)]] leading to:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{a}}}_{h}^{n+1}(\mathit{\boldsymbol{x}})={P}_{h,\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{a}}}_{h}^{n+1}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n}\right) \right\} \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(28)
|}


<span id='_Ref6318173'></span>Using the FE Galerkin method as in [<span id='cite-9'></span>[[#9|9]],<span id='cite-10'></span>[[#10|10]]] to solve the variational formulation of Eqs. <span id='cite-_Ref527460900'></span>[[#_Ref527460900|(22)]] and <span id='cite-_Ref527460073'></span>[[#_Ref527460073|(24)]] we get:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left( \overline{\overline{\mathit{\boldsymbol{M}}}}-\frac{\Delta t\nu }{2}\overline{\overline{\mathit{\boldsymbol{L}}}}\right) {\hat{\boldsymbol{u}}}_{h}^{n+1,i+1}=</math><math>\overline{\overline{\mathit{\boldsymbol{M}}}}{\boldsymbol{u}}_{h}^{n+1/2}-\frac{\Delta t}{2\rho }\overline{\overline{\mathit{\boldsymbol{G}}}}{\boldsymbol{P}}_{h}^{n+1,i}+</math><math>\frac{\Delta t}{2}\overline{\overline{\mathit{\boldsymbol{M}}}}{\boldsymbol{F}}_{h}^{n+1,i}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(29)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\overline{\overline{\mathit{\boldsymbol{L}}}}{\boldsymbol{P}}_{h}^{n+1,i+1}=\frac{2\rho }{\Delta t}\overline{\overline{\mathit{\boldsymbol{D}}}}{\hat{\boldsymbol{u}}}_{h}^{n+1,i+1}+</math><math>\overline{\overline{\mathit{\boldsymbol{D}}}}{\overline{\overline{\mathit{\boldsymbol{M}}}}}^{-1}\overline{\overline{\mathit{\boldsymbol{G}}}}{\boldsymbol{P}}_{h}^{n+1,i}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(30)
|}


Where <math display="inline">\overline{\overline{\mathit{\boldsymbol{M}}}}</math> is the FE mass matrix, <math display="inline">\overline{\overline{\mathit{\boldsymbol{L}}}}</math> is the FE Laplacian matrix, <math display="inline">\overline{\overline{\mathit{\boldsymbol{G}}}}</math> is the FE gradient matrix, and <math display="inline">\overline{\overline{\mathit{\boldsymbol{D}}}}</math> is the FE divergence matrix.  Once the iterative loop over <math display="inline">i</math> converges, then the new acceleration field is obtained from Eq. <span id='cite-_Ref527457971'></span>[[#_Ref527457971|(18)]]. The convergence criteria used for the pressure and velocity are:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\left\| {\boldsymbol{P}}_{h}^{n+1,i+1}-{\boldsymbol{P}}_{h}^{n+1,i}\right\| }{\left\| {P}_{max}^{n+1,i+1}-{P}_{min}^{n+1,i+1}\right\| }<0.001</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(31)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{\left\| {\left| {\mathit{\boldsymbol{U}}}_{h}\right| }^{n+1,i+1}-{\left| {\mathit{\boldsymbol{U}}}_{h}\right| }^{n+1,i}\right\| }{{\left| \mathit{\boldsymbol{U}}\right| }_{max}}<0.001</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(32)
|}


===5. Global least square projector===

Let <math display="inline">\left\{ {\Psi }_{\lambda }\right\}</math>  be a set of values carried by a set of particles and <math display="inline">\left\{ {\psi }_{b}^{\ast }\right\}</math>  be the set of projected values on the mesh nodes. The projector operator used in this work is based on the minimization of the global least square error. The global least square error is given by:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{\psi }=\sum _{\lambda }^{}{\left( {\psi }_{h}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) -{\Psi }_{\lambda }\right) }^{2}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(33)
|}


where <math display="inline">{\psi }_{h}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) =</math><math>\sum _{b}^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{b}^{\ast }</math>

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\epsilon }_{\psi }=\sum _{\lambda }^{}{\left( \sum _{b}^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{b}^{\ast }-{\Psi }_{\lambda }\right) }^{2}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(34)
|}


Minimizing the least square error via <math display="inline">\partial {\epsilon }_{\psi }/\partial {\psi }_{b}^{\ast }=</math><math>0</math>:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\sum _{\lambda }^{}\left( \sum _{c}^{}{{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{c}^{\ast }\right) =</math><math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\Psi }_{\lambda }</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref503780164'></span>(35)
|}


Eq. <span id='cite-_Ref503780164'></span>[[#_Ref503780164|(35)]] represents a linear system of equation with as many unknowns as mesh nodes. This system can be seen as a discrete version of the classic FE projection. Introducing an interpolated field on the particles  <math display="inline">{\Psi }_{\lambda }=</math><math>\sum _{c}^{}{N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{c}</math> :

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math display="inline">\sum _{\lambda }^{}\left( \sum _{c}^{}{{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{c}^{\ast }\right) =</math><math>\sum _{\lambda }^{}\left( \sum _{c}^{}{{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{c}\right)</math> ,
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(36)
|}


leading to <math display="inline">{\psi }_{c}^{\ast }={\psi }_{c}</math> and the fulfilment of the coherence condition given in eq. <span id='cite-_Ref7778240'></span>[[#_Ref7778240|(20)]].

In order to analyse the order of convergence of the projector, let’s assume that the particles values are given by a continuous and smooth function <math display="inline">\psi</math> . Then the equalities <math display="inline">{\Psi }_{\lambda }=</math><math>\psi \left( {\mathit{\boldsymbol{X}}}_{\lambda }\right)</math>  and <math display="inline">\psi \left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) =</math><math>\sum _{c}^{}{N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \psi \left( {\mathit{\boldsymbol{x}}}_{c}\right) +</math><math>O({h}_{e}^{2})</math> hold, where <math display="inline">{h}_{e}</math> is the characteristic length of the element where particle <math display="inline">\lambda</math>  is at. Using Eq. <span id='cite-_Ref503780164'></span>[[#_Ref503780164|(35)]]:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \left( \sum _{c}^{}{N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) {\psi }_{c}^{\ast }\right) =</math><math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \left( \sum _{c}^{}{N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \psi \left( {\mathit{\boldsymbol{x}}}_{c}\right) \right) +</math><math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) O({h}_{e}^{2})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref531339507'></span>(37)
|}


And, the solution to eq. <span id='cite-_Ref531339507'></span>[[#_Ref531339507|(37)]] can be expressed as:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\psi }_{c}^{\ast }=\psi \left( {\mathit{\boldsymbol{x}}}_{c}\right) +\delta {\psi }_{c}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(38)
|}


Where <math display="inline">\delta {\psi }_{c}</math> is the solution of:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \left( \sum _{c}^{}{N}^{c}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) \delta {\psi }_{c}\right) =</math><math>\sum _{\lambda }^{}{N}^{b}\left( {\mathit{\boldsymbol{X}}}_{\lambda }\right) O({h}_{e}^{2})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(39)
|}


Then, the projected value converges with second order accuracy in space:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\psi }_{c}^{\ast }=\psi \left( {\mathit{\boldsymbol{x}}}_{c}\right) +O({h}_{e}^{2})</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(40)
|}


===6. Error accumulation in the SL-PFEM===

SL-PFEM methods are prone to accumulate errors as information is passed from the Lagrangian particles onto the background mesh and vice-versa. Those errors are originated in the projection and interpolation steps. And the accumulation of those errors could lead to several undesired effects such as the degradation of the rate of convergence of the numerical scheme, and the incapability to reach a steady state solution. This section aims at demonstrating that the rate of convergence of the proposed scheme is not degraded due to this phenomenon and that second order convergence can be achieved.

Given a velocity field <math display="inline">\mathit{\boldsymbol{u}}\left( \mathit{\boldsymbol{x}},{t}^{n}\right) =</math><math>{\mathit{\boldsymbol{u}}}^{n}\left( \mathit{\boldsymbol{x}}\right)</math>  and the corresponding acceleration field <math display="inline">\mathit{\boldsymbol{a}}\left( \mathit{\boldsymbol{x}},{t}^{n}\right) =</math><math>{\mathit{\boldsymbol{a}}}^{n}\left( \mathit{\boldsymbol{x}}\right) ={d}_{t}{\mathit{\boldsymbol{u}}}^{n}\left( \mathit{\boldsymbol{x}}\right)</math> , and a set of particles with initial condition <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }^{0}=</math><math>{\mathit{\boldsymbol{u}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math> , the velocity Verlet algorithm estimates the velocity and position at <math display="inline">{t}^{n}=</math><math>n\Delta t</math> with second order convergence:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\hat{\mathit{\boldsymbol{X}}}}_{\lambda }^{n+1}={\hat{\mathit{\boldsymbol{X}}}}_{\lambda }^{n}+</math><math>\Delta t{\hat{\mathit{\boldsymbol{U}}}}_{\lambda }^{n}+\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}^{n}\left( {\hat{\mathit{\boldsymbol{X}}}}_{\lambda }^{n}\right) =</math><math>{\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}+O\left( {t}^{n}\Delta {t}^{2}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref3283373'></span>(41)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\hat{\mathit{\boldsymbol{U}}}}_{\lambda }^{n+1}={\hat{\mathit{\boldsymbol{U}}}}_{\lambda }^{n}+</math><math>\frac{\Delta {t}^{2}}{2}\left( {\mathit{\boldsymbol{a}}}^{n}\left( {\hat{\mathit{\boldsymbol{X}}}}_{\lambda }^{n}\right) +\right. </math><math>\left. {\mathit{\boldsymbol{a}}}^{n\mathit{\boldsymbol{+1}}}\left( {\hat{\mathit{\boldsymbol{X}}}}_{\lambda }^{n+1}\right) \right) =</math><math>{\mathit{\boldsymbol{u}}}^{n+1}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}\right) +</math><math>O\left( {t}^{n}\Delta {t}^{2}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref3283374'></span>(42)
|}


Now, we will evaluate the error accumulation over time, by estimating the error in the Lagrangian integration of a given acceleration field and subsequent mesh projection-interpolation. Let us consider a space discretization using a FE mesh, in which the given acceleration field is approximated by an interpolated field, as follows:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{a}}}_{h}^{n}\left( \mathit{\boldsymbol{x}}\right) =\sum _{b}^{}{N}^{b}\left( \mathit{\boldsymbol{x}}\right) \mathit{\boldsymbol{a}}\mathit{\boldsymbol{(}}{\mathit{\boldsymbol{x}}}_{b}\mathit{\boldsymbol{,}}{t}^{n}\mathit{\boldsymbol{)}}=</math><math>\sum _{b}^{}{N}^{b}\left( \mathit{\boldsymbol{x}}\right) {\mathit{\boldsymbol{a}}}_{b}^{n}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(43)
|}


The velocity at the mesh nodes is initiated from the initial particle velocity <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }^{0}</math>''' '''through the projection:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left\{ {\hat{\mathit{\boldsymbol{u}}}}_{b}^{0}\right\} ={P}_{\left\{ \lambda \right\} }^{0}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{0}\right\} \right] =</math><math>{P}_{\left\{ \lambda \right\} }^{0}\left[ \left\{ {\mathit{\boldsymbol{u}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) \right\} \right] =</math><math>\left\{ {\mathit{\boldsymbol{u}}}_{b}^{0}+{\boldsymbol{\epsilon }}^{P,0}\right\}</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(44)
|}


Where <math display="inline">{\boldsymbol{\epsilon }}^{P,0}</math> is the projection error. The projected velocity is then interpolated back to the particles using the projected values:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\hat{\mathit{\boldsymbol{u}}}}_{}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) =</math><math>\sum _{b}^{}{N}^{b}\left( \mathit{\boldsymbol{x}}\right) {\hat{\mathit{\boldsymbol{u}}}}_{b}^{0}=</math><math>\sum _{b}^{}{N}^{b}\left( \mathit{\boldsymbol{x}}\right) \left( {\mathit{\boldsymbol{u}}}_{b}^{0}+{\boldsymbol{\epsilon }}^{P,0}\right) =</math><math>{\mathit{\boldsymbol{u}}}_{h}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>{\boldsymbol{\epsilon }}_{h}^{P,0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(45)
|}


Then the interpolation error has to be taken into account (both for the velocity and acceleration):

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\mathit{\boldsymbol{u}}}_{h}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) =</math><math>{\mathit{\boldsymbol{u}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>{\boldsymbol{\epsilon }}_{h}^{u,0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math>

<math>{\mathit{\boldsymbol{a}}}_{h}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) =</math><math>{\mathit{\boldsymbol{a}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>{\boldsymbol{\epsilon }}_{h}^{a,0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(46)
|}


Finally, the Verlet integrator is used to obtain the particles’ position and velocities:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>{\mathit{\boldsymbol{X}}}_{\lambda }^{0}+\Delta t{\mathit{\boldsymbol{u}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>\Delta t{\epsilon }^{\mathit{\boldsymbol{u}},0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>\Delta t{\boldsymbol{\epsilon }}_{h}^{P,0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +</math><math>\frac{\Delta {t}^{2}}{2}{\epsilon }^{\mathit{\boldsymbol{a}},0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(47)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{U}}}}_{\lambda }^{1}={\mathit{\boldsymbol{U}}}_{\lambda }^{0}+</math><math>\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{a}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +\right. </math><math>\left. {\mathit{\boldsymbol{a}}}^{1}\left( {\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \right) +</math><math>\frac{\Delta t}{2}\left( {\epsilon }^{\mathit{\boldsymbol{a}},0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +\right. </math><math>\left. {\epsilon }^{\mathit{\boldsymbol{a}},1}\left( {\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(48)
|}


Considering that the projector and the interpolator are second order accurate, and using Eq. <span id='cite-_Ref527385792'></span>[[#_Ref527385792|(8)]]-<span id='cite-_Ref527385793'></span>[[#_Ref527385793|(9)]]:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>\underbrace{{{\mathit{\boldsymbol{X}}}_{\lambda }^{0}\mathit{\boldsymbol{+}}\Delta t{\mathit{\boldsymbol{u}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) +\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) }}_{{\hat{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\mathit{\boldsymbol{=}}{\mathit{\boldsymbol{X}}}_{\lambda }^{1}\mathit{\boldsymbol{+}}O(\Delta {t}^{3})\mathit{\boldsymbol{\, }}}+</math><math>\underbrace{{{\Delta t\epsilon }^{\mathit{\boldsymbol{u}},0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) }}_{O(\Delta t{h}_{e}^{2})}+</math><math>\underbrace{{\Delta t{\boldsymbol{\epsilon }}_{h}^{P,0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) }}_{O(\Delta t{h}_{e}^{2})}+</math><math>\underbrace{{{\frac{\Delta {t}^{2}}{2}\epsilon }^{\mathit{\boldsymbol{a}},0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right) }}_{O(\Delta {t}^{2}{h}_{e}^{2})}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(49)
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>\underbrace{{{\mathit{\boldsymbol{U}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{0}}}\mathit{\boldsymbol{+}}\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{a}}}^{\mathit{\boldsymbol{0}}}\left( {\mathit{\boldsymbol{X}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{0}}}\right) \mathit{\boldsymbol{+}}{\mathit{\boldsymbol{a}}}^{\mathit{\boldsymbol{1}}}\left( {\hat{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \right) }}_{{\hat{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\mathit{\boldsymbol{=}}{\mathit{\boldsymbol{U}}}_{\lambda }^{1}\mathit{\boldsymbol{+}}O(\Delta {t}^{3})}+</math><math>\underbrace{{\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{a}}}^{\mathit{\boldsymbol{1}}}\left( {\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \mathit{\boldsymbol{-}}{\mathit{\boldsymbol{a}}}^{\mathit{\boldsymbol{1}}}\left( {\hat{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \right) }}_{\mathit{\boldsymbol{\sim }}\frac{\mathit{\boldsymbol{\Delta t}}}{\mathit{\boldsymbol{2}}}({h}_{e}^{2}\Delta t)\mathit{\boldsymbol{\nabla }}{\mathit{\boldsymbol{a}}}^{\mathit{\boldsymbol{1}}}}+</math><math>\underbrace{{\frac{\Delta t}{2}\left( {\mathit{\boldsymbol{\epsilon }}}^{\mathit{\boldsymbol{a,0}}}\left( {\mathit{\boldsymbol{X}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{0}}}\right) \mathit{\boldsymbol{+}}{\mathit{\boldsymbol{\epsilon }}}^{\mathit{\boldsymbol{a,1}}}\left( {\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}\right) \right) }}_{O(\Delta t{h}_{e}^{2})}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(50)
|}


Which reads:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>{\hat{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}+</math><math>\Delta tO\left( {h}_{e}^{2}\right) ={\mathit{\boldsymbol{X}}}_{\lambda }^{1}+</math><math>O\left( \Delta {t}^{3}\right) +O({h}_{e}^{2}\Delta t)</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref3284803'></span>(51)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>{\hat{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}+</math><math>\Delta t\mathit{\boldsymbol{O}}\left( {h}_{e}^{2}\right) =</math><math>{\mathit{\boldsymbol{U}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}+</math><math>O\left( \Delta {t}^{3}\right) +O\mathit{\boldsymbol{(}}{h}_{e}^{2}\Delta t\mathit{\boldsymbol{)}}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|<span id='_Ref3284805'></span>(52)
|}


Eqs. <span id='cite-_Ref3284803'></span>[[#_Ref3284803|(51)]]-<span id='cite-_Ref3284805'></span>[[#_Ref3284805|(52)]] show that the error for the first time step is <math display="inline">O\left( \Delta {t}^{3}\right) +</math><math>O({h}_{e}^{2}\Delta t)</math>. It is straightforward to show that the resulting errors for the following time steps are those related to the accumulation, which reads:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{n+1}=</math><math>{\mathit{\boldsymbol{X}}}_{\lambda }^{n}+\left( n+1\right) O\left( \Delta {t}^{3}\right) +</math><math>\left( n+1\right) O\left( {h}_{e}^{2}\Delta t\right) ={\mathit{\boldsymbol{X}}}_{\lambda }^{n}+</math><math>{t}^{n+1}\left( O\left( \Delta {t}^{2}\right) +O\left( {h}_{e}^{2}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(53)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{\tilde{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{n+1}=</math><math>{\mathit{\boldsymbol{U}}}_{\mathit{\boldsymbol{\lambda }}}^{n}+</math><math>(n+1)O\left( \Delta {t}^{3}\right) +(n+1)O\mathit{\boldsymbol{(}}{h}_{e}^{2}\Delta t\mathit{\boldsymbol{)}}=</math><math>{\mathit{\boldsymbol{U}}}_{\lambda }^{n}+{t}^{n+1}\left( O\left( \Delta {t}^{2}\right) +\right. </math><math>\left. O\left( {h}_{e}^{2}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(54)
|}


resulting in a second order scheme in space and time.

It is obvious that, when solving the Navier-Stokes equations, the exact acceleration field is not known. On the contrary, it must be estimated. In this work, a FE Galerkin method is used to solve the pressure and viscous terms, which is a second order accurate method [<span id='cite-_Ref6318173'></span>[[#_Ref6318173|9]]]. Then, this error is also interpolated to the particles and has to be added to the interpolation error considered previously. However, this error does not degrades the rate of convergence of the method since it is second order, too.

As a conclusion, the accumulated error of the projected velocity is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\left\{ {\hat{\boldsymbol{u}}}_{b}^{n+1}\right\} ={P}_{\left\{ \lambda \right\} }^{n}\left[ \left\{ {\tilde{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{n+1}\right\} \right] =</math><math>{P}_{\left\{ \lambda \right\} }^{0}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\mathit{\boldsymbol{\lambda }}}^{n+1}+\right. \right. </math><math>\left. \left. {t}^{n}\left( O\left( \Delta {t}^{2}\right) +O\left( {h}_{e}^{2}\right) \right) \right\} \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(55)
|}


Then:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="margin:auto;width: 100%;"
|-
| <math>\left\{ {\hat{\mathit{\boldsymbol{u}}}}_{b}^{n+1}\right\} =\left\{ {\mathit{\boldsymbol{u}}}_{b}^{n+1}\right\} +</math><math>\underbrace{{\left\{ {\mathit{\boldsymbol{\epsilon }}}_{b}^{P,n+1}\right\} }}_{O\left( {h}_{e}^{2}\right) }+</math><math>\underbrace{{\left\{ {\mathit{\boldsymbol{e}}}_{b}^{n+1}\right\} }}_{{t}^{n+1}\left( O\left( \Delta {t}^{2}\right) +O\left( {h}_{e}^{2}\right) \right) }</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(56)
|}

===7. Convergence analyses===

This section presents different examples showing the rate of convergence of the SL-PFEM presented.

<span id='_Ref3285335'></span>

====1.4 Linear wave trajectories====

As it was pointed out in the introduction, for non-dominant convective flows the explicit streamlines computed at the beginning of the time step are not good approximants to pathlines and therefore small time steps are required to achieve accurate results [<span id='cite-_Ref532897341'></span>[[#_Ref532897341|5]]]. The flow field induced by a linear wave is a good example used to demonstrate this point. And, if information regarding the acceleration is introduced, the required time step can be increased (see <span id='cite-_Ref532897315'></span>[[#_Ref532897315|Figure 1]]).

In this section, the velocity Verlet algorithm is used to estimate the trajectory of a particle moving in the acceleration field induced by the velocities of a linear wave. The 2D velocity potential of an Airy wave with mean free surface located at <math display="inline">y=</math><math>0</math> is given by:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\phi (x,y,t)=\frac{Ag}{\omega }\frac{\mathrm{cosh}\,\left( K\left( H+y\right) \right) }{\mathrm{cosh}\,\left( KH\right) }\mathrm{sin}\,(Kx-</math><math>\omega t)</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(57)
|}


Where <math display="inline">\phi</math>  is the velocity potential, <math display="inline">A</math> is the wave amplitude, <math display="inline">K=</math><math>2\pi /L</math> is the wave number and <math display="inline">L</math> is the wave length, <math display="inline">H</math> is the water depth, <math display="inline">\omega =</math><math>2\pi /T</math> is the wave angular frequency and <math display="inline">T</math> is the wave period, <math display="inline">x</math> and <math display="inline">y</math> are the vertical and horizontal spatial coordinates, <math display="inline">t</math> is time. The velocity field is obtained as the gradient of the velocity potential <math display="inline">({u}_{x},{u}_{y})=</math><math>\left( {\phi }_{x},{\phi }_{y}\right)</math> . Then, the corresponding acceleration field is obtained via <math display="inline">\left( {a}_{x},{a}_{y}\right) =</math><math>{d}_{t}({u}_{x},{u}_{y})</math>. <span id='cite-_Ref532897315'></span>[[#_Ref532897315|Figure 1]] shows the velocity field induced by a linear wave. Using the velocity Verlet integrator we obtain the following equations for the position of the particle and its velocity components:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{X}_{\lambda }^{n+1}={X}_{\lambda }^{n}+\Delta t{\phi }_{x}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n})+</math><math>\frac{\Delta {t}^{2}}{2}\left( {\phi }_{xt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{x}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) {\phi }_{xx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}){\phi }_{xy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(58)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{Y}_{\lambda }^{n+1}={Y}_{\lambda }^{n}+\Delta t{\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n})+</math><math>\frac{\Delta {t}^{2}}{2}\left( {\phi }_{yt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{x}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) {\phi }_{yx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}){\phi }_{yy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(59)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{U}_{\lambda }^{n+1}-{U}_{\lambda }^{n}}{\Delta t}=\frac{1}{2}\left( {\phi }_{xt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{x}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) {\phi }_{xx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}){\phi }_{xy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) \right) +</math><math>\frac{1}{2}\left( {\phi }_{xt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) +\right. </math><math>\left. {\phi }_{x}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) {\phi }_{xx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}){\phi }_{xy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(60)
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>\frac{{V}_{\lambda }^{n+1}-{V}_{\lambda }^{n}}{\Delta t}=\frac{1}{2}\left( {\phi }_{yt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{y}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) {\phi }_{yx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}){\phi }_{yy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n},{t}^{n}\right) \right) +</math><math>\frac{1}{2}\left( {\phi }_{yt}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) +\right. </math><math>\left. {\phi }_{y}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) {\phi }_{yx}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) +\right. </math><math>\left. {\phi }_{y}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}){\phi }_{yy}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{n+1},{t}^{n+1}\right) \right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(61)
|}


The velocity Verlet has been used to estimate the particle trajectory due to a linear wave with <math display="inline">A=</math><math>0.01m</math>, <math display="inline">H=0.1m</math>, <math display="inline">L=1m</math>, and <math display="inline">T=</math><math>1.0726s</math> and with initial position <math display="inline">{x}_{0}=0.5m</math> and <math display="inline">{y}_{0}=</math><math>-0.02m</math>. The trajectory has been estimated for several time steps ( <math display="inline">\Delta t=</math><math>T/10</math>, <math display="inline">\, \Delta t=T/20</math>, <math display="inline">\, \Delta t=</math><math>T/40</math>, and <math display="inline">\, \Delta t=T/80</math>), using the analytical solution for the initial velocity. <span id='cite-_Ref531086640'></span>[[#_Ref531086640|Figure 3]] compares the analytical solution of the trajectory for one wave period with the numerical ones. Comparing the final position of the analytical solution and the numerical ones, an error estimate is obtained. We use this error to carry out a convergence analysis (<span id='cite-_Ref7188370'></span>[[#_Ref7188370|Table '''1''']]).

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image3.jpeg|600px]] </div>

<div id="_Ref531086640" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 3: Particle trajectory integration in the wave problem using the velocity Verlet algorithm. Black solid line represents the analytical solution. Red dots represent de numerical solution for </span> <math display="inline">\boldsymbol{\Delta }\mathit{\boldsymbol{t=T/5}}</math><span style="text-align: center; font-size: 75%;">,</span> <math display="inline">\boldsymbol{\, \Delta }\mathit{\boldsymbol{t=T/10}}</math><span style="text-align: center; font-size: 75%;">,</span> <math display="inline">\boldsymbol{\, \Delta }\mathit{\boldsymbol{t=T/20}}</math><span style="text-align: center; font-size: 75%;">,</span> <math display="inline">\boldsymbol{\, \Delta }\mathit{\boldsymbol{t=T/40}}</math><span style="text-align: center; font-size: 75%;">, and </span> <math display="inline">\boldsymbol{\, \Delta }\mathit{\boldsymbol{t=T/80}}</math><span style="text-align: center; font-size: 75%;">.</span></div>

<span id='_Ref7188370'></span><span id='_Ref7188358'></span><div id="_Ref7188365" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 1: Values and errors for the position and velocity of a particle at </span> <math display="inline">=</math><math>\mathit{\boldsymbol{T}}</math><span style="text-align: center; font-size: 75%;"> .</span></div>

{| style="width: 100%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;"|<math>x(T)</math>
|  style="border: 1pt solid black;text-align: center;"|<math>y(T)</math>
|  style="border: 1pt solid black;text-align: center;"|<math display="inline">Error</math> X
|  style="border: 1pt solid black;text-align: center;"|vx(T)
|  style="border: 1pt solid black;text-align: center;"|vy(T)
|  style="border: 1pt solid black;text-align: center;"|Error V
|-
|  style="border: 1pt solid black;text-align: center;"|Analytical
|  style="border: 1pt solid black;text-align: center;"|0.508308
|  style="border: 1pt solid black;text-align: center;"|-0.020000
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.09850702
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.00238937
|  style="border: 1pt solid black;text-align: center;"|
|-
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=T/10</math>
|  style="border: 1pt solid black;text-align: center;"|0.503153
|  style="border: 1pt solid black;text-align: center;"|-0.020397
|  style="border: 1pt solid black;text-align: center;"|0.00517037
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.09850997
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.00178868
|  style="border: 1pt solid black;text-align: center;"|0.00060070
|-
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=T/20</math>
|  style="border: 1pt solid black;text-align: center;"|0.507039
|  style="border: 1pt solid black;text-align: center;"|-0.020091
|  style="border: 1pt solid black;text-align: center;"|0.00127241
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.09851727
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.00222848
|  style="border: 1pt solid black;text-align: center;"|0.00016122
|-
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=T/40</math>
|  style="border: 1pt solid black;text-align: center;"|0.507992
|  style="border: 1pt solid black;text-align: center;"|-0.020015
|  style="border: 1pt solid black;text-align: center;"|0.00031649
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.09851019
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.00234849
|  style="border: 1pt solid black;text-align: center;"|0.00004100
|-
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=T/80</math>
|  style="border: 1pt solid black;text-align: center;"|0.508229
|  style="border: 1pt solid black;text-align: center;"|-0.019997
|  style="border: 1pt solid black;text-align: center;"|7.9029E-05
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.09850785
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|-0.00237911
|  style="border: 1pt solid black;text-align: center;"|0.00001029
|}


Now, we will integrate Eqs. <span id='cite-_Ref7798987'></span>[[#_Ref7798987|(1)]]-<span id='cite-_Ref7798989'></span>[[#_Ref7798989|(2)]] using the semi-Lagrangian approach given the acceleration field induced by a linear wave. The acceleration and velocities at the mesh nodes ( <math display="inline">{\mathit{\boldsymbol{a}}}_{h}^{0}</math> and <math display="inline">{\mathit{\boldsymbol{u}}}_{h}^{0}</math>)  are initialized with the analytical values. The particle’s initial velocities and accelerations are interpolated from the mesh initial values ( <math display="inline">{\mathit{\boldsymbol{U}}}_{\lambda }^{0}=</math><math>{\mathit{\boldsymbol{u}}}_{h}^{0}\left( {\mathit{\boldsymbol{X}}}_{\lambda }^{0}\right)</math>  and <math display="inline">{\mathit{\boldsymbol{A}}}_{\lambda }^{0}=</math><math>{\mathit{\boldsymbol{a}}}_{h}^{0}({\mathit{\boldsymbol{X}}}_{\lambda }^{0})</math>) . Afterwards, during the simulation the acceleration field is imposed on the mesh nodes and interpolated at the particles, while the velocities are obtained via the semi-Lagrangian approach. The algorithm reads as follows:

Step 1: Move particles from <math display="inline">t=n\Delta t</math> to <math display="inline">t=</math><math>\left( n+1\right) \Delta t</math>:
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>{\mathit{\boldsymbol{X}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{X}}}_{\lambda }^{n}+</math><math>\Delta t{\mathit{\boldsymbol{u}}}_{h}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})+</math><math>\frac{\Delta {t}^{2}}{2}{\mathit{\boldsymbol{a}}}_{h}^{n}({\mathit{\boldsymbol{X}}}_{\lambda }^{n})</math>
|}


Step 2: Interpolate the particle’s acceleration at the new position:
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>{\mathit{\boldsymbol{A}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{a}}}_{h}^{n+1}({\mathit{\boldsymbol{X}}}_{\lambda }^{n+1})</math>
|}


Step 3: The new particle´s velocity is updated:
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>{\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}={\mathit{\boldsymbol{U}}}_{\lambda }^{n}+</math><math>\frac{1}{2}\Delta t\left( {\mathit{\boldsymbol{A}}}_{\lambda }^{n}+{\mathit{\boldsymbol{A}}}_{\lambda }^{n+1}\right)</math> 
|}


Step 4: The particle´s velocities are projected onto the mesh to obtain the new velocity field on the background mesh:
{| class="formulaSCP" style="width: 100%; text-align: center;" 
|-
| <math>{\mathit{\boldsymbol{u}}}_{h}^{n+1}{=P}_{\left\{ \lambda \right\} }^{n+1}\left[ \left\{ {\mathit{\boldsymbol{U}}}_{\lambda }^{n+1}\right\} \right]</math> 
|}


In order to achieve second order convergence, each step must be at least second order in time and space. Step 1 has already been demonstrated that is second order. Step 2 interpolates from the mesh values to the particles using the FE linear shape functions. This interpolation is second order as demonstrated in [<span id='cite-_Ref532985341'></span>[[#_Ref532985341|3]]]. Step 3 integrates the velocity in time with a Crank-Nicolson scheme, which is known to be second order in time. And finally, Step 4 is the projection step using the global least square projector which, in the previous section has been demonstrated to be second order in space and time.

The time step has been setup such that <math display="inline">h/\Delta t=4.661</math>, being <math display="inline">h</math> the characteristic mesh size. <span id='cite-_Ref531605890'></span>[[#_Ref531605890|Table 2]] provides the values of the mean quadratic error (MQE) of the velocity on the mesh and <span id='cite-_Ref531778866'></span>[[#_Ref531778866|Figure 4]] shows that the rate of convergence is of second order, as expected.

<div id="_Ref531605890" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 2: Mean quadratic error of the velocity on the mesh at </span> <math display="inline">\boldsymbol{t=T}</math></div>

{| style="width: 56%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;"|<math display="inline">RMSE=</math><math>\sqrt{\sum _{i=1}^{N}{\left( {\mathit{\boldsymbol{u}}}_{h}\left( {\mathit{\boldsymbol{x}}}_{i}\right) -{\mathit{\boldsymbol{u}}}_{a}\left( {\mathit{\boldsymbol{x}}}_{i}\right) \right) }^{2}/N}</math> 
|-
|  style="border: 1pt solid black;text-align: center;"|<math>h=H\boldsymbol{/}2</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|1.226E-03
|-
|  style="border: 1pt solid black;text-align: center;"|<math>h=H\boldsymbol{/}4</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|2.345E-04
|-
|  style="border: 1pt solid black;text-align: center;"|<math>h=H\boldsymbol{/}6</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|1.176E-04
|-
|  style="border: 1pt solid black;text-align: center;"|<math>h=H\boldsymbol{/}8</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|6.951E-05
|-
|  style="border: 1pt solid black;text-align: center;"|<math>h=H/10</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: bottom;"|4.742E-05
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image4.png|462px]] </div>

<div id="_Ref531778866" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 4: Convergence analysis for the velocity on the mesh for a linear wave using the semi-Lagrangian approach.</span></div>

====1.5 Lid driven cavity flow.====

The particulars for this case study are given in <span id='cite-_Ref6225834'></span>[[#_Ref6225834|Table 3]], and <span id='cite-_Ref6302688'></span>[[#_Ref6302688|Figure 5]] shows the discretization used and the boundary conditions applied. An average of 2 iterations is only required for convergence.

Although the steady state solution is achieved after one thousand time steps, the simulation time has been set to <math display="inline">T=</math><math>{10}^{5}\Delta t</math> to verify whether the steady state solution shows symptoms of error accumulation.

<span id='_Ref6225974'></span>The results obtained with the present second-order Verlet-SLPFEM are compared to those obtained using a high resolution spectral method in [<span id='cite-11'></span>[[#11|11]]]. Also results are shown for a first-order  SL-PFEM based on the Backward Euler integrator for the velocity equation. <span id='cite-_Ref6225882'></span>[[#_Ref6225882|Figure 6]] compares the streamlines and <span id='cite-_Ref6225884'></span>[[#_Ref6225884|Figure 7]] compares the horizontal velocity and pressure at the mid-section. While the second-order SL-PFEM based on the Verlet algorithm provides very similar results to those of [<span id='cite-_Ref6225974'></span>[[#_Ref6225974|11]]], the first-order SL-PFEM based on the Euler scheme produces a highly diffusive solution on the pressure and velocity. Despite being a steady-state problem the first-order method is incapable of achieving the right solution with the current time step. It would require a significant reduction of the time step, which leads to a significant increase of CPU time. On the other hand, the second-order scheme produces an accurate solution for the current time step.

<span id='cite-_Ref6225885'></span>[[#_Ref6225885|Figure 8]] shows the evolution of the horizontal velocity at the node located at (x,y)=(0.5,0.175). It is observed that an average of 0.392 meters per second is obtained with a variance of 0.001 m/s approximately. No symptoms of error accumulation are observed.

<div id="_Ref6225834" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 3: particular of the lid driven cavity flow</span></div>

{| style="width: 57%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;vertical-align: top;"|Lid velocity <math display="inline">V</math>
|  style="border: 1pt solid black;vertical-align: top;"|<math>-1m/s</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Reynolds number <math display="inline">Re</math>
|  style="border: 1pt solid black;vertical-align: top;"|<math>1000</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Domain size
|  style="border: 1pt solid black;vertical-align: top;"|<math>1m\, x\, 1m</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Domain discretization
|  style="border: 1pt solid black;vertical-align: top;"|<math>80x80</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Number of triangle elements
|  style="border: 1pt solid black;vertical-align: top;"|<math>25600</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Number of nodes
|  style="border: 1pt solid black;vertical-align: top;"|<math>12961</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Average number of particles per element
|  style="border: 1pt solid black;vertical-align: top;"|<math>3</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Mesh size <math display="inline">\Delta x=</math><math>\Delta y</math>
|  style="border: 1pt solid black;vertical-align: top;"|<math>0.0125</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Time step <math display="inline">\Delta t</math>
|  style="border: 1pt solid black;vertical-align: top;"|<math>0.1</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Courant number
|  style="border: 1pt solid black;vertical-align: top;"|<math>8</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Simulation time
|  style="border: 1pt solid black;vertical-align: top;"|<math>{10}^{5}\Delta t</math>
|-
|  style="border: 1pt solid black;vertical-align: top;"|Sampling time
|  style="border: 1pt solid black;vertical-align: top;"|<math>{10}^{2}\Delta t</math>
|}


{| style="text-align: center;width: 63%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image5.jpeg|366px]] 
|}


<div id="_Ref6302688" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 5: Structured mesh for lid driven cavity flow and boundary conditions.</span></div>

{| style="width: 74%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="text-align: center;vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image6.jpeg|216px]]  [[Image:Draft_Garcia-Espinosa_450194514-image7.jpeg|216px]] 
|-
|  style="text-align: center;vertical-align: top;"|
|-
|  style="text-align: center;vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image8.jpeg|216px]]  [[Image:Draft_Garcia-Espinosa_450194514-image9.jpeg|216px]] 
|-
|  style="vertical-align: top;"|
|}


<span id='_Ref6225882'></span><div id="_Ref6225876" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 6: Streamlines comparison. Up: Second order Verlet-SLPFEM. Down: First Order Euler SL-PFEM</span></div>

{| style="text-align: center;width: 57%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image10.png|336px]]

[[Image:Draft_Garcia-Espinosa_450194514-image11.png|336px]] 
|}


<span id='_Ref6225884'></span><div id="_Ref6225877" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 7: Horizontal velocity and pressure profiles at the mid-section x=0.5. Red dots: spectral method [</span><span id='cite-_Ref6225974'></span>[[#_Ref6225974|<span style="text-align: center; font-size: 75%;">11</span>]]<span style="text-align: center; font-size: 75%;">]. Solid line: Second Order Verlet-SLPFEM. Dash line: First Order Euler-SLPFEM.</span></div>

{| style="text-align: center;width: 58%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image12.png|336px]] 
|}


<div id="_Ref6225885" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 8: Horizontal velocity evolution at (x,y)=(0.5,0.175) for the second-order Verlet -SLPFEM</span></div>

====1.6 Two-dimensional Taylor-Green vortex.====

The two-dimensional Taylor-Green vortex problem is used to evaluate the convergence order of the numerical strategy solving an incompressible flow problem for which an analytical solution is known. The analytical solution is given by:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{u}_{x}\left( x,y,t\right) =-\mathrm{sin}\,\left( x\right) \mathrm{cos}\,\left( y\right) {e}^{-2\nu t}</math>

<math>{u}_{y}\left( x,y,t\right) =\mathrm{cos}\,\left( x\right) \mathrm{sin}\,\left( y\right) {e}^{-2\nu t}</math>

<math>P\left( x,y,t\right) =\frac{1}{4}\left[ \mathrm{cos}\,\left( 2x\right) +\mathrm{cos}\,\left( 2y\right) \right] {e}^{-4\nu t}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(62)
|}


where <math display="inline">\nu =0.001{m}^{2}/s</math> is the kinematic viscosity. The problem is solved using a square domain <math display="inline">[0,\pi ]x[0,\pi ]</math> and setting slip boundary conditions for the velocity. Also, to prevent the pressure to be undetermined, the following condition is imposed at <math display="inline">\left( x,y\right) =</math><math>\left( 0,0\right)</math> :

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>P\left( 0,0,t\right) =\frac{1}{2}{e}^{-4\nu t}</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(63)
|}


The velocity and pressure field are initialized at <math display="inline">t=0</math> with the analytical solution. The simulation is carried out until a final time of <math display="inline">t=</math><math>10s</math> is reached. <span id='cite-_Ref532827353'></span>[[#_Ref532827353|Table 4]] provides the values of the mesh size, time step, and the number of time steps computed. <span id='cite-_Ref532983488'></span>[[#_Ref532983488|Figure 9]] shows the mesh used for case 2.

<div id="_Ref532827353" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 4: Numerical parameters use in the 2D Taylor-Green Vortex</span></div>

{| style="width: 52%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|Case
|  style="border: 1pt solid black;text-align: center;"|Mesh size
|  style="border: 1pt solid black;text-align: center;"|Time step
|  style="border: 1pt solid black;text-align: center;"|Number of

time steps
|-
|  style="border: 1pt solid black;text-align: center;"|1
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /8</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.2s</math>
|  style="border: 1pt solid black;text-align: center;"|50
|-
|  style="border: 1pt solid black;text-align: center;"|2
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /16</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.1s</math>
|  style="border: 1pt solid black;text-align: center;"|100
|-
|  style="border: 1pt solid black;text-align: center;"|3
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /24</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.0667s</math>
|  style="border: 1pt solid black;text-align: center;"|150
|-
|  style="border: 1pt solid black;text-align: center;"|4
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /32</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.05s</math>
|  style="border: 1pt solid black;text-align: center;"|200
|-
|  style="border: 1pt solid black;text-align: center;"|5
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /48</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.0333s</math>
|  style="border: 1pt solid black;text-align: center;"|300
|-
|  style="border: 1pt solid black;text-align: center;"|6
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /64</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.025s</math>
|  style="border: 1pt solid black;text-align: center;"|400
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image13.png|282px]] </div>

<div id="_Ref532983488" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 9: Structured mesh with </span> <math display="inline">\boldsymbol{\Delta }\mathit{\boldsymbol{x}}=</math><math>\mathit{\boldsymbol{\pi }}\boldsymbol{/16}</math><span style="text-align: center; font-size: 75%;">.</span></div>

<span id='cite-_Ref532890800'></span>[[#_Ref532890800|Figure 10]] shows the convergence of the proposed algorithm. It can be seen that both, pressure and velocity converge with second order convergence rate.

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image14.png|600px]] </div>

<div id="_Ref532890800" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 10: Convergence analysis for the 2D Taylor-Green Vortex. Black dashed line represents first order convergence. Black dotted line represents second order convergence. Green: pressure errors. Red: velocity errors.</span></div>

====1.7 Two-dimensional Steady Taylor-Green vortex.====

A modified version of the two-dimensional Taylor-Green vortex problem is used to evaluate the convergence order of the numerical strategy solving an incompressible flow problem for which an analytical solution is known. In this version, a mass force is introduced to replace the energy dissipated by the viscous terms, so that a steady solution is obtained. The mass force to be introduced is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{f}_{x}\left( x,y\right) =-2\nu \mathrm{sin}\,\left( x\right) \mathrm{cos}\,\left( y\right)</math>

<math>{f}_{y}\left( x,y\right) =2\nu \mathrm{cos}\,\left( x\right) \mathrm{sin}\,\left( y\right)</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(64)
|}


and the analytical solution to the problem is:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>{u}_{x}\left( x,y\right) =-\mathrm{sin}\,\left( x\right) \mathrm{cos}\,\left( y\right)</math>

<math>{u}_{y}\left( x,y\right) =\mathrm{cos}\,\left( x\right) \mathrm{sin}\,\left( y\right)</math>

<math>P\left( x,y\right) =\frac{1}{4}\left[ \mathrm{cos}\,\left( 2x\right) +\mathrm{cos}\,\left( 2y\right) \right]</math> 
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(65)
|}


The problem is solved with <math display="inline">\nu =0.01{m}^{2}/s</math> ( <math display="inline">Re=</math><math>314</math>). The fluid domain is a square of <math display="inline">[0,\pi ]x[0,\pi ]</math>. Slip boundary conditions are used for the velocity. Also, to prevent the pressure to be undetermined, the following Dirichlet condition is imposed:

{| class="formulaSCP" style="width: 100%;width: 100%;text-align: center;" 
|-
| 
{| style="text-align: center; margin:auto;width: 100%;"
|-
| <math>P\left( 0,0\right) =0.5</math>
|}
|  style="width: 5px;text-align: right;white-space: nowrap;"|(66)
|}


The velocity and pressure field are initialized with the analytical solution. And average of 2 iterations is only required for convergence.  The simulation is carried out for <math display="inline">t=</math><math>400s</math>. Based on the maximum velocity ( <math display="inline">{u}_{max}=1</math>), the Courant number used is <math display="inline">Cr=</math><math>\Delta t/\Delta x=2.04</math>. The case is initiated with 3 particles per element, and particles are removed when it exceeds 6 particles per element. At the end of the simulation the average number of particles per element is slightly lower than 3.

<span id='cite-_Ref7188377'></span>[[#_Ref7188377|Table 5]] provides the values of the mesh size, time step, and the number of time steps computed. <span id='cite-_Ref7188407'></span>[[#_Ref7188407|Figure 11]] shows an intermediate mesh used in the analysis. <span id='cite-_Ref7188426'></span>[[#_Ref7188426|Figure 12]] how the pressure on the mesh and the particle’s velocity obtained for case 64. Both results are quite close to the analytical solution. For this case, the root mean square error is in the order of 0.003 for both, velocity and pressure.

<div id="_Ref7188377" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 5: Numerical parameters use in the 2D Steady Taylor-Green Vortex</span></div>

{| style="width: 55%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|Case
|  style="border: 1pt solid black;text-align: center;"|Mesh size
|  style="border: 1pt solid black;text-align: center;"|Time step
|  style="border: 1pt solid black;text-align: center;"|Number of

time steps
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|16
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /16</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.4s</math>
|  style="border: 1pt solid black;text-align: center;"|1000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|24
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /24</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.26667s</math>
|  style="border: 1pt solid black;text-align: center;"|1500
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|32
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /32</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.2s</math>
|  style="border: 1pt solid black;text-align: center;"|2000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|48
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /48</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.13333s</math>
|  style="border: 1pt solid black;text-align: center;"|3000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|64
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /64</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.1s</math>
|  style="border: 1pt solid black;text-align: center;"|4000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|96
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /96</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.06667s</math>
|  style="border: 1pt solid black;text-align: center;"|6000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|128
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /128</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.05s</math>
|  style="border: 1pt solid black;text-align: center;"|8000
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|192
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x=\pi /192</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t=0.03333s</math>
|  style="border: 1pt solid black;text-align: center;"|12000
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image15.png|366px]] </div>

<div id="_Ref7188407" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 11: Structured mesh with </span> <math display="inline">\boldsymbol{\Delta }\mathit{\boldsymbol{x}}=</math><math>\mathit{\boldsymbol{\pi }}\boldsymbol{/64}</math><span style="text-align: center; font-size: 75%;"> and boundary conditions applied</span></div>

{| style="width: 87%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="text-align: center;vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image16.jpeg|252px]]  [[Image:Draft_Garcia-Espinosa_450194514-image17.jpeg|252px]] 
|}


<div id="_Ref7188426" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 12: Case 64 results. Pressure (left) and velocity (right)</span></div>

<span id='cite-_Ref7188443'></span>[[#_Ref7188443|Figure 13]] shows the evolution of the errors along the simulation. It is observed that the errors achieved a constant value.  Hence no symptoms of error accumulation over time is appreciated. <span id='cite-_Ref7188509'></span>[[#_Ref7188509|Figure 14]] provides the rate of convergence obtained using an average of three particles per element. The root mean square error (RMSE) has been used for the analysis. <span id='cite-_Ref7188526'></span>[[#_Ref7188526|Table 6]] shows the convergence of the proposed algorithm for pressure and velocity.

{| style="width: 44%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;vertical-align: top;width: 0%;"|[[Image:Draft_Garcia-Espinosa_450194514-image18.png|258px]]

[[Image:Draft_Garcia-Espinosa_450194514-image19.png|258px]] 
|}


<div id="_Ref7188443" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 13: Error evolution along time.</span></div>

<div id="_Ref7188526" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 6: Rate of convergence. Root mean square error (RMSE) and maximum error (MaxE)</span></div>

{| style="margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  colspan='2'  style="border: 1pt solid black;text-align: center;"|Convergence Rate
|-
|  style="border: 1pt solid black;text-align: center;"|Velocity
|  style="border: 1pt solid black;text-align: center;"|Pressure
|-
|  style="border: 1pt solid black;text-align: center;"|2.59
|  style="border: 1pt solid black;text-align: center;"|2.66
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image20.jpeg|390px]] </div>

<div id="_Ref7188509" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 14: Convergence analysis for the 2D Taylor-Green Vortex. Black dashed line represents first order convergence. Black dotted line represents second order convergence. Blue: pressure errors. Red: velocity errors.</span></div>

While in an Eulerian method all dependent variables remain unchanged once the steady state solution is reached, this is not the case in the SL-PFEM method. From the point of view of the particles, the flow is never steady since the intrinsic variables they transport are continually changing.

If a SL-PFEM is not properly tailored, then the error accumulation will make the method to be incapable of computing a steady-state solution. This is because the accumulation of errors along the simulation time will eventually cause a numerical blow up. The present SL-PFEM formulation is capable of computing steady state solutions without errors concerns.

====1.8 Three-dimensional flow past a cylinder.====

The three dimensional flow around a cylinder is simulated to assess the suitability of the present method for solving three dimensional problems. The particulars of the problem are given in <span id='cite-_Ref535835118'></span>[[#_Ref535835118|Table 7]]. An average of 2-3 iterations are only required for convergence.

<span id='_Ref535836174'></span>For a Reynolds number of 200 and an aspect ratio of the cylinder height (H) to the cylinder diameter (D) of H/D=1, the flow behaves as a two dimensional flow [<span id='cite-12'></span>[[#12|12]]].  A three dimensional structured mesh with five levels of refinement has been used, keeping the Courant number constant. <span id='cite-_Ref535840176'></span>[[#_Ref535840176|Figure 15]] shows the computational domain used, and <span id='cite-_Ref535840212'></span>[[#_Ref535840212|Figure 16]] shows the corresponding meshes used in this analysis. The particulars for each mesh are provided in <span id='cite-_Ref12462928'></span>[[#_Ref12462928|Table 8]], as well as the total number of time steps simulated. In average, two particles per elements are used.

<span id='_Ref535835118'></span><div id="_Ref535835114" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 7: Flow past a cylinder particulars</span></div>

{| style="width: 35%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Cylinder\, diameter\, (D)</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<math>1m</math>
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Cylinder\, height\, (H)</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<math>1m</math>
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Inlet\, velocity\, (V)</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<math>1m/s</math>
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Viscosity\, (\nu )</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<math>0.005</math>
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Re</math>
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|<math>200</math>
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image21.jpeg|384px]] </div>

<div id="_Ref535840176" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 15: Computational domain used to generate the structured mesh.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image22.jpeg|600px]] </div>

<div id="_Ref535840212" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 16: Computational mesh refinement for the 3D flow past a cylinder.</span></div>


<div id="_Ref12462928" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 8: Numerical data used in the 3D flow past a cylinder</span></div>

{| style="width: 57%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Mesh</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta x</math>
|  style="border: 1pt solid black;text-align: center;"|<math>\Delta t</math>
|  style="border: 1pt solid black;text-align: center;"|<math>N\, tetras</math>
|  style="border: 1pt solid black;text-align: center;"|<math>N\, Nodes</math>
|  style="border: 1pt solid black;text-align: center;"|<math>N\, \Delta t</math>
|-
|  style="border: 1pt solid black;text-align: center;"|1
|  style="border: 1pt solid black;text-align: center;"|0.1333
|  style="border: 1pt solid black;text-align: center;"|0.0667
|  style="border: 1pt solid black;text-align: center;"|62208
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|14487
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|1500
|-
|  style="border: 1pt solid black;text-align: center;"|2
|  style="border: 1pt solid black;text-align: center;"|0.1
|  style="border: 1pt solid black;text-align: center;"|0.05
|  style="border: 1pt solid black;text-align: center;"|147456
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|33412
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|2000
|-
|  style="border: 1pt solid black;text-align: center;"|3
|  style="border: 1pt solid black;text-align: center;"|0.08
|  style="border: 1pt solid black;text-align: center;"|0.04
|  style="border: 1pt solid black;text-align: center;"|288000
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|64185
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|2500
|-
|  style="border: 1pt solid black;text-align: center;"|4
|  style="border: 1pt solid black;text-align: center;"|0.0667
|  style="border: 1pt solid black;text-align: center;"|0.0333
|  style="border: 1pt solid black;text-align: center;"|497664
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|109686
|  style="border: 1pt solid black;text-align: center;vertical-align: top;"|3000
|}


<span id='cite-_Ref535840306'></span>[[#_Ref535840306|Table 9]] provides the Strouhal numbers obtained in each case study. It is observed how the lift force amplitude and the Strouhal number converge towards values of 0.27 and 0.196, respectively. <span id='cite-_Ref535854548'></span>[[#_Ref535854548|Figure 17]] plots the logarithmic errors for the Strouhal number versus the logarithmic mesh size for each case study. Comparing with the second order reference line, it is observed that the convergence rate of the Strouhal number is close to second order. <span id='cite-_Ref535854807'></span>[[#_Ref535854807|Figure 18]] and <span id='cite-_Ref536173038'></span>[[#_Ref536173038|Figure 19]] provides a snapshot of the pressure field, particles velocities, streamlines and vorticity for Case 2.

<div id="_Ref535840306" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Table 9: Numerical results for the 3D flow past a cylinder</span></div>

{| style="width: 38%;margin: 1em auto 0.1em auto;border-collapse: collapse;" 
|-
|  style="border: 1pt solid black;text-align: center;"|<math>Case</math>
|  style="border: 1pt solid black;text-align: center;"|<math>St</math>
|  style="border: 1pt solid black;text-align: center;"|<math>Error\, St</math>
|  style="border: 1pt solid black;text-align: center;"|<math>{C}_{L}</math>
|-
|  style="border: 1pt solid black;text-align: center;"|1
|  style="border: 1pt solid black;text-align: center;"|0.133
|  style="border: 1pt solid black;text-align: center;"|0.0630
|  style="border: 1pt solid black;text-align: center;"|0.177
|-
|  style="border: 1pt solid black;text-align: center;"|2
|  style="border: 1pt solid black;text-align: center;"|0.162
|  style="border: 1pt solid black;text-align: center;"|0.0339
|  style="border: 1pt solid black;text-align: center;"|0.255
|-
|  style="border: 1pt solid black;text-align: center;"|3
|  style="border: 1pt solid black;text-align: center;"|0.178
|  style="border: 1pt solid black;text-align: center;"|0.0184
|  style="border: 1pt solid black;text-align: center;"|0.246
|-
|  style="border: 1pt solid black;text-align: center;"|4
|  style="border: 1pt solid black;text-align: center;"|0.186
|  style="border: 1pt solid black;text-align: center;"|0.0101
|  style="border: 1pt solid black;text-align: center;"|0.258
|-
|  style="border: 1pt solid black;text-align: center;"|Ref. [<span id='cite-_Ref535836174'></span>[[#_Ref535836174|12]]]
|  style="border: 1pt solid black;text-align: center;"|0.196
|  style="border: 1pt solid black;text-align: center;"|
|  style="border: 1pt solid black;text-align: center;"|
|}


<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image23.jpg|408px]] </div>

<div id="_Ref535854548" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 17: Rate of convergence for the Strouhal number. Circles represent the log (error) for each case study. Solid line represents the trend of the convergence. Dot line represents second order convergence.</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 
{|
|-
| [[Image:Draft_Garcia-Espinosa_450194514-image24.jpeg|264px]]
| [[Image:Draft_Garcia-Espinosa_450194514-image25.jpeg|center|264px]]
|}
</div>

<div id="_Ref535854807" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 18: Case 2 snapshots: top view of pressure map on mesh (left) and particles velocities (right).</span></div>

<div class="center" style="width: auto; margin-left: auto; margin-right: auto;">
 [[Image:Draft_Garcia-Espinosa_450194514-image26.jpeg|420px]] </div>

<div id="_Ref536173038" class="center" style="width: auto; margin-left: auto; margin-right: auto;">
<span style="text-align: center; font-size: 75%;">Figure 19: Case 2 snapshots: oblique view of vorticity on particles.</span></div>

===Conclusions===

In this paper, a second order semi-Lagrangian (SL) method for solving the incompressible Navier-Stokes equations has been presented within the framework of the Particle Finite Element Method (PFEM). The method is based on the velocity Verlet integrator, which offers an explicit formulation for particle tracking with good numerical properties. Therefore, particles only need to be transported once per time step, which makes this integrator very attractive for the semi-Lagrangian PFEM (SL-PFEM). The algorithm is completed with a predictor-multicorrector scheme for integrating the velocity correction using the Finite Element Method. A second order global least square projector is used to make sure that the coherence condition is fulfilled, which avoids the need of iterating and projecting more than once per time step.  The proposed method has been tested for different types of flows in two and three dimensions, and the second order rate of convergence has been achieved. Moreover, it has been shown that there are no symptoms of error accumulation along time, which is a major concern for semi-Lagrangian schemes because it could compromise the rate of convergence and its capability to compute steady-state problems.

===Acknowledgements===

This work relates to Department of the Navy Grant N62909-16-1-2236 issued by Office of Naval Research Global. The United States Government has a royalty-free licence throughout the world in all copyrightable material contained herein.

===Conflict of interest===

On behalf of all authors, the corresponding author states that there is no conflict of interest

===References===

<div id="1"></div>
[[#cite-1|[1]]] S.R. Idelsohn, J. Marti, P. Becker, E. Oñate: Analysis of multifluid flows with large time steps using the particle finite element method. International Journal for Numerical Methods in Fluids 75(9), 621–644 (2014). DOI 10.1002/fld.3908.

<div id="2"></div>
[[#cite-2|[2]]] S. R. Idelsohn, N. Nigro, A. Limache, E. Oñate: Large time-step explicit integration method for solving problems with dominant convection. Computer Methods in Applied Mechanics and Engineering 217-220, 168–185 (2012). DOI 10.1016/j.cma.2011.12. 008.

<div id="3"></div>
[[#cite-3|[3]]] P. Becker, S. R. Idelsohn, E. Oñate. A unified monolithic approach for multi-fluid flows and fluidâASstructure interaction using the Particle Finite Element Method with fixed mesh. Computational Mechanics (2014). DOI 10.1007/s00466-014-1107-0

<div id="4"></div>
[[#cite-4|[4]]] J. M. Gimenez, H. J. Aguerrea, S. R. Idelsohnd, N. M. Nigro. A space-time second-order particle-based method to solve ow problems on arbitrary meshes. Accepted for publication in Journal of Computational Physics. DOI 10.1016/j.jcp.2018.11.034

<div id="5"></div>
[[#cite-5|[5]]] P. Nadukandi, B. Servan-Camas, P. A. Becker, and J. Garcia-Espinosa. Seakeeping with the semi-Lagrangian Particle Finite Element Method. <span style="text-align: center; font-size: 75%;">Comp. Part. Mech. (2017) 4: 321. DOI 10.1007/s40571-016-0127-2</span>

<div id="6"></div>
[[#cite-6|[6]]] L. Verlet. Computer "Experiments" on Classical Fluids. I. Thermodynamical properties of Lenard-Jones Molecules. Phys. Rev. E. 59(5), 98-103. (1967)

<div id="7"></div>
[[#cite-7|[7]]] W. C. Swope, H. C. Andersen, P. H. Berens, and K. R. Wilson. A computer simulation method for the calculation of equilibrium constants for the formation of physical clusters of molecules: Application to small water clusters. J. Chem. Phys. 76, 637 (1982). DOI 10.1063/1.442716.

<div id="8"></div>
[[#cite-8|[8]]] E. Hairer, C. Lubich and G. Wanner (2003). Geometric numerical integration illustrated by the Störmer–Verlet method. Acta Numerica, 12, pp 399-450. DOI 10.1017/S0962492902000144

<div id="9"></div>
[[#cite-9|[9]]] J. Garcia-Espinosa, A. Valls, and E. Oñate. ODDLS: A new unstructured mesh ﬁnite element method for the analysis of free surface ﬂow problems Int. J. Numer. Meth. Engng 76(9) 1297-1470 (2008).

<div id="10"></div>
[[#cite-10|[10]]] E. Oñate , J. García-Espinosa, S.R. Idelsohn, F. Del Pin. Finite calculus formulations for ﬁnite element analysis of incompressible ﬂows. Eulerian, ALE and Lagrangian approaches. Comput. Methods Appl. Mech. Engrg. 195 3001–3037 (2006).

<div id="11"></div>
[[#cite-11|[11]]] O. Botella and R. Peyret. Benchmark spectral results on the lid-driven cavity flow. Computers & Fluids Vol. 27, No. 4, pp. 421-433, 1998.

<div id="12"></div>
[[#cite-12|[12]]] H. Jiang and L. Cheng. Strouhal Reynolds number relationship for flow past a circular cylinder. J. Fluid Mech. (2017), vol. 832, pp. 170-188.


The [https://link.springer.com/article/10.1007%2Fs40571-019-00258-9 final publication] is available at link.springer.com.