​

​

==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|384px]] </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="text-align: center; margin:auto;width: 100%;"

|-

| <math>{\tilde{\mathit{\boldsymbol{X}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>+++</math>

|}

|  style="width: 5px;text-align: right;white-space: nowrap;"|(49)

|-

| 

{| style="text-align: center; margin:auto;width: 100%;"

|-

| <math>{\tilde{\mathit{\boldsymbol{U}}}}_{\mathit{\boldsymbol{\lambda }}}^{\mathit{\boldsymbol{1}}}=</math><math>++</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="text-align: center; margin:auto;width: 100%;"

|-

| <math>\left\{ {\hat{\mathit{\boldsymbol{u}}}}_{b}^{n+1}\right\} =\left\{ {\mathit{\boldsymbol{u}}}_{b}^{n+1}\right\} +</math><math>+</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''']]).

​

[[Image:Draft_Garcia-Espinosa_450194514-image3.jpeg|600px]]

​

<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>