Latest revision as of 19:31, 6 August 2019

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.

The final publication is available at link.springer.com.

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

Until [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 [1, 2, 3, 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 [5], where a linear wave problem was used as an example. Figure 1 has been extracted from [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 [6, 7, 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 [8].

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 ${\textstyle {\mathit {\boldsymbol {u}}}({\mathit {\boldsymbol {x}}},t)}$ be a velocity field, let ${\textstyle \left\{\lambda \right\}}$ be a set of particles each of them identified with a label ${\textstyle \lambda }$ . The position occupied by each particle at time ${\textstyle t}$ is denoted by ${\textstyle {\mathit {\boldsymbol {X}}}_{\lambda }{\mathit {\boldsymbol {(}}}t{\mathit {\boldsymbol {)}}}}$, and the velocity of the particle by ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }{\mathit {\boldsymbol {(}}}t{\mathit {\boldsymbol {)}}}}$. If the particle velocity is given by the velocity field ${\textstyle {\mathit {\boldsymbol {u}}}}$, then ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }\left(t\right){\mathit {\boldsymbol {=\,}}}{\mathit {\boldsymbol {u}}}({\mathit {\boldsymbol {X}}}_{\lambda }{\mathit {\boldsymbol {(}}}t{\mathit {\boldsymbol {)}}},t)}$, and the particles’s trajectories are the pathlines defined by the velocity field ${\textstyle {\mathit {\boldsymbol {u}}}}$. The particle’s equations of motion are:

where the particle’s acceleration ${\textstyle {\mathit {\boldsymbol {A}}}_{\lambda }\left(t\right)}$ is the rate of change of ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }{\mathit {\boldsymbol {(}}}t{\mathit {\boldsymbol {)}}}}$. An acceleration field ${\textstyle {\mathit {\boldsymbol {a}}}({\mathit {\boldsymbol {x}}},t)}$ can be derived from the velocity field as:

where ${\textstyle {d}_{t}}$ is the total derivative

Given an acceleration field, ${\textstyle {\mathit {\boldsymbol {a}}}}$, derived from a velocity field ${\textstyle {\mathit {\boldsymbol {u}}}}$ as in Eq. (4), the particles’ velocities and positions can be obtained integrating the equations of motion with initial conditions:

3. Numerical integration of particles’ equations of motion

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

In molecular dynamics, a well-known numerical integration scheme for the equations of motion is the velocity Verlet integrator [6,7], which reads:

where eq. (8) is explicit in time and provides the particle position. We can rewrite eqs. (8) and (9) as:

where ${\textstyle {\mathit {\boldsymbol {X}}}_{\lambda }^{n}={\mathit {\boldsymbol {X}}}_{\lambda }\left({t}^{n}\right)}$ , ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }^{n}=}$${\displaystyle {\mathit {\boldsymbol {U}}}_{\lambda }\left({t}^{n}\right)}$ , ${\textstyle {\mathit {\boldsymbol {a}}}^{n}\left({\mathit {\boldsymbol {X}}}_{\lambda }^{n}\right)=}$${\displaystyle {\mathit {\boldsymbol {a}}}\left({\mathit {\boldsymbol {X}}}_{\lambda }^{n},{t}^{n}\right)}$ .

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:

Where ${\textstyle P}$ is the fluid pressure, ${\textstyle \rho }$ is the fluid density, ${\textstyle \nu }$ is the kinematic viscosity and ${\textstyle {\mathit {\boldsymbol {f}}}}$ are mass forces. Eq. (12) can be expressed as ${\textstyle {d}_{t}{\mathit {\boldsymbol {u=}}}{\mathit {\boldsymbol {a}}}}$, where the acceleration field is given by:

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 ${\textstyle \left\{\lambda \right\}}$ be a set of fluid particles. Then the equation of motions can be integrated using the velocity Verlet integrator:

As stated above, this scheme is second order accurate over time. Let´s assume that all variables are known at time ${\textstyle {t}^{n}}$. Then the new position of the fluid particles at time ${\textstyle {t}^{n+1}}$ can be easily calculated using Eq.(15). Reordering Eq. (16):

where ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }^{n+1/2}={\mathit {\boldsymbol {U}}}_{\lambda }^{n}+}$${\displaystyle {\frac {\Delta t}{2}}{\mathit {\boldsymbol {a}}}^{n}({\mathit {\boldsymbol {X}}}_{\lambda }^{n})}$. In order to solve eq. (17), the term ${\textstyle {\mathit {\boldsymbol {a}}}^{n+1}}$ will be solved in an Eulerian framework. Then an equivalent Eulerian equation is imposed:

Where ${\textstyle {\mathit {\boldsymbol {u}}}^{n+1/2}}$ is obtained via projection:

${\textstyle {P}_{\left\{\lambda \right\}}^{n+1}}$ is a projector operator that transfer the values of the variables from the set of particle with positions ${\textstyle \left\{{\mathit {\boldsymbol {X}}}_{\lambda }^{n+1}\right\}}$ onto the mesh.

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

where ${\textstyle \psi }$ 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. (18) and (20):

which is equivalent to say that ${\textstyle \,{\mathit {\boldsymbol {u}}}^{n+1}\left({\mathit {\boldsymbol {x}}}\right)=}$${\displaystyle {P}_{\left\{\lambda \right\}}^{n+1}\left[\left\{{\mathit {\boldsymbol {U}}}_{\lambda }^{n+1}\right\}\right]}$ . There is another outcome from the fulfilment of the coherence condition. 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 ${\textstyle {\frac {\Delta t}{2}}{\mathit {\boldsymbol {a}}}^{n+1}({\mathit {\boldsymbol {X}}}_{\lambda }^{n+1})}$ , will be projected back to the mesh as it is, with no further need of correction.

Eq. (18) can be solved iteratively by further splitting into:

Where ${\textstyle i}$ is the iteration index. Now in order to obtain the pressure ${\textstyle {P}^{n+1,i+1}}$, the continuity condition is imposed. Introducing the continuity equation ( ${\textstyle \nabla \cdot {\mathit {\boldsymbol {u}}}^{n+1,i+1}=}$${\displaystyle 0}$) into Eq.(23) and reordering

Once the iterative process between Eqs. (22)-(23) has converged, from eq. (18) we obtain:

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. (22) and (24). First, we need to discretize it in space using a finite element mesh. Then, any intrinsic dependent variable ${\textstyle \psi }$ can be expressed as:

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

Once the spatial discretization is set with a FE mesh, the projector operator must transfer values from particles onto the background mesh following Eq. (21):

Where ${\textstyle {P}_{h,\left\{\lambda \right\}}}$ must fulfil the condition given in Eq. (20) leading to:

Using the FE Galerkin method as in [9,10] to solve the variational formulation of Eqs. (22) and (24) we get:

Where ${\textstyle {\overline {\overline {\mathit {\boldsymbol {M}}}}}}$ is the FE mass matrix, ${\textstyle {\overline {\overline {\mathit {\boldsymbol {L}}}}}}$ is the FE Laplacian matrix, ${\textstyle {\overline {\overline {\mathit {\boldsymbol {G}}}}}}$ is the FE gradient matrix, and ${\textstyle {\overline {\overline {\mathit {\boldsymbol {D}}}}}}$ is the FE divergence matrix. Once the iterative loop over ${\textstyle i}$ converges, then the new acceleration field is obtained from Eq. (18). The convergence criteria used for the pressure and velocity are:

5. Global least square projector

Let ${\textstyle \left\{{\Psi }_{\lambda }\right\}}$ be a set of values carried by a set of particles and ${\textstyle \left\{{\psi }_{b}^{\ast }\right\}}$ 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:

where ${\textstyle {\psi }_{h}\left({\mathit {\boldsymbol {X}}}_{\lambda }\right)=}$${\displaystyle \sum _{b}^{}{N}^{b}\left({\mathit {\boldsymbol {X}}}_{\lambda }\right){\psi }_{b}^{\ast }}$

Minimizing the least square error via ${\textstyle \partial {\epsilon }_{\psi }/\partial {\psi }_{b}^{\ast }=}$${\displaystyle 0}$:

Eq. (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 ${\textstyle {\Psi }_{\lambda }=}$${\displaystyle \sum _{c}^{}{N}^{c}\left({\mathit {\boldsymbol {X}}}_{\lambda }\right){\psi }_{c}}$ :

leading to ${\textstyle {\psi }_{c}^{\ast }={\psi }_{c}}$ and the fulfilment of the coherence condition given in eq. (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 ${\textstyle \psi }$ . Then the equalities ${\textstyle {\Psi }_{\lambda }=}$${\displaystyle \psi \left({\mathit {\boldsymbol {X}}}_{\lambda }\right)}$ and ${\textstyle \psi \left({\mathit {\boldsymbol {X}}}_{\lambda }\right)=}$${\displaystyle \sum _{c}^{}{N}^{c}\left({\mathit {\boldsymbol {X}}}_{\lambda }\right)\psi \left({\mathit {\boldsymbol {x}}}_{c}\right)+}$${\displaystyle O({h}_{e}^{2})}$ hold, where ${\textstyle {h}_{e}}$ is the characteristic length of the element where particle ${\textstyle \lambda }$ is at. Using Eq. (35):

And, the solution to eq. (37) can be expressed as:

Where ${\textstyle \delta {\psi }_{c}}$ is the solution of:

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

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 ${\textstyle {\mathit {\boldsymbol {u}}}\left({\mathit {\boldsymbol {x}}},{t}^{n}\right)=}$${\displaystyle {\mathit {\boldsymbol {u}}}^{n}\left({\mathit {\boldsymbol {x}}}\right)}$ and the corresponding acceleration field ${\textstyle {\mathit {\boldsymbol {a}}}\left({\mathit {\boldsymbol {x}}},{t}^{n}\right)=}$${\displaystyle {\mathit {\boldsymbol {a}}}^{n}\left({\mathit {\boldsymbol {x}}}\right)={d}_{t}{\mathit {\boldsymbol {u}}}^{n}\left({\mathit {\boldsymbol {x}}}\right)}$ , and a set of particles with initial condition ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }^{0}=}$${\displaystyle {\mathit {\boldsymbol {u}}}^{0}\left({\mathit {\boldsymbol {X}}}_{\lambda }^{0}\right)}$ , the velocity Verlet algorithm estimates the velocity and position at ${\textstyle {t}^{n}=}$${\displaystyle n\Delta t}$ with second order convergence:

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:

The velocity at the mesh nodes is initiated from the initial particle velocity ${\textstyle {\mathit {\boldsymbol {U}}}_{\lambda }^{0}}$ through the projection:

Where ${\textstyle {\boldsymbol {\epsilon }}^{P,0}}$ is the projection error. The projected velocity is then interpolated back to the particles using the projected values:

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

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

Considering that the projector and the interpolator are second order accurate, and using Eq. (8)-(9):

Which reads:

Eqs. (51)-(52) show that the error for the first time step is ${\textstyle O\left(\Delta {t}^{3}\right)+}$${\displaystyle O({h}_{e}^{2}\Delta t)}$. It is straightforward to show that the resulting errors for the following time steps are those related to the accumulation, which reads:

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