FERNÁDEZ et al

Quartier Polytech 1, B52/3, Allee de la Decouverte 13A, 4000 Liege, Belgium. *

# Generalized-$\alpha$ scheme in the PFEM for velocity-pressure and
displacement-pressure formulations of the incompressible Navier-Stokes
equations

Eduardo FERNÁNDEZ1,∗ Simon FÉVRIER1 Martin LACROIX1 Romain BOMAN1 Jean-
Philippe PONTHOT1 Aerospace and Mechanical Engineering Department, University
of Liège, Belgium<EMAIL_ADDRESS>
(26 April 2016; 6 June 2016; 6 June 2016)

###### Abstract

[Summary] Despite the increasing use of the Particle Finite Element Method
(PFEM) in fluid flow simulation and the outstanding success of the
Generalized-$\alpha$ time integration method, very little discussion has been
devoted to their combined performance. This work aims to contribute in this
regard by addressing three main aspects. Firstly, it includes a detailed
implementation analysis of the Generalized-$\alpha$ method in PFEM. The work
recognizes and compares different implementation approaches from the
literature, which differ mainly in the terms that are $\alpha$-interpolated
(state variables or forces of momentum equation) and the type of treatment for
the pressure in the time integration scheme. Secondly, the work compares the
performance of the Generalized-$\alpha$ method against the Backward Euler and
Newmark schemes for the solution of the incompressible Navier-Stokes
equations. Thirdly, the study is enriched by considering not only the
classical velocity-pressure formulation but also the displacement-pressure
formulation that is gaining interest in the fluid-structure interaction field.
The work is carried out using various 2D and 3D benchmark problems such as the
fluid sloshing, the solitary wave propagation, the flow around a cylinder, and
the collapse of a cylindrical water column.

###### keywords:

Time integration, Fluid mechanics, Particle methods, CFD, Generalized-Alpha,
PFEM

††articletype: PREPRINT

##  1 Introduction

The Particle Finite Element Method (PFEM) has drawn attention of the
simulation community due to the possibility to formulate fluid flow equations
in a Lagrangian framework, allowing the use of classical Lagrangian FEM and
easing tracking of fluid boundaries, even in case of large and unpredictable
boundary motions due to efficient remeshing algorithms 1. The method has been
extended to various materials and multi-physics problems with moving domains,
such as plasticity 2, fluid-structure interaction 3, 4 and phase change 5
among others 6, 7.

In the PFEM, Lagrangian-based governing equations are solved using spatial
FEM-Galerkin and temporal discretizations. The first one subdivides the
computational domain in a number of finite elements to approximate the
solution of partial differential equations governing the body motion. This
procedure, applied to the momentum equation of a Newtonian fluid, leads to the
following expression:

$\bm{f}^{\mathrm{dyn}}(\accentset{\mbox{\large.}}{\mathbf{v}},\mathbf{x},t)+\bm{f}^{\mathrm{int}}(\mathbf{v},\mathbf{p},\mathbf{x},t)=\bm{f}^{\mathrm{ext}}(\mathbf{x},t)$
(1)

where $\bm{f}^{\mathrm{dyn}}$, $\bm{f}^{\mathrm{int}}$ and
$\bm{f}^{\mathrm{ext}}$ represent, respectively, dynamic (inertial or
d’Alembert forces), internal and external forces acting on the nodes defining
the finite elements. These forces depend on the nodal acceleration
($\accentset{\mbox{\large.}}{\mathrm{v}}$), velocity ($\mathrm{v}$), position
($\mathrm{x}$), pressure (${\rm Phys.~{}Rev.~{}E}$), and time ($t$). The
momentum equation is complemented with the continuity equation to form the set
of Navier-Stokes equations. In addition, the transient problem requires a
temporal discretization to approximate time variation of state variables. For
example, starting from a known configuration of the body at time $t_{n}$, the
updated condition of the fluid at $t_{n+1}$ is governed by the discretized
momentum equation:

$\bm{f}^{\mathrm{dyn}}(\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1},\mathbf{x}_{n+1})+\bm{f}^{\mathrm{int}}(\mathbf{v}_{n+1},\mathbf{p}_{n+1},\mathbf{x}_{n+1})=\bm{f}^{\mathrm{ext}}(\mathbf{x}_{n+1})$
(2)

where the subscript $n$+$1$ indicates that variables are defined at time
$t_{n+1}$. To relate state variables with their time derivatives
($\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$, $\mathrm{v}_{n+1}$ and
$\mathrm{x}_{n+1}$), a time integration scheme must be used. For example,
common schemes in the PFEM literature are the implicit Backward Euler 3, 5, 8,
9, 10, 11, 12, Trapezoidal 13, 14, Newmark 15, 16 and Newmark-Bossak 17, 18,
19, 20. Although explicit schemes can also be found 21. However, one of the
best performing time integration schemes reported in the computational
dynamics literature, and which has not been assessed in the PFEM context yet,
is the Generalized-$\alpha$ method. This time integration scheme was proposed
by Chung and Hulbert 22 for solving dynamic equations in solid mechanics.
Later, Jansen et.al.23 extended the idea to fluid mechanics. Since then,
several works have demonstrated outstanding performance of the
Generalized-$\alpha$ method for solving the incompressible Navier-Stokes
equations, mostly using Eulerian finite element discretization 24. The
essential feature of the Generalized-$\alpha$ (GA) method is that it writes
state equations at times $t_{n+\alpha_{f}}$ and $t_{n+\alpha_{m}}$, instead of
doing so at time $t_{n+1}$. For this method, the momentum equation is as
follows:

$\bm{f}^{\mathrm{dyn}}_{n+\alpha_{m}}+\bm{f}^{\mathrm{int}}_{n+\alpha_{f}}-\bm{f}^{\mathrm{ext}}_{n+\alpha_{f}}=0$
(3)

For solving Eq. (3), two approaches can be identified in the literature, which
are denoted as GA-I and GA-II in this work. The first one (GA-I) is the most
reported one for fluid dynamics and follows the original contribution of Chung
and Hulbert 22 and Jansen et.al.23. GA-I computes nodal acceleration,
velocity, position and pressure at $t_{n+\alpha_{m}}$ or $t_{n+\alpha_{f}}$
assuming a linear combination between those at $t_{n}$ and $t_{n+1}$, as
follows:

$\displaystyle\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}$
$\displaystyle=(1-\alpha_{m})\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+\alpha_{m}\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$
(4a) $\displaystyle\mathrm{x}_{n+\alpha_{m}}$
$\displaystyle=(1-\alpha_{m})\>{\mathrm{x}}_{n}+\alpha_{m}\>{\mathrm{x}}_{n+1}$
(4b) $\displaystyle\mathrm{v}_{n+\alpha_{f}}$
$\displaystyle=(1-\alpha_{f})\>{\mathrm{v}}_{n}+\alpha_{f}\>{\mathrm{v}}_{n+1}$
(4c) $\displaystyle\mathrm{x}_{n+\alpha_{f}}$
$\displaystyle=(1-\alpha_{f})\>{\mathrm{x}}_{n}+\alpha_{f}\>{\mathrm{x}}_{n+1}$
(4d) $\displaystyle\mathrm{p}_{n+\alpha_{f}}$
$\displaystyle=(1-\alpha_{f})\>{\mathrm{p}}_{n}+\alpha_{f}\>{\mathrm{p}}_{n+1}$
(4e)

where $\alpha_{m}$ and $\alpha_{f}$ are user-defined parameters. In this way,
the momentum Eq. (3) is written as:

$\bm{f}^{\mathrm{dyn}}(\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha_{m}},\mathbf{x}_{n+\alpha_{m}})+\bm{f}^{\mathrm{int}}(\mathbf{v}_{n+\alpha_{f}},\mathbf{p}_{n+\alpha_{f}},\mathbf{x}_{n+\alpha_{f}})=\bm{f}^{\mathrm{ext}}(\mathbf{x}_{n+\alpha_{f}})$
(5)

Unavoidably, the GA-I approach (Eq. 5) results in a system of equations with
constitutive matrices that must be computed at times $t_{n+\alpha_{m}}$ and
$t_{n+\alpha_{f}}$. The second implementation approach (GA-II) follows the
rationale of Hilber-Hughes-Taylor (HHT-$\alpha$)25 and Wood-Bossak-Zienkiewicz
(WBZ-$\alpha$)26, who paved the way to the Generalized-$\alpha$ method. Their
idea is to formulate dynamic, internal and external forces of the momentum
equation at $t_{n+\alpha_{m}}$ and $t_{n+\alpha_{f}}$ using an interpolation
of forces between $t_{n}$ and $t_{n+1}$, as follows:

$\displaystyle\bm{f}^{\mathrm{dyn}}_{n+\alpha_{m}}$
$\displaystyle=(1-\alpha_{m})\>\bm{f}^{\mathrm{dyn}}_{n}+\alpha_{m}\>\bm{f}^{\mathrm{dyn}}_{n+1}$
(6a) $\displaystyle\bm{f}^{\mathrm{int}}_{n+\alpha_{f}}$
$\displaystyle=(1-\alpha_{f})\>\bm{f}^{\mathrm{int}}_{n}+\alpha_{f}\>\bm{f}^{\mathrm{int}}_{n+1}$
(6b) $\displaystyle\bm{f}^{\mathrm{ext}}_{n+\alpha_{f}}$
$\displaystyle=(1-\alpha_{f})\>\bm{f}^{\mathrm{ext}}_{n}+\alpha_{f}\>\bm{f}^{\mathrm{ext}}_{n+1}$
(6c)

When replacing Eq. (6) into the momentum balance Eq. (3), GA-II yields a
system of equations with state variables and constitutive matrices at $t_{n}$
and $t_{n+1}$. This offers advantages over the implementation approach GA-I.
For instance, it avoids computationally expensive algorithms or constitutive
matrices at $t_{n+\alpha_{f}}$ and $t_{n+\alpha_{m}}$, such as in problems
involving contact and/or plasticity 27, 28. Possibly for this reason, the GA-
II strategy can be easily found in solid mechanics, but hardly in Newtonian-
based fluid mechanics. From now on, various implementation approaches will be
introduced, which are diagrammed in Fig. 1 for the sake of clarity.

Figure 1: Classification of the different implementation approaches of the
Generalized-$\alpha$ method for solving the incompressible Navier-Stokes
equations, which results in six different implementations of the algorithm.

In addition, it can be identified that Eqs. (4a)-(4e) are used in two ways in
the literature. The first one consists of substituting Eq. (4) in the momentum
equation to obtain a system whose unknowns are the velocity and pressure at
$t_{n+1}$ (see, for example, Valdés Vázquez 29). This implementation approach
is denoted GA-I${}^{q_{n+1}}$ in this work. The second approach consists of
solving state equations directly at $t_{n+\alpha_{f}}$ to obtain the velocity
and pressure at $t_{n+\alpha_{f}}$. Then, using Eqs. (4a)-(4e) to retrieve
kinematic variables and pressure at $t_{n+1}$ (see for example Jansen
et.al.23). This implementation approach is denoted GA-I${}^{q_{n+\alpha}}$ in
this work, as illustrated in the second classification level of Fig. 1.

Another implementation detail that deserves attention is the way in which
pressure is treated in the time integration scheme. For example, two different
approaches can be identified in the literature for the incompressible Navier-
Stokes equations. The first one follows the principle of Jansen et.al.23, in
which the pressure is set at $t_{n+\alpha_{f}}$ in the momentum equation using
Eq. (4e). The second approach uses the pressure at $t_{n+1}$ (or sets
$\alpha_{f}=1$ in Eq. 4e while keeping $\alpha_{f}\neq 1$ in Eqs. 4c-4d).
These approaches are denoted as GA-I${}_{\color[rgb]{0,0,1}p_{n+\alpha}}$ and
GA-I${}_{\color[rgb]{1,0,0}{p_{n+1}}}$, respectively. Combining these with the
other above-mentioned implementation approaches, 6 different schemes of the
alpha-generalized method are eventually derived from the literature for a
velocity-pressure formulation:
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$, $\text{GA-
II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$, which are illustrated in the third
classification level of Fig. 1. For instance,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ is used in 30,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ in 23, 31, 32,
and $\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ in 24, 29, 33. The
other schemes, to the best of our knowledge, could not be recognized in the
fluid mechanics literature but they were derived accordingly to become part of
this study.

Despite the popularity of the Generalized-$\alpha$ method in solid and fluid
dynamics, the different approaches mentioned above have not been explored in
such detail as to recognize in advance which one suits better the PFEM. One of
the few works addressing such a variety of approaches is the recent work of
Liu et.al.34. The authors use stable B-spline (NURBS) elements for Eulerian
spatial discretization of Navier-Stokes equations to analyse the effect of
using the pressure at $t_{n+1}$ instead of $t_{n+\alpha}$
(GA-I${}_{\color[rgb]{1,0,0}p_{n+1}}$ versus
GA-I${}_{\color[rgb]{0,0,1}p_{n+\alpha}}$). The authors show that
GA-I${}_{\color[rgb]{1,0,0}p_{n+1}}$ reduces the performance of the
Generalized-$\alpha$ method since only first-order accuracy is achieved, at
least for the pressure. However, since the present work uses a spatial
Lagrangian PFEM discretization stabilized with the Pressure-Stabilizing
Petrov-Galerkin (PSPG) formulation, it is not possible to guarantee a priori
that using the pressure at $p_{n+\alpha_{f}}$ will lead to better performance
than using it at $p_{n+1}$. The argument is that, depending on the
implementation approach or on the formulation (velocity-based, or
displacement-based), it is possible to obtain a pressure term, $p_{n}$, acting
on the right-hand side of the momentum equation. This could result in spurious
oscillations of the external force vector due to numerical instabilities
inherent to the spatial discretization, and would be a reason to avoid the
$\alpha$-interpolation of the pressure.

One of the advantages provided by the PFEM is the possibility to identify
evolving fluid boundaries easily, which facilitates numerical schemes for
fluid-structure interaction, among others. In this context, advantages have
been observed in formulating the Navier-Stokes equations with respect to nodal
displacements and pressures instead of the classical velocity-pressure
approach, since it facilitates unified formulations for fluids and solids.
Some of the few authors introducing a displacement-pressure formulation in
PFEM for the incompressible Navier-Stokes equations are Ryzhakov et.al.20, who
use the Newmark-Bossak scheme when producing results but the Backward Euler
for description of their methodology. In the Lagrangian-based FEM for fluid
dynamics, the works of Radovitzky and Ortiz 35 and Avancini and Sanches36
present displacement-pressure and position-pressure formulations,
respectively. Notably, both works use the Newmark’s method as time integration
scheme. Avancini and Sanches36 present a total Lagrangian formulation without
remeshing while Radovitzky and Ortiz35 use an advancing front method as
remeshing technique, which constraints fluid topology. Thus, both works
experience difficulties in capturing wave breaking and splashing phenomena,
which is an essential feature of PFEM. To the best of our knowledge, no
authors in Lagrangian-based fluid mechanics, let alone in PFEM, include the
Generalized-$\alpha$ method in a displacement-based formulation.

This work aims to contribute to the PFEM performance for fluid flow
simulations by introducing the Generalized-$\alpha$ time integration scheme in
the formulation. In view of the foregoing discussion, this work is focused on
five aspects. Firstly, to incorporate the Generalized-$\alpha$ time
integration scheme into PFEM for solving the incompressible Navier-Stokes
equations in an updated Lagrangian framework. Secondly, to include all the
implementation approaches in order to identify differences or equivalences.
Thirdly, to verify whether or not omitting the intermediate pressure
($p_{n+\alpha}$) in the time integration scheme leads to a worse performance
of the Generalized-$\alpha$ method. Fourthly, to compare the performance of
the Generalized-$\alpha$ method with the popular choices in the PFEM
literature, such as the Backward Euler and Newmark schemes. Fifthly and
lastly, to carry out the analysis for both velocity-based and displacement-
based formulations, given the interest of the community in the latter for
fluid-structure interaction simulation. The study is carried out using
incompressible Newtonian fluids and well-reported benchmarks as the flow
around a cylinder, the fluid sloshing, the solitary wave propagation, and the
collapse of a 3D cylindrical water column.

Consistent with the CFD literature, results indicate that the
Generalized-$\alpha$ method in PFEM outperforms the Backward Euler and Newmark
schemes as it does not suffer from excessive numerical damping for large time
steps and exhibits less spurious oscillations than the classical Trapezoidal
rule. In turn, it is observed that implementation approaches assuming an
interpolation of state variables (GA-I) or of equilibrium forces (GA-II) lead
to similar results in the PFEM, at least for incompressible Newtonian fluids.
Regarding the treatment of the pressure, a reduction in performance of the
Generalized-$\alpha$ scheme was observed when the pressure is imposed at
$t_{n+1}$ instead of $t_{n+\alpha_{f}}$, which is in line with the
observations of Liu et.al34. Although two cases must be distinguished. The
first one leads to a large degradation of the time integration scheme and
results from integrating the pressure gradient in the current configuration
(at $t_{n+1}$). In contrast, the second case integrates the pressure gradient
in the intermediate configuration (at $t_{n+\alpha_{f}}$), which leads to a
minor degradation of the Generalized-$\alpha$ performance. Finally, findings
are likewise valid for the displacement-based formulation.

The reminder of this manuscript is organized as follows. Section 2 introduces
the discretized Navier-Stokes equations for an incompressible Newtonian fluid.
Section 3 presents the studied time integration schemes including all
implementation approaches of the Generalized-$\alpha$ method and describes the
PFEM implementation of this work. The benchmark problems, results and
discussions are provided in section 4. Finally, section 5 gathers the
conclusions of this work.

##  2 Space discretization

In a simplified fashion, the Particle Finite Element Method (PFEM) is composed
of two phases. Given a spatio-temporal FEM discretization of the fluid at time
tn, the first step is to solve the Lagrangian system of equations using FEM to
obtain the nodal velocity and pressure at tn+1. Then, nodal acceleration and
position are updated consistently with the time integration scheme. The second
stage in PFEM consists of remeshing the fluid domain in its updated
configuration to improve mesh quality when necessary. In the following, the
FEM part is presented for an incompressible Newtonian fluid. For a global
overview of PFEM on the simulation of other materials, the reader is referred
to7.

###  2.1 Governing equations

Incompressible Navier-Stokes equations can be obtained from the momentum and
mass conservation equations for a continuum. For a given time $t$ and current
configuration $\Omega(t)$, these equations read as follows:

$\displaystyle\rho\>\accentset{\mbox{\large.}}{\mathrm{v}}-\mu\Delta\mathrm{v}+\nabla{\rm
Phys.~{}Rev.~{}E}$ $\displaystyle=\mathrm{f}\;\;\;,\;\;\text{in \;}\Omega(t)$
(7a) $\displaystyle\nabla\cdot\mathrm{v}$ $\displaystyle=0\;\;\;,\;\;\text{in
\;}\Omega(t)$ (7b)

where $\rho$ is the density,
$\mathrm{v}=[\mathrm{v}_{x},\mathrm{v}_{y},\mathrm{v}_{z}]$ is the velocity
vector, $\Delta$ is the Laplacian operator, and $\mathrm{f}$ is the body force
vector. The point over the velocity
($\accentset{\mbox{\large.}}{\text{\small{$\square$}}}$) represents the
Lagrangian time derivative, thus $\accentset{\mbox{\large.}}{\mathrm{v}}$ is
the acceleration vector. The pressure and dynamic viscosity that define the
Cauchy stress tensor of a Newtonian fluid are denoted, respectively, as ${\rm
Phys.~{}Rev.~{}E}$ and $\mu$ in Eq. (7a).

The system of Eqs. (7) must be completed with Dirichlet and Neumann boundary
conditions to obtain a velocity-pressure formulation of the Navier-Stokes
equations. This formulation arises due to the constitutive law that defines
the Cauchy stress tensor as a function of pressure and velocity. Note that
such tensor can also be defined in terms of displacement and pressure, knowing
that:

$\mathrm{v}=\accentset{\mbox{\large.}}{\mathrm{d}}\;\;\;\;\text{and}\;\;\;\;\accentset{\mbox{\large.}}{\mathrm{v}}=\accentset{\mbox{\large..}}{\mathrm{d}}$
(8)

where $\mathrm{d}=[\mathrm{d}_{x},\mathrm{d}_{y},\mathrm{d}_{z}]$ is the
displacement vector defined as the difference between the current
($\mathrm{x}$) and reference ($\mathrm{x}_{0}$) position
($\mathrm{d}=\mathrm{x}-\mathrm{x}_{0}$). Using Eq. (8) in (7), the following
displacement-based formulation is obtained:

$\displaystyle\rho\>\accentset{\mbox{\large..}}{\mathrm{d}}-\mu\Delta\accentset{\mbox{\large.}}{\mathrm{d}}+\nabla{\rm
Phys.~{}Rev.~{}E}$ $\displaystyle=\mathrm{f}\;\;\;,\;\;\text{in \;}\Omega(t)$
(9a) $\displaystyle\nabla\cdot\accentset{\mbox{\large.}}{\mathrm{d}}$
$\displaystyle=0\;\;\;,\;\;\text{in \;}\Omega(t)$ (9b)

Since equations are written in a Lagrangian formalism, no convective effects
have to be taken into account. Thus, the same FEM-Galerkin procedure can be
applied to Eq. (7) or Eq. (9) to obtain the semi-discrete Navier-Stokes
equations

###  2.2 Finite element formulation

To facilitate the remeshing process required in PFEM, equal order linear
interpolation is used for both the pressure and kinematic variables. For this,
triangular (2D) and tetrahedral (3D) elements are used in the spatial
discretization of the governing equations. Inside finite elements, state
variables are interpolated as:

$\mathrm{d}=\mathbf{N}_{v}\;\bar{\mathrm{d}}\;\;\;,\;\;\;\mathrm{v}=\mathbf{N}_{v}\;\bar{\mathrm{v}}\;\;\;,\;\;\;{\rm
Phys.~{}Rev.~{}E}=\mathbf{N}_{p}\;\bar{{\rm Phys.~{}Rev.~{}E}}$ (10)

where $\mathbf{N}_{v}$ is a matrix and $\mathbf{N}_{p}$ is a vector containing
linear shape functions for kinematic variables and pressure, respectively. Bar
symbol in Eq. (10) denotes variables at the nodes of finite elements. For a
triangular (2D) element with nodes labeled as 1, 2 and 3, the aforementioned
entities are defined as:

$\displaystyle\begin{aligned}
\mathbf{N}_{v}=\left[\begin{matrix}N_{1}&N_{2}&N_{3}&0&0&0\\\\[4.30554pt]
0&0&0&N_{1}&N_{2}&N_{3}\end{matrix}\right]\;\;\;\;\;,\;\;\;\;\;\mathbf{N}_{p}=\left[\begin{matrix}N_{1}&N_{2}&N_{3}\end{matrix}\right]\end{aligned}$
(11a)
$\displaystyle\bar{\mathrm{d}}=\left[\begin{matrix}\mathrm{d}_{x,1}&\>\mathrm{d}_{x,2}&\>\mathrm{d}_{x,3}&\>\mathrm{d}_{y,1}&\>\mathrm{d}_{y,2}\>&\mathrm{d}_{y,3}\end{matrix}\right]^{\intercal}\;\;\;,\;\;\;\bar{\mathrm{v}}=\left[\begin{matrix}\mathrm{v}_{x,1}&\>\mathrm{v}_{x,2}&\>\mathrm{v}_{x,3}&\>\mathrm{v}_{y,1}&\>\mathrm{v}_{y,2}\>&\mathrm{v}_{y,3}\end{matrix}\right]^{\intercal}\;\;\;,\;\;\;\bar{{\rm
Phys.~{}Rev.~{}E}}=\left[\begin{matrix}{\rm Phys.~{}Rev.~{}E}_{1}&\>{\rm
Phys.~{}Rev.~{}E}_{2}&\>{\rm
Phys.~{}Rev.~{}E}_{3}\end{matrix}\right]^{\intercal}$ (11b)

The space-discretized momentum and continuity equations are written as
follows:

$\displaystyle\mathbf{M}\>\accentset{\mbox{\large.}}{\bar{\mathbf{v}}}+\mathbf{K}\>\bar{\mathbf{v}}-\mathbf{D}^{\intercal}\>\bar{\mathbf{p}}$
$\displaystyle=\bar{\mathbf{f}}$ (12a)
$\displaystyle\mathbf{C}\>\accentset{\mbox{\large.}}{\bar{\mathbf{v}}}+\mathbf{D}\>\bar{\mathbf{v}}+\mathbf{L}\>\bar{\mathbf{p}}$
$\displaystyle=\bar{\mathbf{h}}$ (12b)

and for a displacement-based formulation, the system becomes:

$\displaystyle\mathbf{M}\>\accentset{\mbox{\large..}}{\bar{\mathbf{d}}}+\mathbf{K}\>\accentset{\mbox{\large.}}{\bar{\mathbf{d}}}-\mathbf{D}^{\intercal}\>\bar{\mathbf{p}}$
$\displaystyle=\bar{\mathbf{f}}$ (13a)
$\displaystyle\mathbf{C}\>\accentset{\mbox{\large..}}{\bar{\mathbf{d}}}+\mathbf{D}\>\accentset{\mbox{\large.}}{\bar{\mathbf{d}}}+\mathbf{L}\>\bar{\mathbf{p}}$
$\displaystyle=\bar{\mathbf{h}}$ (13b)

In the momentum equation (Eqs. 12a and 13a), $\mathbf{M}$, $\mathbf{K}$ and
$\mathbf{D}$ are, respectively, the mass matrix, the matrix containing the
viscous terms, and the gradient matrix, all computed on the current
configuration $\Omega(t)$. $\bar{\mathbf{f}}$ is a vector containing the body
forces and surface tractions, denoted as $\bar{b}$ and $\bar{\mathrm{t}}$,
respectively. At elemental level, these entities are defined as:

$\mathbf{M}_{e}=\int_{\Omega_{e}}{\rho\,\mathbf{N}^{\intercal}_{v}\,\mathbf{N}_{v}}\,d\Omega\;\;\;,\;\;\;\mathbf{K}_{e}=\int_{\Omega_{e}}{\mathbf{B}^{\intercal}\,\mathbf{m}_{v}\,\mathbf{B}}\,d\Omega\;\;\;,\;\;\;\mathbf{D}_{e}=\int_{\Omega_{e}}{\mathbf{N}_{p}^{\intercal}\,\mathbf{m}_{p}\,\mathbf{B}}\,d\Omega\;\;\;,\;\;\;\bar{\mathbf{f}}_{e}=\int_{\Omega_{e}}{\rho\,\mathbf{N}_{v}^{\intercal}\,\bar{\bm{b}}}\,d\Omega+\int_{\Gamma_{t}}{\mathbf{N}_{v}^{\intercal}\bar{\mathbf{t}}}\,d\Gamma$
(14)

where $\mathbf{B}$ is a matrix containing the gradient of shape functions for
the velocity, $\mathbf{m}_{v}$ is a matrix with material viscous properties,
and $\mathbf{m}_{p}$ is a vector that sorts matrix B for the pressure
gradient. In 2D, these entities are defined as:

$\mathbf{B}=\left[\begin{matrix}\displaystyle\frac{\partial N_{1}}{\partial
x}\quad&\displaystyle\frac{\partial N_{2}}{\partial
x}\quad&\displaystyle\frac{\partial N_{3}}{\partial
x}\quad&0&0&0\\\\[8.61108pt] 0&0&0&\displaystyle\frac{\partial N_{1}}{\partial
y}&\displaystyle\frac{\partial N_{2}}{\partial y}&\displaystyle\frac{\partial
N_{3}}{\partial y}\\\\[8.61108pt] \displaystyle\frac{\partial N_{1}}{\partial
y}&\displaystyle\frac{\partial N_{2}}{\partial y}&\displaystyle\frac{\partial
N_{3}}{\partial y}&\displaystyle\frac{\partial N_{1}}{\partial
x}\quad&\displaystyle\frac{\partial N_{2}}{\partial
x}\quad&\displaystyle\frac{\partial N_{3}}{\partial x}\\\
\end{matrix}\right]\;\;\;,\;\;\;\;\;\;\mathbf{m}_{v}=\mu\,\left[\begin{matrix}\;\displaystyle\frac{4}{3}\quad&\quad-\displaystyle\frac{2}{3}\quad&\quad
0\\\\[8.61108pt]
-\displaystyle\frac{2}{3}\quad&\quad\;\displaystyle\frac{4}{3}\quad&\quad
0\\\\[8.61108pt] \;0\quad&\quad\;0\quad&\quad 1\\\
\end{matrix}\right]\;\;\;,\;\;\;\;\;\;\mathbf{m}_{p}=[1\quad 1\quad 0]$ (15)

The continuity equations (Eqs. 12b and 13b) include the Pressure-Stabilizing
Petrov-Galerkin (PSPG) stabilization 37. Its purpose is to circumvent the
saddle point problem, where pressure acts as a Lagrange multiplier of the
incompressibility constraint 38, and to satisfy the LBB condition
(Ladyzhenskaya-Babu$\check{\text{s}}$ka-Brezzi) since equal order linear
elements are used to discretize velocity (or displacement) and pressure. At
elemental level, the stabilization terms, $\mathbf{C}$, $\mathbf{L}$ and
$\bar{\mathbf{h}}$ are defined as:

$\mathbf{C}_{e}=\int_{\Omega_{e}}{\tau_{\mathrm{PSPG}}\,\nabla\mathbf{N}_{p}^{\intercal}\,\mathbf{N}_{v}}\,d\Omega\;\;\;,\;\;\;\mathbf{L}_{e}=\int_{\Omega_{e}}{\tau_{\mathrm{PSPG}}\,\frac{1}{\rho}\,\nabla\mathbf{N}_{p}^{\intercal}\,\nabla\mathbf{N}_{p}}\,d\Omega\;\;\;,\;\;\;\bar{\mathbf{h}}_{e}=\int_{\Omega_{e}}{\tau_{\mathrm{PSPG}}\,\rho\,\nabla\mathbf{N}_{p}^{\intercal}\,\bar{\bm{b}}}\,d\Omega$
(16)

with $\tau_{\mathrm{PSPG}}$ computed as39:

$\tau_{\mathrm{PSPG}}=\displaystyle\frac{1}{\left(\frac{4}{\Delta
t^{2}}+\frac{4\mathrm{v}^{\intercal}\mathrm{v}}{\tilde{h}^{2}}+\left(\frac{4\mu}{\rho\>\tilde{h}^{2}}\right)^{2}\right)^{\frac{1}{2}}}$
(17)

where $\Delta t$ is the time step and $\tilde{h}$ is a characteristic size
defined as the circumcircle diameter of the element.

For further implementation details concerning the terms in Eqs. (14) - (17),
the reader is referred to 40. To associate the semi-discretized equations with
the introductory discussion of this manuscript, it is pertinent to define the
forces involved in the momentum equation. The dynamic
($\bm{f}^{\mathrm{dyn}}$), internal ($\bm{f}^{\mathrm{int}}$) and external
($\bm{f}^{\mathrm{ext}}$) forces are defined as:

$\displaystyle\bm{f}^{\mathrm{dyn}}(t)$
$\displaystyle=\mathbf{M}\>\accentset{\mbox{\large.}}{\bar{\mathbf{v}}}=\mathbf{M}\>\accentset{\mbox{\large..}}{\bar{\mathbf{d}}}$
(18a) $\displaystyle\bm{f}^{\mathrm{int}}(t)$
$\displaystyle=\mathbf{K}\>\bar{\mathbf{v}}-\mathbf{D}^{\intercal}\>\bar{\mathbf{p}}=\mathbf{K}\>\accentset{\mbox{\large.}}{\bar{\mathbf{d}}}-\mathbf{D}^{\intercal}\>\bar{\mathbf{p}}$
(18b) $\displaystyle\bm{f}^{\mathrm{ext}}(t)$ $\displaystyle=\bar{\mathbf{f}}$
(18c)

Equations (12) and (13) are the final formulas of the momentum and continuity
equations that are solved in this work. However, to complete the system,
equations for the velocity and displacement time derivatives are needed, which
are given by the time integration methods presented below.

##  3 Time discretization and integration

To complete the system composed of momentum and continuity equations (Eqs. 12
and 13), time $t$ is discretized and time derivatives are integrated over the
time step [$t_{n}$ , $t_{n+1}$] as follows:

$\bar{\mathrm{v}}_{n+1}=\bar{\mathrm{v}}_{n}+\displaystyle\int_{t_{n}}^{t_{n+1}}\accentset{\mbox{\large.}}{\bar{\mathrm{v}}}(t)\>dt\;\;\;\;,\;\;\;\;\;\;\bar{\mathrm{d}}_{n+1}=\bar{\mathrm{x}}_{n+1}-\bar{\mathrm{x}}_{n}=\displaystyle\int_{t_{n}}^{t_{n+1}}\accentset{\mbox{\large.}}{\bar{\mathrm{d}}}(t)\>dt$
(19)

where $\bar{\mathrm{x}}_{n+1}$ and $\bar{\mathrm{x}}_{n}$ are the nodal
positions at $t_{n+1}$ and $t_{n}$, respectively. However, as
$\accentset{\mbox{\large.}}{\bar{\mathrm{v}}}(t)$ and
$\accentset{\mbox{\large.}}{\bar{\mathrm{d}}}(t)$ are unknown functions, they
must be approximated leading to different time integration schemes. The
approximate acceleration and velocity functions are denoted by
$\accentset{\mbox{\large.}}{v}$ and $v$, respectively. These hold for the time
step and are integrated as:

$\bar{\mathrm{v}}_{n+1}=\bar{\mathrm{v}}_{n}+\displaystyle\int_{0}^{\Delta
t}\accentset{\mbox{\large.}}{v}(\tau)\>d\tau\;\;\;\;,\;\;\;\;\;\;\bar{\mathrm{d}}_{n+1}=\displaystyle\int_{0}^{\Delta
t}v(\tau)\>d\tau$ (20)

where $\Delta t=t_{n+1}-t_{n}$ is the time step and $\tau$ is time during a
time step. In the following, two set of approximation formulas are presented,
the Backward Euler and Newmark. For now on, equations deal with nodal
variables only, so the bar symbol will be skipped for the sake of simplicity.

###  3.1 Backward Euler

The Backward Euler time integration scheme defines the unknown acceleration
function $\accentset{\mbox{\large.}}{v}(\tau)$ as constant and equal to the
acceleration of the current time step
$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$, leading to :

$\begin{matrix}\mathrm{v}_{n+1}&=\mathrm{v}_{n}+\displaystyle\int_{0}^{\Delta
t}\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\>d\tau=\mathrm{v}_{n}+\Delta
t\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\end{matrix}$ (21)

Similarly, Backward Euler can be applied for integrating the velocity and
obtain the displacement as $\mathrm{d}_{n+1}=\Delta t\>\mathrm{v}_{n+1}$.
Nevertheless, for consistency with Newmark’s method that is presented later,
displacement is obtained from integrating the acceleration twice, as follows:

$\displaystyle\mathrm{v}(\tau)$
$\displaystyle=\mathrm{v}_{n}+\displaystyle\int_{0}^{\tau}\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}dt=\mathrm{v}_{n}+\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$
(22) $\displaystyle\mathrm{d}_{n+1}$
$\displaystyle=\displaystyle\int_{0}^{\Delta
t}\mathrm{v}(\tau)\>d\tau=\displaystyle\int_{0}^{\Delta
t}\left[\mathrm{v}_{n}+\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\right]\>d\tau=\mathrm{v}_{n}\>\Delta
t+\frac{\Delta t^{2}}{2}\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ (23)

Kinematic relationships provided by the time integration scheme must be used
in the discretized momentum and continuity equations, either in the velocity-
pressure system (Eq. 12) or displacement-pressure system (Eq. 13). In this
work, both cases are considered. However, the algebraic substitution is not
detailed to avoid overextending the manuscript. The resulting momentum
equations are summarized in the first row of tables provided in appendices A
and B.

###  3.2 Newmark’s Method

This method41 can be thought of as one that approximates the unknown
acceleration function $\accentset{\mbox{\large.}}{v}(\tau)$ by a constant
acceleration equal to a linear combination between the current and previous
acceleration, as follows:

$\accentset{\mbox{\large.}}{v}(\tau)=(1-\gamma)\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+\gamma\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$
(24)

where $\gamma$ is an interpolation parameter. Using Eq. (24) in (20), the
well-known Newmark equation for nodal velocities is obtained:

$\mathrm{v}_{n+1}=\mathrm{v}_{n}+\displaystyle\int_{0}^{\Delta
t}\left[(1-\gamma)\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+\gamma\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\right]\>d\tau=\mathrm{v}_{n}+(1-\gamma)\>\Delta
t\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+\gamma\>\Delta
t\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ (25)

To compute nodal displacements, the Newmark’s method similarly approximates
the acceleration function $\accentset{\mbox{\large.}}{v}(\tau)$ by an
interpolation between the current and previous acceleration, but now
controlled by a parameter $2\beta$, as follows:

$\accentset{\mbox{\large.}}{v}(\tau)=(1-2\beta)\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+2\beta\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$
(26)

This function (Eq. 26) is integrated twice to obtain the well-known Newmark
equation for nodal displacements:

$\displaystyle\mathrm{v}(\tau)$
$\displaystyle=\mathrm{v}_{n}+\displaystyle\int_{0}^{\tau}\left[(1-2\beta)\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+2\beta\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\right]\>dt=\mathrm{v}_{n}+(1-2\beta)\>\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+2\beta\>\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$
(27) $\displaystyle\mathrm{d}_{n+1}$
$\displaystyle=\displaystyle\int_{0}^{\Delta
t}\left[\mathrm{v}_{n}+(1-2\beta)\>\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+2\beta\>\tau\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}\right]\>d\tau=\Delta
t\>\mathrm{v}_{n}+\frac{1-2\beta}{2}\Delta
t^{2}\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n}+\beta\Delta
t^{2}\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ (28)

The expressions given by Newmark’s method allow to complete the system
composed by the momentum and continuity equation, either in the velocity-based
(Eq. 12) or in the displacement-based (Eq. 13) formulation. To avoid
overextending the manuscript, the procedure is shown below only for the
displacement-based formulation. The reader is referred to Appendix A (second
row of table) to obtain the momentum equation for the velocity-based
formulation.

Using Eq. (25) to isolate the acceleration:

$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}=\frac{\mathrm{v}_{n+1}-\mathrm{v}_{n}}{\gamma\Delta
t}-\frac{1-\gamma}{\gamma}\accentset{\mbox{\large.}}{\mathrm{v}}_{n}$ (29)

and doing so also in Eq. (28):

$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}=\frac{\mathrm{d}_{n+1}}{\beta\Delta
t^{2}}-\frac{\mathrm{v}_{n}}{\beta\Delta
t}-\frac{1-2\beta}{2\beta}\accentset{\mbox{\large.}}{\mathrm{v}}_{n}$ (30)

velocity in terms of displacement can be obtained using Eqs. (29) and (30):

$\mathrm{v}_{n+1}=\frac{\gamma}{\beta\Delta
t}\mathrm{d}_{n+1}+\frac{\beta-\gamma}{\beta}\mathrm{v}_{n}+\frac{2\beta-\gamma}{2\beta}\Delta
t\accentset{\mbox{\large.}}{\mathrm{v}}_{n}$ (31)

Having expressed nodal acceleration and velocity in terms of displacement,
these are replaced in Eq. (13) to give:

$\displaystyle\mathbf{M}_{n+1}\left(\frac{\mathbf{d}_{n+1}}{\beta\Delta
t^{2}}-\frac{\mathbf{v}_{n}}{\beta\Delta
t}-\frac{1-2\beta}{2\beta}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{K}_{n+1}\left(\frac{\gamma}{\beta\Delta
t}\mathbf{d}_{n+1}+\frac{\beta-\gamma}{\beta}\mathbf{v}_{n}+\frac{2\beta-\gamma}{2\beta}\Delta
t\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)-\mathbf{D}^{\intercal}_{n+1}\>{\mathbf{p}}_{n+1}$
$\displaystyle={\mathbf{f}}_{n+1}\;\;,$ (32a)
$\displaystyle\mathbf{C}_{n+1}\left(\frac{\mathbf{d}_{n+1}}{\beta\Delta
t^{2}}-\frac{\mathbf{v}_{n}}{\beta\Delta
t}-\frac{1-2\beta}{2\beta}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{D}_{n+1}\left(\frac{\gamma}{\beta\Delta
t}\mathbf{d}_{n+1}+\frac{\beta-\gamma}{\beta}\mathbf{v}_{n}+\frac{2\beta-\gamma}{2\beta}\Delta
t\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{L}_{n+1}\>{\mathbf{p}}_{n+1}$
$\displaystyle={\mathbf{h}}_{n+1}\;\;,$ (32b)

Arranging terms in a matrix form leads to:

$\left[\begin{matrix}\mathbf{M}_{n+1}\>\displaystyle\frac{1}{\beta\Delta
t^{2}}+\mathbf{K}_{n+1}\>\displaystyle\frac{\gamma}{\beta\Delta t}&\hskip
14.22636pt-\mathbf{D}^{\intercal}_{n+1}\\\\[12.91663pt]
\mathbf{C}_{n+1}\>\displaystyle\frac{1}{\beta\Delta
t^{2}}+\mathbf{D}_{n+1}\>\displaystyle\frac{\gamma}{\beta\Delta t}&\hskip
14.22636pt\mathbf{L}_{n+1}\end{matrix}\right]\;\left[\begin{matrix}\mathbf{d}_{n+1}\\\\[17.22217pt]
\mathbf{p}_{n+1}\end{matrix}\right]=\left[\begin{matrix}\mathbf{f}_{n+1}+\mathbf{M}_{n+1}\left(\displaystyle\frac{\mathbf{v}_{n}}{\beta\Delta
t}+\frac{1-2\beta}{2\beta}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{K}_{n+1}\left(\displaystyle\frac{\gamma-\beta}{\gamma}\mathbf{v}_{n}+\displaystyle\frac{\gamma-2\beta}{2\beta}\Delta
t\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)\\\\[12.91663pt]
\mathbf{h}_{n+1}+\mathbf{C}_{n+1}\left(\displaystyle\frac{\mathbf{v}_{n}}{\beta\Delta
t}+\frac{1-2\beta}{2\beta}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{D}_{n+1}\left(\displaystyle\frac{\gamma-\beta}{\gamma}\mathbf{v}_{n}+\displaystyle\frac{\gamma-2\beta}{2\beta}\Delta
t\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)\end{matrix}\right]$ (33)

Eq. (33) is the final expresion of the discretized Navier-Stokes equations
solved in this work using the Newmark’s Method in a displacement-pressure
formulation. The momentum equation is also given in Appendix B (second row of
table) to ease comparison with equations obtained with the other time
integration schemes.

The $\gamma$ and $\beta$ parameters of Newmark’s method allow to control
accuracy and numerical damping, two competing criteria that require the user
to find a compromise. Notably, $\gamma=1$ and $\beta=0.5$ lead to the highly-
damped Backward Euler method defined previously, and $\gamma=0.5$ and
$\beta=0.25$ results in the classical undamped (but second order accurate)
Trapezoidal rule. When comparing tables from Appendices A and B, it can be
observed that the velocity-pressure formulation does not contain the $\beta$
parameter, since this is assigned to the displacement term. For this reason,
velocity-pressure formulations in Eulerian-based frameworks do not invoke such
a parameter, as can be appreciated in the manuscript of Jansen et.al.23. Thus,
to interpolate the Backward Euler and Trapezoidal rule through the Newmark’s
method, we resort to the linear stability analysis developed for solid
mechanics 42, 43. This establishes that, for a linear system, unconditional
stability is ensured if:

$\gamma\geq\frac{1}{2}\;\;\;\;\text{and}\;\;\;\;\beta\geq\frac{1}{4}\left(\gamma+\frac{1}{2}\right)^{2}$
(34)

For the numerical examples of this work, $\gamma$ is chosen and $\beta$ is
obtained as the smallest value satisfying Eq. (34). This applies to both the
Newmark and the Generalized-$\alpha$ methods.

As discussed in section 1, Backward Euler dominates the PFEM literature
presumably for its ease of implementation, and few works can be found using
the Newmark’s method44. However, one of the most successful time integration
schemes in solid and fluid mechanics, which has not been assessed in the PFEM
context to the best of our knowledge, is the Generalized-$\alpha$ presented
below.

###  3.3 Generalized-$\alpha$ method

Unlike the Backward Euler and Newmark formulas that complete the system
composed of the momentum and continuity equations, the Generalized-$\alpha$
method seeks to increase accuracy of the time integration scheme without
compromising numerical damping for high frequencies. This method is based on
the Newmark’s method and was proposed by Chung and Hulbert22 in the context of
solid mechanics and then extended to fluid mechanics by Jansen et.al23. The
Generalized-$\alpha$ method consists of posing the system of equations for a
time between $t_{n}$ and $t_{n+1}$. The time of the dynamic, internal and
external forces are defined differently, as this provides control over
accuracy and numerical dissipation of high frequency modes. Commonly, dynamic
forces are writen at $t_{n+\alpha_{m}}$, while internal and external forces at
$t_{n+\alpha_{f}}$. These intermediate times are expressed in terms of
weighting parameters $\alpha_{m}$ and $\alpha_{f}$, as follows:

$\begin{matrix}t_{n+\alpha_{m}}=(1-\alpha_{m})\>t_{n}+\alpha_{m}\>t_{n+1}=t_{n}+\alpha_{m}\>\Delta
t\\\\[8.61108pt]
t_{n+\alpha_{f}}=(1-\alpha_{f})\>t_{n}+\alpha_{f}\>t_{n+1}=t_{n}+\alpha_{f}\>\Delta
t\end{matrix}$ (35)

As discussed in the introduction of this manuscript, there are two
implementation approaches of Eq. (3) in the literature, which are denoted as
GA-I and GA-II. The first assumes that state variables scale linearly between
$t_{n}$ and $t_{n+1}$, which leads to the Eqs. (4a)-(4e). The second assumes a
linear variation for the equilibrium forces, leading to Eqs. (6a)-(6c). It can
be easily demonstrated that under specific conditions, both approaches are
equivalent. For example, if mass conservation is assumed:
$\mathbf{M}_{n}=\mathbf{M}_{n+\alpha_{m}}=\mathbf{M}_{n+1}$, then,

$\displaystyle\bm{f}^{\mathrm{dyn}}_{n+\alpha_{m}}$
$\displaystyle=\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha}$
(36a)
$\displaystyle=\mathbf{M}_{n+\alpha_{m}}((1-\alpha_{m})\accentset{\mbox{\large.}}{\mathbf{v}}_{n}+\alpha_{m}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1})\hskip
65.44142pt\text{(GA-I assumption)}$ (36b)
$\displaystyle=(1-\alpha_{m})\>\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}+\alpha_{m}\>\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1}$
(36c)
$\displaystyle=(1-\alpha_{m})\>\mathbf{M}_{n}\mathbf{v}_{n}+\alpha_{m}\>\mathbf{M}_{n+1}\mathbf{v}_{n+1}\hskip
59.75095pt\text{(because of mass conservation)}$ (36d)
$\displaystyle=(1-\alpha_{m})\>\bm{f}^{\mathrm{dyn}}_{n}+\alpha_{m}\>\bm{f}^{\mathrm{dyn}}_{n+1}\hskip
85.35826pt\text{(GA-II assumption)}$ (36e)

This shows that scaling the acceleration (GA-I) or the dynamic force (GA-II)
at $t_{n+\alpha_{m}}$ is equivalent under mass conservation. Regarding the
internal force, similar development can be performed to show that there is a
similarity between GA-I and GA-II. Taking the viscous forces for analysis:

$\displaystyle\overbrace{\mathbf{K}_{n+\alpha_{f}}\>((1-\alpha_{f})\>\mathbf{v}_{n}+\alpha_{f}\mathbf{v}_{n+1})}^{\text{GA-I}}$
$\displaystyle\approx$
$\displaystyle\overbrace{(1-\alpha_{f})\>\mathbf{K}_{n}\>\mathbf{v}_{n}+\alpha_{f}\>\mathbf{K}_{n+1}\>\mathbf{v}_{n+1}}^{\text{GA-
II}}$ (37a)
$\displaystyle\text{\lx@texthl@color{${(1-\alpha_{f})}$}}\>\mathbf{K}_{n+\alpha_{f}}\>\text{\lx@texthl@color{${\mathbf{v}_{n}}$}}+\text{\lx@texthl@color{${\alpha_{f}}$}}\mathbf{K}_{n+\alpha_{f}}\>\text{\lx@texthl@color{${\mathbf{v}_{n+1}}$}}$
$\displaystyle\approx$
$\displaystyle\text{\lx@texthl@color{${(1-\alpha_{f})}$}}\>\mathbf{K}_{n}\>\text{\lx@texthl@color{${\mathbf{v}_{n}}$}}+\text{\lx@texthl@color{${\alpha_{f}}$}}\>\mathbf{K}_{n+1}\>\text{\lx@texthl@color{${\mathbf{v}_{n+1}}$}}$
(37b)

where highlighted terms are those appearing on both sides of Eq. (37b). Note
that GA-I computes matrix $\mathbf{K}$ in $\Omega(t_{n+\alpha_{f}})$, while in
GA-II does it in $\Omega(t_{n})$ and $\Omega(t_{n+1})$. These tend to the same
matrix in case of small deformation where
$\Omega(t_{n})\approx\Omega(t_{n+\alpha_{f}})\approx\Omega(t_{n+1})$, or when
time step tends to zero, since if $t_{n+1}\to t_{n}\Rightarrow
t_{n+\alpha_{f}}\to t_{n}$. Thus, the Generalized-$\alpha$ time integration
method should produce similar results under GA-I or GA-II schemes in problems
featuring small deformations or using small time steps. Also, for linear
problems, scaling state variables (GA-I) or forces in the momentum equation
(GA-II) is equivalent. As Newtonian fluids assume a linear relationship
between shear stress and shear rate, GA-I and GA-II are expected to provide
similar results in this work. Equations using GA-I for solving the
incompressible Navier-Stokes equations (Eqs. 12 and 13) are presented below.

####  3.3.1 GA-I

To implement GA-I, it suffices to replace into Eq. (3) the equations developed
for the incompressible Newtonian fluid considered in this work. In a velocity-
pressure formulation, Eq. (3) reads as follows:

$\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha_{m}}+\mathbf{K}_{n+\alpha_{f}}\mathbf{v}_{n+\alpha_{f}}-\mathbf{D}^{\intercal}_{n+\alpha_{f}}\>\mathbf{p}_{n+\alpha_{f}}=\mathbf{f}_{n+\alpha_{f}}$
(38)

State variables $\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}$,
$\mathrm{v}_{n+\alpha_{f}}$ and ${\rm Phys.~{}Rev.~{}E}_{n+\alpha_{f}}$ are
defined by the Eqs. (4a), (4c) and (4e), which depend on the state variables
at time $t_{n}$ and $t_{n+1}$. The Generalized-$\alpha$ method uses the
Newmark formulas (Eqs. 25 and 28) to relate state variables at $t_{n+1}$ and
their time derivatives. As mentioned in the introduction, this set of
equations (Eqs. (4a)-(4e), (25) and (28)) are used in two ways in the
literature, which were denoted by GA-I${}^{q_{n+\alpha}}$ and
GA-I${}^{q_{n+1}}$. The system of equations obtained with these two approaches
are presented below.

##### GA-I${}^{q_{n+\alpha}}$

This approach writes the unknowns of the momentum equation at times
$t_{n+\alpha_{f}}$. In a velocity-pressure formulation, Eq. (38) is expressed
in terms of $\mathrm{v}_{n+\alpha_{f}}$ and $p_{n+\alpha_{f}}$. To do so, Eqs.
(25), (4a) and (4c) can be used to get the relationship between
$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}$ and
$\mathrm{v}_{n+\alpha_{f}}$ while eliminating
$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ and $\mathrm{v}_{n+1}$, which
is:

$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}=\frac{\alpha_{m}}{\gamma\>\alpha_{f}\>\Delta
t}\left(\mathrm{v}_{n+\alpha_{f}}-\mathrm{v}_{n}\right)+\left(1-\frac{\alpha_{m}}{\gamma}\right)\accentset{\mbox{\large.}}{\mathrm{v}}_{n}$
(39)

Replacing (39) into (38) and arranging terms, the following momentum equation
is obtained:

$\left(\mathbf{M}_{n+\alpha_{m}}\>\frac{\alpha_{m}}{\gamma\>\alpha_{f}\>\Delta
t}+\mathbf{K}_{n+\alpha_{f}}\right)\mathbf{v}_{n+\alpha_{f}}-\mathbf{D}^{\intercal}_{n+\alpha_{f}}\mathbf{p}_{n+\alpha_{f}}=\mathbf{f}_{n+\alpha_{f}}+\mathbf{M}_{n+\alpha_{m}}\left(\frac{\alpha_{m}}{\gamma\>\alpha_{f}\>\Delta
t}\>\mathrm{v}_{n}-\frac{\gamma-\alpha_{m}}{\gamma}\>\dot{\mathrm{v}}_{n}\right)$
(40)

Applying PSGP stabilization, the system of equations to be solved becomes:

$\begin{matrix}\left[\begin{matrix}\mathbf{M}_{{n+\alpha_{m}}}\>\displaystyle\frac{{\alpha_{m}}}{{\alpha_{f}}\gamma\Delta
t}+\mathbf{K}_{{n+\alpha_{f}}}&\hskip
11.38109pt-\mathbf{D}^{\intercal}_{{n+\alpha_{f}}}\\\\[17.22217pt]
\mathbf{C}_{{n+\alpha_{m}}}\>\displaystyle\frac{{\alpha_{m}}}{{\alpha_{f}}\gamma\Delta
t}+\mathbf{D}_{{n+\alpha_{f}}}&\hskip
11.38109pt\mathbf{L}_{{n+\alpha_{f}}}\end{matrix}\right]\;\left[\begin{matrix}\mathbf{v}_{{n+\alpha_{f}}}\\\\[21.52771pt]
\mathbf{p}_{{n+\alpha_{f}}}\end{matrix}\right]=\left[\begin{matrix}\mathbf{f}_{{n+\alpha_{f}}}+\mathbf{M}_{{n+\alpha_{m}}}\left(\displaystyle\frac{{\alpha_{m}}}{\gamma\>{\alpha_{f}}\>\Delta
t}\mathbf{v}_{\mathrm{n}}+\displaystyle\frac{{\alpha_{m}}-\gamma}{\gamma}\dot{\mathbf{v}}_{\mathrm{n}}\right)\\\\[12.91663pt]
\mathbf{h}_{{n+\alpha_{f}}}+\mathbf{C}_{{n+\alpha_{m}}}\left(\displaystyle\frac{{\alpha_{m}}}{\gamma\>{\alpha_{f}}\>\Delta
t}\mathbf{v}_{\mathrm{n}}+\displaystyle\frac{{\alpha_{m}}-\gamma}{\gamma}\dot{\mathbf{v}}_{\mathrm{n}}\right)\end{matrix}\right]\end{matrix}$
(41)

After solving Eq. (41), state variables at $t_{n+1}$ are computed with
equations given by the GA-I assumption (Eq. 4) and Newmark’s formulas.
Expressions for updating state variables once the system of equations is
solved are listed in Appendix C (third row of table).

For a displacement-based formulation, acceleration
$\accentset{\mbox{\large..}}{\mathrm{d}}_{n+\alpha_{m}}$ and velocity
$\accentset{\mbox{\large.}}{\mathrm{d}}_{n+\alpha_{f}}$ must be expressed in
terms of displacements ${\mathrm{d}}_{n+\alpha_{f}}$. However, Newmark’s
formulas relating state variables and their time derivatives are valid for
$t_{n+1}$, and not for $t_{n+\alpha}$. As a consequence, algebraic
substitution from the available formulas leads to a system whose unknowns are
displacement at $t_{n+1}$ and pressure at $t_{n+\alpha_{f}}$, as shown in
Appendix B (third row of table). To arrive at such a system of equations, the
acceleration $\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}$ (equal to
$\accentset{\mbox{\large..}}{\mathrm{d}}_{n+\alpha_{m}}$) in terms of
displacements $\mathrm{d}_{n+1}$ is obtained using Eqs. (30) and (4a), which
leads to:

$\accentset{\mbox{\large.}}{\mathrm{v}}_{n+\alpha_{m}}=\left(1-\frac{\alpha_{m}}{2\beta}\right)\accentset{\mbox{\large.}}{\mathrm{v}}_{n}-\frac{\alpha_{m}}{\beta\Delta
t}\mathrm{v}_{n}+\frac{\alpha_{m}}{\beta\Delta t^{2}}\mathrm{d}_{n+1}$ (42)

Then, velocities $\mathrm{v}_{n+\alpha_{f}}$ (equal to
$\accentset{\mbox{\large.}}{\mathrm{d}}_{n+\alpha_{f}}$) in terms of
displacements $\mathrm{d}_{n+1}$ are obtained using Eqs. (31) and (4c), which
leads to:

$\mathrm{v}_{n+\alpha_{f}}=\left(1-\frac{\alpha_{f}\gamma}{\beta}\right)\mathrm{v}_{n}+\alpha_{f}\Delta
t\left(1-\frac{\gamma}{2\beta}\right)\dot{\mathrm{v}}_{n}+\frac{\alpha_{f}\gamma}{\beta\Delta
t}\>\mathrm{d}_{n+1}$ (43)

Replacing Eqs. (42) and Eqs. (43) in the momentum equation (38) and adding the
PSPG stabilization, the system of equations becomes:

$\begin{matrix}\left[\begin{matrix}\mathbf{M}_{{n+\alpha_{m}}}\>\frac{{\alpha_{m}}}{\beta\Delta
t^{2}}+\mathbf{K}_{{n+\alpha_{f}}}\frac{{\alpha_{f}}\gamma}{\beta\Delta
t}&\>\>\>\>\>\>-\mathbf{D}^{\intercal}_{{n+\alpha_{f}}}\\\\[17.22217pt]
\mathbf{C}_{{n+\alpha_{m}}}\>\frac{{\alpha_{m}}}{\beta\Delta
t^{2}}+\mathbf{D}_{{n+\alpha_{f}}}\frac{{\alpha_{f}}\gamma}{\beta\Delta
t}&\>\>\>\>\>\>\mathbf{L}_{{n+\alpha_{f}}}\end{matrix}\right]\;\left[\begin{matrix}\mathbf{d}_{n+1}\\\\[21.52771pt]
\mathbf{p}_{{n+\alpha_{f}}}\end{matrix}\right]=\left[\begin{matrix}\mathbf{f}_{{n+\alpha_{f}}}+\mathbf{M}_{{n+\alpha_{m}}}\left(\frac{{\alpha_{m}}}{\beta\Delta
t}\mathbf{v}_{n}+\frac{{\alpha_{m}}-2\beta}{2\beta}\dot{\mathbf{v}}_{n}\right)-\mathbf{K}_{{n+\alpha_{f}}}\left(\frac{\beta-{\alpha_{f}}\gamma}{\beta}\mathbf{v}_{n}+\text{\footnotesize{${\alpha_{f}}\Delta
t$}}\frac{2\beta-\gamma}{2\beta}\dot{\mathbf{v}}_{n}\right)\\\\[12.91663pt]
\mathbf{h}_{{n+\alpha_{f}}}+\mathbf{C}_{{n+\alpha_{m}}}\left(\frac{{\alpha_{m}}}{\beta\Delta
t}\mathbf{v}_{n}+\frac{{\alpha_{m}}-2\beta}{2\beta}\dot{\mathbf{v}}_{n}\right)-\mathbf{D}_{{n+\alpha_{f}}}\left(\frac{\beta-{\alpha_{f}}\gamma}{\beta}\mathbf{v}_{n}+\text{\footnotesize{${\alpha_{f}}\Delta
t$}}\frac{2\beta-\gamma}{2\beta}\dot{\mathbf{v}}_{n}\right)\end{matrix}\right]\end{matrix}$
(44)

##### GA-I${}^{q_{n+1}}$

Instead of solving the system at $t_{n+\alpha_{f}}$, the GA-I parameterization
can be used to replace terms in the momentum equation to avoid states at
$t_{n+\alpha_{f}}$, as follows:

$\mathbf{M}_{n+\alpha_{m}}\overbrace{((1-\alpha_{m})\>\accentset{\mbox{\large.}}{\mathbf{v}}_{n}+\alpha_{m}\>\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1})}^{\text{\normalsize{$\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha_{m}}$}}}\>+\>\mathbf{K}_{n+\alpha_{f}}\overbrace{((1-\alpha_{f})\>\mathbf{v}_{n}+\alpha_{f}\>\mathbf{v}_{n+1})}^{\text{\normalsize{$\mathbf{v}_{n+\alpha_{f}}$}}}\>-\>\mathbf{D}^{\intercal}_{n+\alpha_{f}}\>\overbrace{((1-\alpha_{f})\>\mathbf{p}_{n}+\alpha_{f}\>\mathbf{p}_{n+1})}^{\text{\normalsize{$\mathbf{p}_{n+\alpha_{f}}$}}}\>=\>\mathbf{f}_{n+\alpha_{f}}$
(45)

Arranging states at $t_{n}$ in the right-hand side:

$\mathbf{M}_{n+\alpha_{m}}\>\alpha_{m}\>\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1}+\mathbf{K}_{n+\alpha_{f}}\>\alpha_{f}\>\mathbf{v}_{n+1}-\mathbf{D}^{\intercal}_{n+\alpha_{f}}\>\alpha_{f}\>\mathbf{p}_{n+1}=\mathbf{f}_{n+\alpha_{f}}-\mathbf{M}_{n+\alpha_{m}}(1-\alpha_{m})\>\accentset{\mbox{\large.}}{\mathbf{v}}_{n}-\mathbf{K}_{n+\alpha_{f}}(1-\alpha_{f})\>\mathbf{v}_{n}+\mathbf{D}^{\intercal}_{n+\alpha_{f}}\>(1-\alpha_{f})\>\mathbf{p}_{n}$
(46)

Resorting to Newmark’s formulas:

$\begin{matrix}\hskip
14.22636pt\mathbf{M}_{n+\alpha_{m}}\alpha_{m}\left(\displaystyle\frac{\mathbf{v}_{n+1}-\mathbf{v}_{n}}{\gamma\Delta
t}-\displaystyle\frac{1-\gamma}{\gamma}\>\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\right)+\mathbf{K}_{n+\alpha_{f}}\alpha_{f}\mathbf{v}_{n+1}-\mathbf{D}^{\intercal}_{n+\alpha_{f}}\alpha_{f}\>\mathbf{p}_{n+1}\\\\[8.61108pt]
&=\mathbf{f}_{n+\alpha_{f}}-\mathbf{M}_{n+\alpha_{m}}(1-\alpha_{m})\dot{\mathbf{v}}_{n}-\mathbf{K}_{n+\alpha_{f}}(1-\alpha_{f}){\mathbf{v}}_{n}+\mathbf{D}^{\intercal}_{n+\alpha_{f}}(1-\alpha_{f})\mathbf{p}_{n}\end{matrix}$
(47)

Repeating the procedure for the continuity equation and rearranging terms, the
following velocity-based system of equations is obtained:

$\begin{matrix}\left[\begin{matrix}\mathbf{M}_{{n+\alpha_{m}}}\>\displaystyle\frac{{\alpha_{m}}}{\gamma\Delta
t}+\mathbf{K}_{{n+\alpha_{f}}}{{\alpha_{f}}}&\>\>\>\>\>\>-\mathbf{D}^{\intercal}_{{n+\alpha_{f}}}{\alpha_{f}}\\\\[17.22217pt]
\mathbf{C}_{{n+\alpha_{m}}}\>\displaystyle\frac{\alpha_{m}}{\gamma\Delta
t}+\mathbf{D}_{{n+\alpha_{f}}}{\alpha_{f}}&\>\>\>\>\>\>\mathbf{L}_{{n+\alpha_{f}}}{\alpha_{f}}\end{matrix}\right]\;\left[\begin{matrix}\mathbf{v}_{n+1}\\\\[21.52771pt]
\mathbf{p}_{n+1}\end{matrix}\right]=\left[\begin{matrix}\mathbf{f}_{{n+\alpha_{f}}}+\mathbf{M}_{{n+\alpha_{m}}}\left(\displaystyle\frac{\alpha_{m}}{\gamma\Delta
t}\mathbf{v}_{\mathrm{n}}+\displaystyle\frac{{\alpha_{m}}-\gamma}{\gamma}\accentset{\mbox{\large.}}{\mathbf{v}}_{\mathrm{n}}\right){-(1-\alpha_{f})\left(\mathbf{K}_{n+\alpha_{f}}\mathbf{v}_{\mathrm{n}}-\mathbf{D}^{\intercal}_{n+\alpha_{f}}\mathbf{p}_{\mathrm{n}}\right)}\\\\[12.91663pt]
\mathbf{h}_{{n+\alpha_{f}}}+\mathbf{C}_{{n+\alpha_{m}}}\left(\displaystyle\frac{\alpha_{m}}{\gamma\Delta
t}\mathbf{v}_{\mathrm{n}}+\displaystyle\frac{{\alpha_{m}}-\gamma}{\gamma}\accentset{\mbox{\large.}}{\mathbf{v}}_{\mathrm{n}}\right){-(1-\alpha_{f})\left(\mathbf{D}_{n+\alpha_{f}}\mathbf{v}_{\mathrm{n}}+\mathbf{L}_{n+\alpha_{f}}\mathbf{p}_{\mathrm{n}}\right)}\end{matrix}\right]\end{matrix}$
(48)

The procedure for the displacement-based equation is omitted, although the
reader can refer to Appendix B (fourth row of table) to obtain the momentum
equation in the implementation approach GA-I${}^{q_{n+1}}$. Note that in both
approaches, GA-I${}^{q_{n+\alpha}}$ and GA-I${}^{q_{n+1}}$, and in both
formulations (velocity-based and displacement-based), matrices are computed at
$t_{n+\alpha}$, but GA-I${}^{q_{n+1}}$ defines the unknown state variables at
$t_{n+1}$. However, the change of variables in GA-I${}^{q_{n+1}}$ introduces
other terms on the right-hand side, in particular, the pressure term
$\mathrm{p}_{n}$, which is not present in the GA-I${}^{q_{n+\alpha}}$, Newmark
nor Backward Euler schemes.

##### GA-I${}_{p_{n+\alpha}}$ and GA-I${}_{p_{n+1}}$

As explained in the introduction of this manuscript, another implementation
approach commonly seen in the literature is to avoid the pressure at
$t_{n+\alpha_{f}}$ by writing it at $t_{n+1}$. Although placing the pressure
at $t_{n+1}$ is not consistent with the principle of the Generalized-$\alpha$
scheme, it is recognized as the most popular choice in the literature34. An
argument is that, as the pressure acts as a Lagrange multiplier, it should not
be subjected to time integration33. However, in the way in which time
integration schemes were presented in this section, such argument can be
questionable. It should be recalled that Backward Euler and Newmark have been
introduced as methods to approximate time derivatives of state variables
appearing in the momentum equation. However, the incompressible Navier-Stokes
equations do not invoke a time derivative of pressure, unless other
formulations are used such as the quasi-incompressible Newtonian fluid13. In
other words, per se, there is no need to apply time integration to the
pressure. On the other hand, the Generalized-$\alpha$ has been introduced as a
method that aims at improving Newmark’s accuracy by posing state equations at
intermediate times $t_{n+\alpha_{m}}$ and $t_{n+\alpha_{f}}$, but it does not
present a new approximation formula to the time derivatives of state
variables. So, rather than a time integration method, the Generalized-$\alpha$
method can be seen as a numerical scheme that upgrades Newmark’s accuracy from
first to second order.

The approach that places pressure at $t_{n+1}$ is likewise included in this
work for two reasons. First, for the sake of completeness since it represents
a popular choice in the literature. Second, depending on the implementation
approach, a pressure term ${\rm Phys.~{}Rev.~{}E}_{n}$ is obtained on the
right-hand side of the system of equations (compare GA-I${}^{q_{n+\alpha}}$
and GA-I${}^{q_{n+1}}$). This could result in spurious oscillations of the
external force vector due to numerical instabilities inherent to the spatial
discretization. If so, taking the pressure out of Generalized-$\alpha$ scheme
could prove beneficial in FEM formulations requiring stabilization.

The implementation approaches using pressure at $t_{n+1}$ and $t_{n+\alpha}$
are denoted by the subscripts $p_{n+1}$ and $p_{n+\alpha}$, respectively. In
combination with the other two approaches, four implementations of the
Generalized-$\alpha$ scheme in a velocity-pressure formulation are obtained:
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$. For a displacement-
pressure formulation, the combination reduces to three approaches:
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, and
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$, since
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ leads to the same
system of equations than
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$. Development of
equations for the cases imposing pressure at $t_{n+1}$ is omitted. Instead,
the reader is provided with the momentum equations in Appendix A (fifth and
sixth rows of table) for a velocity-based formulation, and in Appendix B
(fifth row of table) for a displacement-based formulation. Note that in such
tables, gradient matrix multiplying the pressure ${\rm
Phys.~{}Rev.~{}E}_{n+1}$ is integrated in the updated configuration
$\Omega(t_{n+1})$ to be consistent with the pressure time ($t_{n+1}$) obtained
by setting $\alpha_{f}=1.0$. That is, discretized momentum equations listed in
Appendices A and B for approaches
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ are obtained from the
following momentum equation:

$\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha_{m}}+\mathbf{K}_{n+\alpha_{f}}\mathbf{v}_{n+\alpha_{f}}-\mathbf{D}^{\intercal}_{\text{\lx@texthl@color{${n+1}$}}}\>\mathbf{p}_{n+1}=\mathbf{f}_{n+\alpha_{f}}$
(49)

An alternative implementation for
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ is to compute the
gradient matrix in the intermediate configuration $\Omega(t_{n+\alpha_{f}})$,
as done in Eulerian-based works 29. This alternative would lead to the
following momentum equation:

$\mathbf{M}_{n+\alpha_{m}}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+\alpha_{m}}+\mathbf{K}_{n+\alpha_{f}}\mathbf{v}_{n+\alpha_{f}}-\mathbf{D}^{\intercal}_{\text{\lx@texthl@color{${n+\alpha_{f}}$}}}\>\mathbf{p}_{n+1}=\mathbf{f}_{n+\alpha_{f}}$
(50)

where the highlighted term is to stand out the difference between Eqs. (49)
and (50). For the sake of simplicity, the time-discretized equations for Eq.
(50) are omitted, but they can be easily obtained from Table A1 (fifth and
sixth rows) and B2 (fifth row) by using a gradient matrix computed in
$\Omega(t_{n+\alpha_{f}})$.

In the following, equations are developed for the GA-II approach that assumes
linear scaling between times $t_{n}$ and $t_{n+1}$ for the dynamic, internal
and external forces.

####  3.3.2 GA-II

The formulas for computing the intermediate forces (Eq. 6) are replaced in the
momentum equation (Eq. 3), which leads to:

$\alpha_{m}\>\bm{f}^{\mathrm{dyn}}_{n+1}+\alpha_{f}\>\bm{f}^{\mathrm{int}}_{n+1}=\alpha_{f}\>\bm{f}^{\mathrm{ext}}_{n+1}+(1-\alpha_{f})\>\bm{f}^{\mathrm{ext}}_{n}-(1-\alpha_{m})\bm{f}^{\mathrm{dyn}}_{n}-(1-\alpha_{f})\bm{f}^{\mathrm{int}}_{n}$
(51)

Replacing terms for an incompressible Newtonian fluid, the following
discretized momentum equation is obtained:

$\begin{matrix}\mathbf{M}_{n+1}\alpha_{m}\accentset{\mbox{\large.}}{\mathbf{v}}_{n+1}+\mathbf{K}_{n+1}\alpha_{f}\mathbf{v}_{n+1}-\mathbf{D}^{\intercal}_{n+1}\alpha_{f}\>\mathbf{p}_{n+1}=\alpha_{f}\>\mathbf{f}_{n+1}+(1-\alpha_{f})\>(\mathbf{f}_{n}-\mathbf{K}_{n}\>{\mathbf{v}}_{n}+\mathbf{D}^{\intercal}_{n}\>\mathbf{p}_{n})-\mathbf{M}_{n}(1-\alpha_{m})\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\end{matrix}$
(52)

Then, Newmark formulas are used to get a velocity-pressure formulation.
Expressing the nodal acceleration in terms of the nodal velocities:

$\begin{matrix}\mathbf{M}_{n+1}\alpha_{m}\left(\displaystyle\frac{\mathbf{v}_{n+1}-\mathbf{v}_{n}}{\gamma\>\Delta
t}{-\frac{(1-\gamma)}{\gamma}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}}\right)+\mathbf{K}_{n+1}\alpha_{f}\mathbf{v}_{n+1}-\mathbf{D}^{\intercal}_{n+1}\alpha_{f}\>\mathbf{p}_{n+1}=\\\\[12.91663pt]
&\alpha_{f}\>\mathbf{f}_{n+1}+(1-\alpha_{f})\>(\mathbf{f}_{n}-\mathbf{K}_{n}\>{\mathbf{v}}_{n}+\mathbf{D}^{\intercal}_{n}\>\mathbf{p}_{n})-\mathbf{M}_{n}(1-\alpha_{m})\accentset{\mbox{\large.}}{\mathbf{v}}_{n}\end{matrix}$
(53)

Including the PSPG stabilization and arranging terms:

$\begin{matrix}\left[\begin{matrix}\mathbf{M}_{n+1}\>\displaystyle\frac{\alpha_{m}}{\gamma\Delta
t}+\mathbf{K}_{n+1}{\alpha_{f}}&\quad\quad-\mathbf{D}^{\intercal}_{n+1}\alpha_{f}\\\\[17.22217pt]
\mathbf{C}_{n+1}\>\displaystyle\frac{\alpha_{m}}{\gamma\Delta
t}+\mathbf{D}_{n+1}{\alpha_{f}}&\quad\quad\mathbf{L}_{n+1}\alpha_{f}\end{matrix}\right]\;\left[\begin{matrix}\mathbf{v}_{n+1}\\\\[21.52771pt]
\mathbf{p}_{n+1}\end{matrix}\right]=\\\\[34.44434pt]
&\left[\begin{matrix}\alpha_{f}\>\mathbf{f}_{n+1}+(1{\scriptstyle-}\alpha_{f})\>(\overbrace{\color[rgb]{0,0,1}\mathbf{f}_{n}-\mathbf{K}_{n}\mathbf{v}_{\mathrm{n}}+\mathbf{D}^{\intercal}_{n}\mathbf{p}_{\mathrm{n}}}^{\text{\normalsize{\color[rgb]{0,0,1}$\bm{f}^{\mathrm{ext}}_{n}-\bm{f}^{\mathrm{int}}_{n}$}}})-(1{\scriptstyle-}\alpha_{m})\overbrace{\color[rgb]{0,0,1}\mathbf{M}_{n}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}}^{\text{\normalsize{\color[rgb]{0,0,1}$\bm{f}^{\mathrm{dyn}}_{n}$}}}+\mathbf{M}_{n+1}\alpha_{m}\left(\frac{1}{\gamma\Delta
t}\mathbf{v}_{\mathrm{n}}+\frac{1-\gamma}{\gamma}\accentset{\mbox{\large.}}{\mathbf{v}}_{\mathrm{n}}\right)\\\\[12.91663pt]
\alpha_{f}\>\mathbf{h}_{n+1}+(1{\scriptstyle-}\alpha_{f})\>({\color[rgb]{0,0,1}\mathbf{h}_{n}-\mathbf{D}_{n}\mathbf{v}_{\mathrm{n}}-\mathbf{L}_{n}\mathbf{p}_{\mathrm{n}}})-(1{\scriptstyle-}\alpha_{m})\>\>{\color[rgb]{0,0,1}\mathbf{C}_{n}\accentset{\mbox{\large.}}{\mathbf{v}}_{n}}\>+\mathbf{C}_{n+1}\alpha_{m}\left(\frac{1}{\gamma\Delta
t}\mathbf{v}_{\mathrm{n}}+\frac{1-\gamma}{\gamma}\accentset{\mbox{\large.}}{\mathbf{v}}_{\mathrm{n}}\right)\end{matrix}\right]\end{matrix}$
(54)

Notably, Eq. (54) does not contain terms at $t_{n+\alpha}$, which is the most
distinctive feature with respect to the GA-I equations. The reader can easily
verify this by comparing the seventh row of Table A1 (appendix A) with the
previous rows of that table. At first glance, this could be seen as an
advantage for computing programming, since there is no need to calculate and
store intermediate state variables. However, in the PFEM context, this does
not present a major advantage over the GA-I approach. The reason is that the
finite element mesh is reconstructed at each time step, so dynamic, internal
and external forces from the previous time step must be recalculated in the
new mesh (blue terms in Eq. 54). This requires to assemble the $\mathbf{K}$,
$\mathbf{D}$, $\mathbf{L}$ and $\mathbf{C}$ matrices not only at $t_{n+1}$ but
also at $t_{n}$. If there is no remeshing, then forces at $t_{n}$ can be
stored to avoid additional matrix assembling.

As done for the GA-I approach, equations can be developed for a displacement-
pressure formulation, whose momentum equation can be found in the sixth row of
Table B2 (appendix B). In addition, equations resulting from setting the
pressure at $t_{n+1}$ instead of $t_{n+\alpha_{f}}$ are presented in the last
rows of Tables A1 and B2\. Those equations assume that gradient of pressure
${\rm Phys.~{}Rev.~{}E}_{n+1}$ is integrated in the current configuration, as
in Eq. (49). To obtain the expression that integrates the pressure gradient in
the intermediate configuration, as in Eq. (50), gradient matrix
$\mathbf{D}^{\intercal}_{n+1}$ must be replaced by
$\mathbf{D}^{\intercal}_{n+\alpha_{f}}$, which is computed as:

$\mathbf{D}^{\intercal}_{n+\alpha_{f}}=(1-\alpha_{f})\>\mathbf{D}^{\intercal}_{n}+\alpha_{f}\>\mathbf{D}^{\intercal}_{n+1}$
(55)

####  3.3.3 Generalized-$\alpha$ : Set of parameters

Being based on Newmark’s method, the Generalized-$\alpha$ method requires 4
parameters to be defined by the user: $\gamma$, $\beta$, $\alpha_{m}$, and
$\alpha_{f}$. However, these can be expressed as a function of the spectral
radius $\rho_{\infty}$. According to the stability analysis of Jansen et.al.23
developed for a linear system, the parameterization that ensures stability and
second-order accuracy is as follows:

$\alpha_{m}=\frac{1}{2}\left(\frac{3-\rho_{\infty}}{1+\rho_{\infty}}\right)\hskip
7.11317pt,\hskip 7.11317pt\alpha_{f}=\frac{1}{1+\rho_{\infty}}\hskip
7.11317pt,\hskip 7.11317pt\gamma=\frac{1}{2}+\alpha_{m}-\alpha_{f}\hskip
7.11317pt,\hskip 7.11317pt0\leq\rho_{\infty}\leq 1$ (56)

Note that formulas in Eq. (56) differ from those of Chung and Hulbert 22. On
one side because works define $\alpha$ in opposite ways ($\alpha$ weights
$t_{n+1}$ in Jansen.et.al.23 while $t_{n}$ in Chung and Hulbert22). Taking
this into account, expressions for $\alpha_{f}$ and $\gamma$ in Eq. (56) turn
out to be equivalent in both works. On the other hand, because amplification
matrices are different, since Chung and Hulbert22 consider displacement,
velocity and acceleration, while Jansen et.al.23 only accounts for velocity
and acceleration. This results in different equations for $\alpha_{m}$, from
which we chose that of Jansen et.al.23 since our momentum equation contains
only two kinematic variables.  Given the use of a Lagrangian scheme, the
$\beta$ variable from Newmark’s method must be included. As stated in Section
3.2, $\beta$ is parameterized according to Eq. (34). This would be one
distinctive aspect among Lagrangian and Eulerian implementations of the
Generalized-$\alpha$ method.

As highlighted by Jansen et.al.23 and Lovrić et.al.32, $\rho_{\infty}=0.0$
leads to a scheme with strong numerical damping for high frequencies. In terms
of spectral stability and for a linear problem, $\rho_{\infty}=0.0$ is
equivalent to the backward differentiation scheme of second order. To make an
analogy with the way time integration formulas have been presented, the second
order backward differentiation scheme is obtained by defining the approximate
acceleration function as:

$\accentset{\mbox{\large.}}{v}=\frac{1}{3}\>\frac{\mathrm{v}_{n}-\mathrm{v}_{n-1}}{\Delta
t}+\frac{2}{3}\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ (57)

which time-integrated results in the following (compare with equation 2.7 of
reference 32):

$\mathrm{v}_{n+1}=\mathrm{v}_{n}+\frac{1}{3}\>(\mathrm{v}_{n}-\mathrm{v}_{n-1})+\frac{2}{3}\>\Delta
t\>\accentset{\mbox{\large.}}{\mathrm{v}}_{n+1}$ (58)

On the other hand, $\rho_{\infty}=1.0$ preserves all high frequencies and it
is equivalent to the midpoint rule. For a linear problem, it is also
equivalent to the Trapezoidal rule obtained by setting $\gamma=0.5$ and
$\beta=0.25$ in Newmark’s method.

For the numerical examples of this manuscript, $\rho_{\infty}$ is imposed and
$\alpha_{m}$, $\alpha_{f}$ and $\gamma$ are obtained from Eq. (56), while
$\beta$ from Eq. (34). It is noteworthy that $\rho_{\infty}=0.0$ leads to
$\alpha_{m}=1.5$ and $\alpha_{f}=1.0$. In such a case, the GA-I and GA-II
approches are identical since external and internal forces are computed at the
current time, $\bm{f}^{\mathrm{int}}_{n+1}$ and $\bm{f}^{\mathrm{ext}}_{n+1}$,
in addition, dynamic forces in these approaches are equivalent when mass
conservation holds (as shown in Eq. 36). Furthermore, if $\rho_{\infty}=0.0$,
approaches that impose pressure at $t_{n+1}$ agree with those that impose it
at $t_{n+\alpha_{f}}$. Meaning, all of the aforementioned implementation
approaches of the Generalized-$\alpha$ method should produce the same result
when $\rho_{\infty}=0.0$.

Before reporting numerical comparisons of the different time integration
schemes, the following subsection briefly describes our PFEM implementation.

###  3.4 PFEM implementation

The PFEM implementation of this work follows the original idea of Idelsohn
et.al.8, which is based on a remeshing procedure using the Alpha Shape
algorithm45. Namely, the fluid is initially discretized using a cloud of
particles, as shown in Fig. 2(a). Then a Delaunay triangulation is applied to
generate a mesh on the convex hull of the cloud (Fig. 2(b)). Next, the Alpha
Shape algorithm is used to remove those elements whose circumcircle radius
exceeds a global characteristic size of the mesh scaled by a
$\alpha_{\mathrm{shape}}$ parameter. In this way, depending on the value of
$\alpha_{\mathrm{shape}}$, large or highly distorted elements are removed from
the triangulation, leaving a discretization of the fluid with a clear
definition of its boundaries and free surface, as shown in Fig. 2(c). This
mesh is used to discretize the Lagrangian Navier-Stokes equations, whose
solution gives the nodal velocity (or displacement) and pressure (Fig. 2(f)).
The new fluid position is obtained by the time integration scheme (Fig. 2(e)).
Finally, the mesh quality is checked to determine if re-meshing is necessary
(Fig. 2(d)) or if the triangulation can be preserved for the next time step. A
recent survey on PFEM and more details on the above steps can be found in 7.

((a))

((b))

((c))

((d))

((e))

((f))

(a) (b) (c)
(d) (e) (f)

Figure 2: Illustration of the PFEM procedure using the collapse of a
cylindrical water column problem. (a) Cloud of particles with known state
variables at $t_{n}$. (b) Delaunay triangulation on the entire set of
particles. (c) Large and distorted elements are removed from the triangulation
using Alpha Shape algorithm45. (f) Norm of nodal velocities for illustrating a
FEM. (e) Illustration of updating state variables and nodal position. (d) The
updated cloud of particles to be used in the next time step.

The nonlinear system of equations, i.e. Eq. (33), (41), (44) or (48), is
assembled in matrix form as $\mathbf{A}\>\mathbf{q}=\mathbf{b}$, where
$\mathbf{A}$ is a large matrix storing matrix terms of the momentum and
continuity equations computed at time $t_{n+1}$, $\mathbf{q}$ is a vector
containing the solution the for state variables, and $\mathbf{b}$ is a vector
storing the external forces. Note that all mass matrices are "consistent" and
computed at the time specified by the subscript ($t_{n}$, $t_{n+1}$ or
$t_{n+\alpha_{m}}$). The system $\mathbf{A}\>\mathbf{q}=\mathbf{b}$ is solved
using the Picard (or fixed-point) algorithm, which sets as predictor the
values of state variables of previous time $t_{n}$. Regarding the
stabilization term, $\tau_{\mathrm{PSPG}}$ is computed using the velocity of
the previous iteration of the nonlinear algorithm. Once the solution of
$\mathbf{q}$ is obtained, particle kinematics is updated according to the time
integration scheme. Updating equations are summarized in Appendix C.

The success of the PFEM method relies on the use of efficient remeshing
algorithms, such as the Delaunay triangulation followed by the Alpha Shape.
This process involves a negligible computational cost in 2D and a minor one in
3D. However, the Alpha Shape algorithm brings some issues during remeshing
operations. The most notable is the mass variation of the system produced by
the discrete addition and removal of finite elements, which can be treated
with mesh refinement or tuning of the $\alpha_{\mathrm{shape}}$ parameter46.
However, examples in the present work are free of such a problem either
because they do not require remeshing (such as the fluid sloshing), possess an
Eulerian domain (such as the flow around a cylinder), or because they have
slip boundary conditions and do not feature fluid splashing (such as the
solitary wave propagation and collapse of a water column). The parameter
$\alpha_{\mathrm{shape}}$ is set to 1.25 in all examples of this work.

PFEM can handle problems with non-uniform mesh size and Eulerian domains, such
as the flow around a cylinder problem introduced in the following section.
Mesh refinement is achieved by adding and removing particles into the domain
according to a size field. For the flow around a cylinder problem (Section
4.3), the size field is a function of the distance to the cylinder, while for
the collapse of a water column problem (Section 4.4), the size field is
defined in terms of the distance to the free surface. If the size of a finite
element is larger than the prescribed size, a particle is added in the center
of the element and mesh refinement occurs. The un-refinement process operates
in a similar way, but now a particle is removed if the distance to the nearest
particle is less than that prescribed in the size field. This addition and
removal process imposes upper and lower size limits on the elements depending
on the desired size distribution. An illustration of this process for the flow
around a cylinder problem is shown in Fig. 3. In this context, the Alpha Shape
must be applied locally, i.e. the upper bound of the Alpha-Shape algorithm
must be defined by a parameter $\alpha_{\mathrm{shape}}$ scaled by the desired
local element size and not by a global reference size as previously mentioned.

To simulate an Eulerian domain (the flow around a cylinder), particles in the
inlet and outlet zone are treated in a Lagrangian fashion, i.e., they are
allowed to move. However, before proceeding with the remeshing process for the
next time step, these particles are repositioned on the inlet and outlet
lines. For further details regarding the simulation of Eulerian domains using
PFEM, readers are referred to 47, 48.

Importantly, the remeshing process in PFEM does not affect the use of the
Generalized-$\alpha$ method. Therefore, implementation aspects of the
Generalized-$\alpha$ presented for the standard FEM are equally valid in PFEM.
The only practical shortcoming is with respect to the GA-II approach, as it
requires to recompute the dynamic, internal and external forces of the
previous time $t_{n}$ due to the new nodal connectivity at time $t_{n+1}$.
Likewise, drawbacks could arise with material models that require historical
variables, however such cases are outside the scope of this work.

Figure 3: Illustration of particle addition and removal for mesh refinement in
PFEM. The flow around a cylinder problem is illustrated. On the left side of
the cylinder, particles are added in the center of elements as they approach
the cylinder. On the right side, particles are removed if the distance to the
nearest particle is less than an imposed threshold.

The time integration schemes analyzed in this work have been implemented in an
in-house code named PFEM3D 40. This is a C++ code that uses the Eigen49
library for linear algebra operations, the CGAL50 library for triangulation
and Alpha Shape, OpenMP51 for parallel threading, and Gmsh52 for the
input/output of finite element meshes.

##  4 Numerical Examples

Four distinctive benchmark problems are used in this section for comparing
time integration schemes. The first is the sloshing of a fluid characterized
by small displacements so no remeshing is needed, which allows to validate the
FEM implementation of the time integration schemes. The second is the solitary
wave propagation in a constant depth container, which uses remeshing for minor
improvements to the mesh discretizing the wave. The third is the flow around a
cylinder, characterized by an Eulerian computational domain that features mesh
refinement and significant velocity gradients around the cylinder, so that
remeshing is applied at each time step. The fourth is the collapse of a
cylindrical water column which, unlike the previous ones, is a 3D problem
where time step $\Delta t$ is defined in terms of a maximum CFL
(Courant–Friedrichs–Lewy) number.

This section mainly reports simulation data for comparison effects and does
not present snapshots of simulations. Instead, the reader is guided to
reference 53 for animations of each problem.

###  4.1 Small amplitude sloshing

This free-surface problem consists of simulating the sloshing of a fluid
caused by imposing a non-uniform initial elevation of the free surface 54. The
geometry of the problem is illustrated in Fig. 4(a). The width of the
container is $b=1$ m and the initial free surface elevation is defined as:

$h(x)=H+a\>\mathrm{sin}\left(\frac{\pi(x-x_{0})}{b}-\frac{\pi}{2}\right)$ (59)

where $H$ is the mean height equal to 1.0 m and $a$ is the initial amplitude
equal to 0.01 m. A unit gravity acceleration is used, $g$ = 1.0 m/s2, as well
as fictitious fluid parameters, $\mu=0.01$ Pa s and $\rho$ = 1.0 kg/m3. The
fluid domain is discretized with 3073 particles and 5942 elements of 0.02 m
average size, as shown in Fig. 4(b). No external pressure is imposed in the
free-surface, a free-slip condition is defined on the container walls and non-
slip condition in the base. Given the small initial amplitude $a$, elements
are subjected to small deformations so that remeshing process is not needed in
this problem.

((a))

((b)) Figure 4: Small amplitude sloshing. (a) Geometry illustration. (b)
Initial finite element discretization

(a) (b)

Several sloshing problems are solved with variations in the time step $\Delta
t$, which are 0.400, 0.100, 0.025, and 0.00625 s. The criterion for comparison
is free surface elevation at coordinate $x_{b}$ (see Fig. 4(a)). Fig. 5(b)
plots $h(x_{b})$ versus time for different time steps using Backward Euler
(Fig. 5(a)) and Generalized-$\alpha$ (Fig. 5(b)) in a velocity-pressure
formulation. Implementation
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ with a
spectral radius $\rho_{\infty}=0.0$ is used for Generalized-$\alpha$. Graphs
also show a solution obtained with the Trapezoidal rule ($\Delta t=0.00625$
s). For validating the simulation, extreme points of the oscillatory curve
reported by Avancini and Sanches 36 are also displayed.

The strong dependence of Backward Euler on the time step is clear when
comparing Figs. 5(a) and 5(b). Even when using a small time step $\Delta
t=0.00625$ it is not possible to match the curve of the Trapezoidal rule.
Instead, the Generalized-$\alpha$ method produces curves very close to that of
the Trapezoidal rule, even when using a large $\Delta t=0.100$. To facilitate
the comparison between time integration schemes, an error measure is plotted
next. The error criterion considers the area comprised between a reference
elevation curve $h_{\mathrm{ref}}(t)$ and that of the simulation $h(x_{b},t)$,
that is:

$A_{\mathrm{diff}}=\int_{0}^{20}\left|h_{\mathrm{ref}}(t)-h(x_{b},t)\right|dt$
(60)

where $A_{\mathrm{diff}}$ is the error measure. The reference elevation
$h_{\mathrm{ref}}(t)$ represents a solution obtained with a time step
$\Delta_{t}=0.0015625$ using the Trapezoidal rule.

((a))

((b)) Figure 5: Results of the fluid sloshing. Elevation of free surface with
respect to time. Results obtained using the Trapezoidal rule, (a) Backward
Euler and (b) Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ and
$\rho_{\infty}=0.0$).

(a) (b)

((a))

((b)) Figure 6: Results of the fluid sloshing. Error between reference
elevation $h_{\mathrm{ref}}$ and simulated elevation $h$ at coordinate
$x_{b}$. Simulations obtained using (a) Backward Euler, Newmark’s Method (NM),
Trapezoidal rule, and (b) Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$).

(a) (b)

Now, the Backward Euler, Newmark, Generalized-$\alpha$ and Trapezoidal rule
are considered in a velocity-based formulation. For Newmark, 6 problems are
solved that differ in the value of $\gamma$, which are set to 0.9, 0.8, 0.7,
0.6, 0.55 and 0.52. The Generalized-$\alpha$ method considers 10 values of
$\rho_{\infty}$ ranging from 0.0 to 0.9 with an increment of 0.1. Each time
integration scheme is solved for 7 time steps $\Delta t$, ranging from 0.00625
to 0.4 s. That is, 126 problems are solved in total using different time
integration schemes and time steps. Each one is compared against the reference
curve to calculate the $A_{\mathrm{diff}}$ error. Results are shown in Figs.
6(a) and 6(b). Performance of time integration schemes is consistent with the
literature. Among all the schemes, Backward Euler exhibits the lowest accuracy
in the free surface height (black curve in Fig. 6(a)). Accuracy is improved by
reducing time step $\Delta t$ but with an order of convergence approaching 1.
The Trapezoidal rule is the most accurate and its error decreases with second
order with respect to time step (red curve in Fig. 6(a)). Newmark’s Method
(NM) interpolates both Backward Euler and Trapezoidal schemes. However,
$\gamma$ must be chosen very close to 0.5 to exhibit second-order convergence.
In the case of Generalized-$\alpha$, all curves exhibit second order
convergence for large time steps (Fig. 6(b)). This expected result is due to
the parameterization of $\gamma$, $\beta$, $\alpha_{m}$ and $\alpha_{f}$ in
terms of the spectral radius $\rho_{\infty}$, which ensures stability and
second-order accuracy for a linear problem.

After validating our implementations of time integration schemes in PFEM with
the results of Avancini and Sanches 36, and obtaining an anticipated behavior
of these for a problem with small displacements and deformations, the
following examples compare the different implementation approaches of the
Generalized-$\alpha$ method and incorporate the displacement-based
formulation.

###  4.2 Solitary wave propagation

The geometry of the free-surface problem is illustrated in Fig. 7(a). This
consists of simulating a wave of amplitud $A$ propagating in a rectangular
container of constant depth $H$. Instead of displacing a fluid volume at the
beginning of the simulation to generate the wave, another common approach 35,
36, 54, 55 is adopted here: to impose as initial condition the analytical
approximation proposed Wehausen and Laitone 56. Taking their system of
equations with second order approximation (equation 31.37 of reference 56),
the free surface height, velocity and pressure are defined as follows:

$\displaystyle h(x,t)$
$\displaystyle=\mathrm{H}+A\>\mathrm{sech}^{2}(\xi)-\frac{3A^{2}}{4H}\mathrm{sech}^{2}(\xi)\left(1-\mathrm{sech}^{2}(\xi)\right)$
(61a) $\displaystyle\frac{\mathrm{v}_{x}(x,y,t)}{\sqrt{gH}}$
$\displaystyle=\frac{h-H}{H}-\frac{A^{2}}{2H^{2}}\left(1+6\>\frac{y-H}{H}+3\>\frac{(y-H)^{2}}{H^{2}}\right)\>\mathrm{sech}^{2}(\xi)+\frac{A^{2}}{2H^{2}}\left(1+9\>\frac{y-H}{H}+9\>\frac{(y-H)^{2}}{2H^{2}}\right)\mathrm{sech}^{4}(\xi)$
(61b) $\displaystyle\frac{\mathrm{v}_{y}(x,y,t)}{\sqrt{gH}}$
$\displaystyle=\frac{y\sqrt{3}A^{3/2}}{H^{5/2}}\mathrm{sech}^{2}(\xi)\mathrm{tanh}(\xi)\left[1+\frac{A}{H}\left(\frac{1}{8}-\frac{y}{H}-\frac{(y-H)^{2}}{2H^{2}}\right)-\frac{A}{\mathrm{2H}}\left(1-6\>\frac{y-H}{H}-3\>\frac{(y-H)^{2}}{H^{2}}\right)\mathrm{sech}^{2}(\xi)\right]$
(61c) $\displaystyle\frac{{\rm Phys.~{}Rev.~{}E}(x,y,t)}{\rho gH}$
$\displaystyle=1+\frac{h-H-y}{H}-\frac{3A^{2}}{4\mathrm{H^{2}}}\left(\frac{y^{2}}{H^{2}}-1\right)\left(2\>\mathrm{sech}^{2}(\xi)-3\>\mathrm{sech}^{4}(\xi)\right)$
(61d)

where the time-dependent parameter $\xi$ is defined as:

$\xi(x,t)=\frac{x-c\>t}{H}\left(\frac{3A}{4H}\right)^{1/2}\left(1-\frac{5A}{8H}\right)$
(62)

and the wave celerity given by:

$\frac{c}{\sqrt{gH}}=1+\frac{A}{2H}-\frac{3A^{2}}{20H^{2}}$ (63)

((a))

((b)) Figure 7: Solitary wave propagation in constant water depth. (a)
Geometry illustration and (b) initial Finite element discretization around the
wave. Subfigure (b) represents a section of the model that is indicated by
dashed lines in (a).

(a) (b)

Due to the position of the coordinate system (see Fig. 7(a)), initial
condition for velocity, pressure and free surface height is obtained by
setting $t=0$ in Eq. (62). This produces a wave displacement to the right,
whose position in time is defined by Eq. (61a). The chosen fluid has viscosity
$\mu$ = 0.001 Pa s and density $\rho$ = 1000 kg/m3. The container depth is $H$
= 10 m, the amplitude $A$ is 2 m and the gravity acceleration is $g=9.81$m/s2.
The initial domain is discretized using 2400 particles and 4090 finite
elements of 2 m average size, as shown in Fig. 7(b).

((a))

((b)) Figure 8: Wave profile after 40 s using (a) Backward Euler and (b)
Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$) and
different time steps.

(a) (b)

##### Performance of time integration schemes

Wave profiles after 40 s of simulation are shown in Figs. 8(a) and 8(b) using
Backward Euler and Generalized-$\alpha$, respectively. Both in a velocity-
based formulation. Again, results show a strong dependence of Backward Euler
on the time step $\Delta t$, whose excessive numerical damping for large time
steps results in a decrease of wave amplitude. Figs. 10(a) and 8(b) also
include the analytical approximation of Eqs. (61)-(63), which is in good
agreement with the simulations. However, the maximum height does not
completely match. This is because analytical equations are limited to linear
conditions, such as waves with small amplitude $A$. As shown in 55 and 36, the
analytical approximations are rather accurate if $A/H\leq$ 0.1, but less if
$A/H\geq$ 0.2. Therefore, for error analysis of time integration schemes, the
maximum wave height at 40 seconds of simulation is compared with that obtained
from a simulation using a time step $\Delta t=0.00078125$, the
Generalized-$\alpha$ method $(\rho_{\infty}=0.0)$ and a discretization with
finite elements of 1 m size. The obtained reference height is
$h_{\mathrm{ref}}=$ 12.0255.

((a)) ((b)) Figure 9: Maximum height of the solitary wave normalized with
respect to the reference height. (a) Curves obtained with Backward Euler,
Trapezoidal, Newmark and Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$). Subfigure
(b) is an enlarged view of (a). Numbers next to black curves (dashed and
solid) indicate the order of convergence of these curves.

(a) (b)

The maximum wave height divided by the reference height is plotted next. The
wave height is obtained from simulations using different time integration
schemes and time steps. Specifically, seven time steps $\Delta t$ are
considered, ranging from 0.003125 to 0.2. The Backward Euler, Trapezoidal,
Newmark and Generalized-$\alpha$ schemes are considered. For Newmark, four
combinations of parameters are chosen from defining $\gamma$ as 0.9, 0.8, 0.7
and 0.6. For Generalized-$\alpha$, implementation approach
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ is used and
ten combinations of parameters are chosen from defining $\rho_{\infty}$ as 0.0
to 0.9 with an increment of 0.1. Thus, 112 solitary wave propagation problems
are solved in total. Results are summarized in Fig. 9(a), which show a pattern
similar to the previous example. Backward Euler exhibits the lowest wave
height accuracy, which improves as time step is reduced but with an order of
convergence close to 1.0 (as shown by numbers next to the black dashed curve
in Fig. 9(a)). The Trapezoidal scheme is the most accurate while Newmark’s
Method (NM) interpolates Backward Euler and Trapezoidal accuracy. On the other
hand, the Generalized-$\alpha$ method shows wave heights very close to the
reference value, which gathers the GA-I curves in the upper zone of Fig. 9(a).
For better visualization, such a zone is enlarged and is shown in Fig. 9(b).
Numbers next to the black curve indicate the order of convergence for
$\rho_{\infty}$ = 0.0. It can be seen that for large time steps ($\Delta t\geq
0.05$), the maximum wave height converges with 2nd order with respect to
$\Delta t$. For smaller time steps, the order of convergence is reduced. This
is attributed to a spatial discretization error resulting from a coarse
discretization, which was chosen to extend the range of values of $\Delta t$
without compromising convergence of the nonlinear algorithm. Space
discretization error becomes predominant as the maximum height of the
simulation approaches the reference value, which is reflected in the
Generalized-$\alpha$ and Trapezoidal curves when $\Delta t\leq 0.05$.

##### Implementation approaches of the Generalized-$\alpha$ method

The next analysis incorporates results obtained with
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ and $\text{GA-
II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, i.e., it compares the
Generalized-$\alpha$ implementation approaches that write the pressure at
$t_{n+\alpha_{f}}$. The same 7 time step sizes and 10 values of
$\rho_{\infty}$ of previous analysis are considered. Results are summarized in
Fig. 10(a). Color scale indicates the value of $\rho_{\infty}$ while marker
type distinguishes implementation schemes. The relative difference between the
3 schemes is shown in Fig. 10(b). There, percentage difference is computed
with respect to the scheme
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, taken as a
reference.

It is observed that for $\rho_{\infty}$ = 0.0, all three schemes produce the
same result when $\Delta t\leq 0.1$ (difference of $\approx$ 0$\%$ in Fig.
10(b), black color). In addition, the
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ and
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ schemes are
coincident for all $\rho_{\infty}$ (all circle markers are at 0$\%$ in Fig.
10(b)). Thus, at least in this problem, the algebraic treatment of the linear
interpolation of state variables (Eq. 4) is not relevant to the results. In
case of $\text{GA-II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, a difference with
respect to GA-I is observed. The biggest difference is observed for $\Delta t$
= 0.2 and $\rho_{\infty}\geq 0.7$, which is $\approx$ 0.007 $\%$. The
difference between GA-I and GA-II reduces rapidly when decreasing
$\rho_{\infty}$ or the time step $\Delta t$, and becomes negligible for
$\Delta t\leq 0.05$. This observation is in line with previous comments of
section (Eq. 37), where it is stipulated that GA-I and GA-II should converge
to the same result in the presence of small deformations or small time steps
if mass conservation holds.

((a))

((b)) Figure 10: (a) Normalized maximum height of the solitary wave using
different implementation approaches of the Generalized-$\alpha$ method. (b)
Relative difference between approaches, with
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ used as
reference.

(a) (b)

Next, results from implementation approaches of the Generalized-$\alpha$ that
write the pressure at $t_{n+1}$ instead of $t_{n+\alpha_{f}}$ are presented.
The three implementation alternatives are considered,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$. In addition, two scenarios are taken into
account, one that integrates the gradient matrix at $t_{n+1}$ (as in Eq. 49)
and another that integrates it at $t_{n+\alpha}$ (as in Eq. Eq. 50). Results
are shown in Figs. 11(a) and 11(b), respectively (plots appear on different
scales). Figures also include scheme
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ that follows
the rationale of the generalized alpha scheme and writes the pressure at
$t_{n+\alpha_{f}}$.

Results reveal that there is no significant difference between the GA-I and
GA-II implementation approaches that write the pressure in $t_{n+1}$, since in
both figures, $\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$ are coincident. However, there is a big
difference for when the gradient matrix is integrated. If
$\mathbf{D}^{\intercal}_{n+1}$ is used (Fig. 11(a)), then performance of
Generalized-$\alpha$ is reduced, especially for large values of
$\rho_{\infty}$ since performance becomes similar to that of Backward Euler.
In contrast, if $\mathbf{D}_{n+\alpha_{f}}^{\intercal}$ is used (Fig. 11(b)),
then no detriment is seen on the time integration scheme, regardless of the
value of $\rho_{\infty}$ and the implementation approach
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ or $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$). Remarkably, all implementation schemes
lead to the same result if $\rho_{\infty}$ = 0.

((a))

((b)) Figure 11: Normalized maximum height of the solitary wave. Subfigures
compare implementation approaches of the Generalized-$\alpha$ method that
write pressure at $t_{n+1}$. Schemes
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$ compute the gradient matrix
$\mathbf{D}^{\intercal}$ in (a) $\Omega(t_{n+1})$ and (b)
$\Omega(t_{n+\alpha_{f}})$.

(a) (b)

##### Velocity-based and Displacement-based formulations

A final analysis using the solitary wave propagation problem compares results
between velocity-based and displacement-based formulations. This is conducted
for Backward Euler, Newmark and Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$). To avoid
overloading graph results, only two sets of parameters are used for Newmark
($\gamma$ equal to 0.8 and 0.6), and five for Generalized-$\alpha$
($\rho_{\infty}$ equal to 0.0, 0.2, 0.4, 0.6, and 0.8). Results are given in
Fig. 12. It is observed that both formulations lead to practically the same
results. That is, the change of unknown variable in the system of equations
does not seem to have an impact on the accuracy of results. Furthermore, we
did not observe significant changes in computation time between the two
formulations (although computation time was not precisely tracked in the
simulations) nor in the number of iterations used by the nonlinear algorithm.
Therefore, besides adding additional terms on the right-hand side of the
system of equations, it seems that there is no significant computational
impact from switching to a displacement-pressure formulation.

Observations up to this point have been obtained from two problems that
feature small deformations and exhibit quasi linear behavior. Furthermore,
they are analyzed from a kinematic point of view while pressure has been
excluded from the analysis. Therefore, the following example of the flow
around a cylinder is chosen due to the presence of large velocity gradients
and because comparison criterion is based on the force exerted on the
cylinder, which is highly dependent on the pressure field.

Figure 12: Normalized maximum wave height using Backward Euler, Newmark and
Generalized-$\alpha$ in velocity-based and displacement-based formulations.

###  4.3 Flow around a cylinder

This example involves the flow past a circular cylinder. The computational
domain is designed as in Idelsohn et.al57, 58. Therefore, fluid particles
cross a rectangular computational domain, as illustrated in Fig. 13(a). A
horizontal velocity $\mathrm{v}_{0}$ is imposed on particles located at the
left, top and bottom boundaries. Particles that cross the right boundary are
deleted from the model. A rigid cylinder of 1 m diameter is spaced 5 m from
the left, top and bottom boundary, which leaves the cylinder at 15 m from the
right boundary.

((a))

((b)) Figure 13: Flow around a cylinder. (a) Geometry and (b) initial finite
element discretization.

(a) (b)

Results for two Reynolds (Re) numbers are reported in this example, 100 and
1000. The first is included in the study of Dettmer et.al.24, who solve the
flow around a cylinder problem using the Generalized-$\alpha$ scheme with
various values of $\rho_{\infty}$ in an Eulerian framework. The second, Re =
1000, is included for comparing results with Idelsohn et al.57, 58 who are, to
date, the only ones to present results for the flow around a cylinder using
PFEM, to the best of our knowledge.

Reynolds numbers 100 and 1000 are imposed by defining the viscosity $\mu$ as
0.01 and 0.001, respectively. Whereas, fluid density and inlet velocity are
set to $\rho=1$ kg/m3 and $\mathrm{v}_{0}=1$ m/s, respectively. Problems with
Reynolds 100 and 1000 are simulated for 160 and 90 s of physical time,
respectively. For both problems, a non-uniform mesh size is used, which is
defined by the following function:

$h_{\mathrm{elem}}=\mathrm{min}\left[\left(0.05\>(x+5)^{2}+0.05\>y^{2}+0.01\right)\>,\>\left(0.001\>(x-1)^{2}+0.02*y^{2}+0.12\right)\right]$
(64)

where $h_{\mathrm{elem}}$ is the element size at coordinates $x$ and $y$. The
adopted mesh generator (Gmsh52) produces 8018 particles with Eq. (64), places
140 on the cylinder surface and generates 15777 finite elements in the initial
computational domain, as shown in Fig. 13(b). During PFEM simulation, each
remeshing process ensures that spatial distribution of particles agrees with
Eq. (64), otherwise particles are added or removed.

For this example, the Strouhal number, the amplitude of the lift coefficient
and the mean value of the drag coefficient are reported. These are
respectively denoted as $S_{t}$, Amp($C_{L}$) and Mean($C_{D}$). The Strouhal
number ($S_{t}$), lift ($C_{L}$) and drag ($C_{D}$) coefficients are defined
as:

$S_{t}=\frac{f_{L}\>D}{\mathrm{v}_{0}}\;\;\;\;,\;\;\;\;\;C_{L}(t)=\frac{F_{\mathrm{lift}}(t)}{\frac{1}{2}\rho\>\mathrm{v}_{0}^{2}\>D}\;\;\;\;,\;\;\;\;\;C_{D}(t)=\frac{F_{\mathrm{drag}}(t)}{\frac{1}{2}\rho\>\mathrm{v}_{0}^{2}\>D}$
(65)

where $F_{\mathrm{lift}}$ and $F_{\mathrm{drag}}$ are the lift and drag
forces, $D$ is the cylinder diameter (equal to 1 m) and $f_{L}$ is the
frequency of the lift force. Typical lift and drag coefficients curves, for
problems with Reynolds 100, are shown in Figs. 14(a) and 14(b), respectively.
Curves for problems with Reynolds 1000 are shown in Figs. 14(c) and 14(d). The
frequency of the lift force $f_{L}$, equivalent to the frequency of the lift
coefficient, is obtained using a Fast Fourier Transform (FFT) on signal
segments shown in Figs. 14(a) and 14(c).

((a))

((b))

((c))

((d)) Figure 14: Simulation of the flow around a cylinder for Reynolds 100
(subfigures a,b) and 1000 (subfigures c,d). Evolution of the lift (subfigures
a,c) and drag (subfigures b,d) coefficients. Results obtained with Backward
Euler and $\Delta t=0.016$.

(a) (b)
(c) (d)

##### Verification against the literature

Several problems of the flow around a cylinder with variations in the time
integration scheme and time step are solved. Specifically, for problems with
Re = 100, six time steps are considered ranging from 0.002 to 0.064. Backward
Euler, Trapezoidal, Newmark and Generalized-$\alpha$ are used as time
integration schemes. For Newmark, four combinations of parameters are used by
defining $\gamma$ as 0.9, 0.8, 0.7 and 0.6. For Generalized-$\alpha$,
implementation $\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$
is used and six values of $\rho_{\infty}$, which are 0.0, 0.1, 0.2, 0.5, 0.7,
and 0.9. Therefore, 72 problems of the flow around a cylinder are solved for
Re = 100. The coefficients $S_{t}$, Amp($C_{L}$) and Mean($C_{D}$) are
summarized in Fig. 15(f). On the left side (Figs. 15(a), 15(c) and 15(e)) are
results with Backward Euler, Newmark’s Method (NM), and Trapezoidal rule,
while on the right side (Figs. 15(b), 15(d) and 15(f)) are those obtained with
Generalized-$\alpha$. In addition, some references are included in the graphs,
which appear as horizontal lines for ease of comparison, but it does not
indicate that references are insensitive to time step $\Delta t$.

((a))

((b))

((c))

((d))

((e))

((f)) Figure 15: Coefficients of the flow around a cylinder at Re = 100. (a,b)
Strouhal number (c,d) mean value of the drag coefficient, and (e,f) amplitude
of the lift coefficient. Results obtained with (a,c,d) Backward Euler,
Newmark’s Method (NM) and Trapezoidal rule, and (b,d,f) Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$).

(a) (b)
(c) (d)
(e) (f)

((a))

((b))

((c))

((d))

((e))

((f)) Figure 16: Coefficients of the flow around a cylinder at Re = 1000.
(a,b) Strouhal number (c,d) mean value of the drag coefficient, and (e,f)
amplitude of the lift coefficient. Results obtained with (a,c,d) Backward
Euler, Newmark’s Method (NM) and Trapezoidal rule, and (b,d,f)
Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$).

(a) (b)
(c) (d)
(e) (f)

For problems with Re = 1000, five time steps are used ranging from 0.002 to
0.032. Similarly as for Re = 100, Backward Euler, Trapezoidal, Newmark and
Generalized-$\alpha$ are used. For Newmark, $\gamma$ is set to 0.9, 0.8, 0.7
and 0.6. For Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$),
$\rho_{\infty}$ is set to 0.0, 0.1, 0.2, 0.5. Thus, 50 problems are solved in
total for Re = 1000. Results are summarized in Fig. 16(f) using same layout as
for Re = 100 (Fig. 15(f)).

Results in Figs. 15(f) and 16(f) show that $S_{t}$, Amp($C_{L}$) and
Mean($C_{D}$) obtained with stable time integration schemes converge to the
same value as time step $\Delta t$ is reduced. Converged values are listed in
Table 1. These are obtained from averaging coefficients obtained with $\Delta
t$ = 0.002. Table 1 also includes values from references that use Eulerian24,
32 or semi-Lagrangian approaches (PFEM-2) 57, 58. Notably, coefficients
obtained in this work are in agreement with the literature.

##### Performance of time integration schemes

As in previous examples, Backward Euler is the least accurate at large time
steps $\Delta t$. Its excessive numerical damping for high frequencies reduces
both the magnitude and frequency of the lift force. In contrast, Newmark’s
Method (NM) offers better accuracy as $\gamma$ is reduced, as long as $\gamma$
and $\beta$ do not match the Trapezoidal rule, which leads to some problems in
the present case. Noteworthy, the Trapezoidal curve (red) does not completely
appear in Figs. 15(f) and 16(f) because either values are far outside the
range of the plot (values for $\Delta t$ < 0.008) or because convergence was
not achieved with the imposed time step ($\Delta t$ > 0.016). The instability
of the trapezoidal rule in the absence of numerical damping at high
frequencies is illustrated in Fig. 17(a). This figure depicts 4 cycles of the
lift coefficient (for Re = 100) obtained with different time integration
schemes and $\Delta t=0.004$. The curves are aligned on the left side and
abscissa is tabulated with respect to the converged frequency (equal to the
converged $S_{t}$ number). Due to the lack of numerical damping, the lift
coefficient shows huge values in the Trapezoidal scheme, thus it had to be
scaled by 1/25 to make it fit into Fig. 17(a).

Regarding the Generalized-$\alpha$ scheme, $\rho_{\infty}\leq 0.5$ exhibits
high accuracy with respect to lift and drag coefficients, even for large time
steps $\Delta t$. As numerical damping for high frequencies is reduced,
specifically for $\rho_{\infty}>0.5$, accuracy is lost because spurious
oscillations become dominant in the system. This can be seen in Fig. 17(a) and
17(b).

((a)) Obtained with $\Delta t=0.004$.

((b)) Obtained with $\Delta t=0.064$.

Figure 17: Flow around a cylinder at Re = 100. Four cycles from the lift
coefficient curve obtained with Backward Euler, Newmark’s Method (NM),
Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$) and
Trapezoidal rule. Abscissa is scaled with respect to the reference Strouhal
number $S_{t}=0.174$. In (a), lift coefficient curves of GA-I
($\rho_{\infty}=0.90$) and Trapezoidal rule are scaled by 1/2.5 and 1/25,
respectively.

A remarkable aspect of the Generalized-$\alpha$ method is that results are
very similar for $\rho_{\infty}\leq 0.5$. This holds for both Reynolds (100
and 1000) and for both forces (lift and drag), as shown in Figs. 15(b), 15(d),
15(f), 16(b), 16(d), 16(f). On one hand, this is in line with observations of
Jansen et.al.23 who state "… the period and amplitude of both the lift and the
drag are very weak functions of $\rho_{\infty}$ (which might be expected from
the observation that 30 points per wavelength is adequate for second-order
accurate method)". On the other hand, it is in contradiction with results by
Dettmer and Perić 24, who report a notorious influence of the lift force with
respect to the spectral radius $\rho_{\infty}$. A possible explanation is that
time step used in the present work and in that of Jansen et.al.23 is
significantly smaller than the one used by Dettmer and Perić 24. At the
largest time step ($\Delta t=0.064$), the period of the lift force is
discretized in $\approx$90 points in our work, in $\approx$60 points in the
work of Jansen et.al.23 (with $\Delta t=0.1$), while in $\approx$18 points in
that of Dettmer and Perić 24 (with $\Delta t=0.3$). In the latter case, a
coarse time discretization could be standing out a dependence of results on
$\rho_{\infty}$.

Time steps used in this work are much smaller than those used by reference24
because a finer mesh around the cylinder is used here in order to reduce
instabilities of the pressure gradient caused by the spatial discretization
and the remeshing. If a coarser mesh were used, the lift and drag force would
contain bigger spurious oscillations (see figure 10.11 of reference 57) and it
would be difficult to identify numerical damping associated to the time
integration schemes.

Table 1: Results from numerical simulations of the flow around a cylinder for
Reynolds 100 and 1000. An em dash symbol indicates that value is not found in
the reference.

Source | Re = 1000 | Re = 100
---|---|---
$S_{t}$ | Amp($C_{L}$) | Mean($C_{D}$) | $S_{t}$ | Amp($C_{L}$) | Mean($C_{D}$)
This work | 0.241 | 1.427∗ | 1.512∗ | 0.174 | 0.331∗ | 1.455∗
Dettmer-Perić ($Sr_{9}$) 24 | 0.252 | 1.498∗ | 1.625∗ | 0.181 | 0.358∗ | 1.495∗
Dettmer-Perić ($Sr_{19}$) 24 | 0.243 | — | — | 0.171 | — | —
Lovrić et.al.32 | — | — | — | 0.165 | 0.292∗ | 0.998∗
Idelsohn et.al.57 | 0.2415 | 1.330 | 1.276 | — | — | —
Idelsohn et.al.58 | 0.2475 | 1.630 | 1.639 | — | — | —

Obtained from figure

##### Velocity-based and Displacement-based formulations

So far all results for the flow around a cylinder have been obtained using the
velocity-pressure formulation. Now these are compared with results using the
displacement-pressure formulation at Re = 100. Results for both formulations
are shown in Fig. 18(b), using the Backward Euler, Newmark (Fig. 18(a)) and
Generalized-$\alpha$ (Fig. 18(b)) time integration schemes.

Subtle differences can be seen in Figs. 18(a) and 18(b). The largest
difference between formulations in the Newmark method is $0.09\%$, obtained
with $\Delta t=0.016$ and $\gamma=0.6$. In the Generalized-$\alpha$, the
largest difference is likewise $0.09\%$ and is found for $\Delta t=0.032$ and
$\rho_{\infty}$ = 0.5. Typically, the larger the numerical damping for high
frequencies, the smaller the difference between velocity- and displacement-
based formulations. We believe, that this subtle difference between the two
formulations is mainly caused by high-frequency perturbations, which stand out
in problems with long simulation times and several remeshing, such as the flow
around a cylinder.

((a))

((b)) Figure 18: Strouhal number of the flow around a cylinder at Re = 100.
Results for velocity-based and displacement-based formulations using (a)
Backward Euler and Newmark’s Method (NM) and (b) Generalized-$\alpha$ method.

(a) (b)

##### Implementation approaches of the Generalized-$\alpha$ method

The last analysis using the flow around a cylinder considers the different
implementation approaches of the Generalized-$\alpha$ scheme. The comparison
is performed for Re = 100 and using five time steps $\Delta t$ ranging from
0.004 to 0.064. First, schemes that write the pressure at $t_{n+\alpha_{f}}$
are compared, i.e.,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, and $\text{GA-
II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$. Two values for the spectral radius
are considered, 0.0 and 0.5. Results for the Strouhal number are summarized in
Fig. 19(a). This shows a small difference between implementation approaches,
which is about $0.1\%$ as indicated in Fig. 19(a). This is also seen in the
amplitude of the lift coefficient and the mean of the drag coefficient,
although no plots are given here to avoid overextending the manuscript.
Thereby, as expected, the hypotheses of scaling state variables (GA-I) or
forces of momentum balance (GA-II) lead to practically the same result,
especially when small time steps are used.

((a))

((b)) Figure 19: Strouhal number of the flow around a cylinder at Re = 100.
Results using different implementation approaches of the Generalized-$\alpha$
scheme. In (a) are implementation approaches that follow the principle of the
Generalized-$\alpha$ method and impose the pressure at $t_{n+\alpha_{f}}$. In
(b) are implementation approaches that write the pressure at $t_{n+1}$ and
integrate the discrete gradient operator ($\mathbf{D}^{\intercal}$) in
$\Omega(t_{n+1})$.

(a) (b)

A comparison is now made between schemes that write the pressure at $t_{n+1}$
in the momentum equation, i.e.,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$. It is worth recalling that two scenarios
are presented, one that integrates the gradient matrix in $\Omega(t_{n+1})$
(Eq. 49) and another that does so in $\Omega(t_{n+\alpha_{f}})$ (Eq. 50).
Results for the first case are presented in Fig. 19(b). The figure also
includes results with
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ (solid line)
that are used as reference. It is apparent from Fig. 19(b) that results are
similar for all implementation approaches when $\rho_{\infty}$ = 0.0. However,
difference between the reference
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ and schemes
that write the pressure at $t_{n+1}$ increases with the value of
$\rho_{\infty}$, in detriment of the schemes that do not follow the principle
of the Generalized-$\alpha$ method. Similar observations are obtained when
analyzing the amplitude of the lift coefficient or the mean value of the drag
coefficient (results not reported here). This behavior was observed for the
kinematics of the solitary wave propagation problem and it is now found in the
pressure field of the flow around a cylinder.

((a))

((b))

((c))

((d))

((e))

((f))

((g))

((h))

((i)) Figure 20: Strouhal number ($S_{t}$) from the flow around a cylinder at
Re = 100. Results obtained with different implementation approaches of the
Generalized-$\alpha$ method. Continuous and dashed lines represent approaches
that write pressure at $t_{n+\alpha_{f}}$ and $t_{n+1}$, respectively.
Different spectral radii are considered: (a-c) $\rho_{\infty}=0.0$, (d-f)
$\rho_{\infty}=0.2$, and (g-i) $\rho_{\infty}=0.5$.

(a) (b) (c)
(d) (e) (f)
(g) (h) (i)

((a))

((b))

((c))

((d))

((e))

((f))

((g))

((h))

((i)) Figure 21: Mean value of drag coefficient ($C_{D}$) from the flow around
a cylinder at Re = 100. Results obtained with different implementation
approaches of the Generalized-$\alpha$ method. Continuous and dashed lines
represent approaches that write pressure at $t_{n+\alpha_{f}}$ and $t_{n+1}$,
respectively. Different spectral radii are considered: (a-c)
$\rho_{\infty}=0.0$, (d-f) $\rho_{\infty}=0.2$, and (g-i) $\rho_{\infty}=0.5$.

(a) (b) (c)
(d) (e) (f)
(g) (h) (i)

The second scenario for
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$ integrates the gradient matrix in
$\Omega(t_{n+\alpha_{f}})$, as done in Eulerian-based formulations that adopt
the modified approach for the pressure. Three values of $\rho_{\infty}$ are
considered, 0.0, 0.2 and 0.5. Results for the Strouhal number are shown in
Fig. 20(i). This summarizes all the implementation approaches and plots on
each graph the schemes that write the pressure at $t_{n+\alpha}$ (solid line)
and $t_{n+1}$ ( dashed line). Unlike the previous scenario (using
$\mathbf{D}_{n+1}$) here the Strouhal number is not significantly influenced.
The difference with respect to
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ is about 0.1
$\%$ for $\rho_{\infty}$ = 0.5. Something similar is observed for the
amplitude of the lift coefficient (not reported here). However, a different
behavior is observed in the drag force. Fig. 21(i) shows the mean of the drag
coefficient. As expected, for $\rho_{\infty}$ = 0.0 all 6 approaches
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, $\text{GA-
II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$,
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$) yield very similar results. However, as
$\rho_{\infty}$ increases, schemes start to differ. The trend is that the mean
drag coefficient grows in the schemes that write the pressure at $t_{n+1}$. In
addition, one observes a notorious difference between implementation
approaches $\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{1,0,0}{p_{n+1}}}$,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$, and $\text{GA-
II}_{\color[rgb]{1,0,0}{p_{n+1}}}$, which has not been exhibited in the
previous example. Interestingly,
$\text{GA-I}^{q_{n+1}}_{\color[rgb]{1,0,0}{p_{n+1}}}$ features a high
dependence on the value of $\rho_{\infty}$, even more, the converged value
($\Delta t\to 0$) seems to be dependent on $\rho_{\infty}$. Although this
behaviour is observed only in the drag force, it is an inconsistency of the
Generalized-$\alpha$ method that arises from writing the pressure at $t_{n+1}$
instead of $t_{n+\alpha_{f}}$, and integrating the gradient matrix in
$\Omega(t_{n+\alpha_{f}})$.

###  4.4 Collapse of a cylindrical water column

The final problem consists of simulating the collapse of a water column whose
initial condition follows the experimental setting of Martin and Moyce 59. The
purpose is to validate performance of the analyzed time integration schemes on
a problem with a different simulation setting. For this reason, the problem is
modelled in 3D, as shown in Fig. 22(a). Slip boundary condition is used in the
base (plane $x$-$y$) and in the symmetry plane ($y$-$z$). The initial water
column radius is $r_{0}=0.05715$ m and its height is $h_{0}=2\>r_{0}$. The
acceleration of gravity is $g=9.81$ m/s2, the density $\rho=1000$ kg/m3 and
the viscosity $\mu=0.001$ Pa s. A non-uniform discretization is chosen for
refining elements on the free surface, as shown in Fig. 22(b). The
characteristic element size is 5 mm at the free surface and 10 mm at a
distance $r_{0}$ from the free surface, while a linear interpolation defines
the size of elements in between. For some problems, the time step size
($\Delta t$) is set as a function of a maximum CFL number set at the beginning
of the simulation. The CFL number of an element is computed as:

$\text{cfl}=\frac{\mathrm{max}(\>\lVert\mathbf{v}\rVert\>)\>\Delta
t}{\mathrm{d}_{\mathrm{eq}}}$ (66)

where $\mathrm{max}(\>\lVert\mathbf{v}\rVert\>)$ is the maximum velocity norm
of the element nodes and $\mathrm{d}_{\mathrm{eq}}$ is the equivalent
spherical diameter of the tetrahedral element. The maximum CFL number is
denoted as CFL∗ and is computed as CFL∗ = max(cfl), max(cfl) being the vector
that gathers the cfl of all elements.

((a))

((b)) Figure 22: Collapse of a cylindrical water column problem. (a) Geometry
and (b) initial finite element discretization.

(a) (b)

Six collapse of a water column problems are solved. One using a fixed time
step $\Delta t=0.00005$ and Backward Euler. The other five consider a time
step set by a maximum CFL∗ = 0.5. These five problems use Backward Euler,
Trapezoidal, and Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$) with
$\rho_{\infty}$ = 0.0, 0.25 and 0.50. From each problem, height ($h$) and base
radius ($r_{b}$) of the water column versus dimensionless time ($t^{*}$) are
reported. The latter is defined as $t^{*}=t\>(2g/r_{0})^{0.5}$. The base
radius is computed using particles located in the basal semicircle and taking
the average distance from these particles to the origin [0 0 0]. The height is
measured at the center of the water column (along the $z$-axis). Results are
summarized in Fig. 23(b) along with references using 3D Eulerian-based 60, 61,
62, 3D Lagrangian FEM 63, and 2D Lagrangian PFEM 48 formulations. Legend of
Fig. 23(a) is shown in Fig. 23(b).

((a))

((b)) Figure 23: Collapse of a water cylindrical column problem. Evolution of
(a) base radius $r_{b}$ and (b) height $h$ of the water column. The legend of
curves in (a) is found in subfigure (b).

(a) (b)

Fig. 23(a) shows that the experimental reference (Martin and Moyce 59)
exhibits a similar pattern to that of the present work, but with the
experimental points shifted to the right. The fluid front delay of Martin and
Moyce has also been noted in experimental studies (e.g., see Fig. 12 in
reference 64), which suggest that the discrepancy with Martin and Moyce may be
due to their experimental technique used for dam gate removal. In particular,
this work is closer to the numerical results of Battaglia et.al.62 and Akin
et.al.60 for both fluid front position ($r_{b}$) and height column ($h$).  The
small differences with respect to those references may be due to the fact that
they use a formulation based on a fixed regular mesh and larger time steps
than the present work (approximately one order of magnitude). Due to the
chosen space discretization, we cannot use larger time steps since tetrahedral
elements could be skewed after updating nodal position. Enlarging element size
in favor of larger time steps would impair the geometric resolution of the
fluid front, since the Alpha Shape algorithm eliminates those highly distorted
elements in the fluid front. In other words, the element and time step sizes
were chosen by balancing the trade-off between computation time and
resolution.

As for the time integration schemes, no significant differences are observed.
At least from a kinematic point of view, results are insensitive to the time
integration scheme, presumably because the problem is solved with time steps
small enough to avoid numerical damping. Looking in more detail, left and
right close-up views in Fig. 23(a) show, respectively, that Backward Euler
with fixed and small time step (green curve), and Trapezoidal (red curve) lead
the front position ($r_{b}$). This alternation between curves leading the
front position is seen everywhere in the graph, although Backward Euler with
large time step (CFL∗ = 0.5, black curve) always exhibits the major delay,
which is a sign of numerical damping.

As it was not possible to find references presenting pressure results for this
example, it is decided to report the pressure at the origin point [x y z] = [0
0 0]. Results for each time integration scheme are shown in Figs. 24(a) and
24(b). These represent the same plot but from different perspectives. Results
show that all time integration schemes deliver the same order of magnitude for
the pressure, but different quantities of spurious oscillations. A simple look
at Fig. 24(a) allows to classify results by the amount of spurious
oscillations. The trapezoidal scheme (CFL∗ = 0.5) has the largest number of
oscillations, which grow significantly from $t^{*}$ = 2. In second place is
the generalized alpha (CFL∗ = 0.5) with $\rho_{\infty}$ = 0.5. Third is
Backward Euler with small time step $\Delta t=0.00005$. Note that this $\Delta
t$ is much smaller than the one imposed by CFL∗ = 0.5, as seen in Figs. 24(c)
and 24(d). In fourth place is generalized alpha (CFL∗ = 0.5) with
$\rho_{\infty}$ = 0.25. Finally, schemes with the smallest amount of spurious
oscillation are Backward Euler (CFL∗ = 0.5) and Generalized-$\alpha$ (CFL∗ =
0.5) with $\rho_{\infty}$ = 0.0.

((a))

((b))

((c))

((d)) Figure 24: Collapse of a cylindrical water column problem. (a) Pressure
at [x y z] = [0 0 0], (b) side view of subfigure (a), (c) maximum CFL number
and (d) time step $\Delta t$ during simulation.

(a) (b)
(c) (d)

Importantly, the choice of the best time integration scheme depends on a
number of criteria that must be weighed by the user. For example, if high
accuracy is required for fluid kinematics (position, velocity, acceleration),
then the trapezoidal or Generalized-$\alpha$ scheme with $\rho_{\infty}\approx
1$ is suggested, as long as convergence of the nonlinear algorithm is
achieved. On the other hand, if precision and stability in the pressure field
is desired, for example for fluid-structure simulation, then schemes with more
numerical damping are recommended such as Backward Euler and
Generalized-$\alpha$ with $\rho_{\infty}\approx 0$. If time step size is
relevant to the computational time, then Generalized-$\alpha$ with
$\rho_{\infty}\approx 0$ is a good choice, since it does not present excessive
numerical damping at large time steps in comparison to Backward Euler. In
particular for the authors of this work, the recommendation is
Generalized-$\alpha$
($\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$) with
$\rho_{\infty}\leq 0.5$.

The last remark concerns the computation time. Notably, the
Generalized-$\alpha$ method allows using larger time steps than Backward Euler
or Newmark and still achieve the same order of error. The ratio of time step
sizes between time integration schemes depends on each problem. For instance,
the reader can verify from the flow around a cylinder problem that the
Generalized-$\alpha$ reaches the same accuracy than Backward Euler, even with
a $\Delta t$ that is 32 times larger than that used with Backward Euler.
Therefore, a reduction in computation time is expected with
Generalized-$\alpha$ due to the use of larger time steps.

To provide the reader with an approximate computation cost of our
Generalized-$\alpha$ implementation, the computation time of three collapse of
cylindrical water column problems are reported below. The problems consider
the Backward Euler (BE),
$\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$, and
$\text{GA-II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$ as time integration schemes.
All three problems use the same time step $\Delta t=0.00005s$.

The average computation times are compared in Fig. 25(a). These are computed
as the total computation time divided by the total amount of iterations of the
nonlinear algorithm (time / $n$). This average computation time is normalized
by that obtained with Backward Euler (timeBE / $n_{BE}$). Thus, Fig. 25(a)
shows that the three problems take the same computation time per iteration of
the nonlinear solver. However, the nonlinear solver requires more iterations
in the GA schemes (about 7$\%$ more iterations), which is reflected in the
total computation time shown in Fig. 25(b). The increased computation time in
problems using GA-I and GA-II is due to the Finite Element Analysis (FEA), as
depicted in Fig. 25(c). On the contrary, the time spent in remeshing and other
tasks (e.g. data input/output) is the same as that required by Backward Euler.
This is reasonable since all three problems use the same time step size, so
the amount of remeshing (and other tasks per time step) are the same.

Fig. 25(d) shows the time spent in different tasks of the FEA code. There, the
time to build the system of equations (matrix $\mathbf{A}$ and vector
$\mathbf{b}$) is slightly longer in problems using Generalized-$\alpha$. This
is because GA-I and GA-II present more terms in the governing equations, and
because there are about 7$\%$ more system builds in problems using
Generalized-$\alpha$ schemes. However, most of the computation time is spent
in solving the system $\mathbf{A}\mathbf{q}=\mathbf{b}$, which is a well known
drawback of monolithic formulations.

In summary, for the same set of simulation parameters, there is no significant
difference in computation time and memory usage between the Backward Euler and
Generalized-$\alpha$, because the bottleneck is the routine that solves the
system $\mathbf{A}\mathbf{q}=\mathbf{b}$. This observation also holds for the
2D problems reported in this work.

((a))

((b))

((c)) ((d)) Figure 25: Breakdown of the computation time in the collapse of a
water column problem. Three time integration schemes are considered: Backward
Euler (BE), $\text{GA-I}^{q_{n+\alpha}}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$
and $\text{GA-II}_{\color[rgb]{0,0,1}{p_{n+\alpha}}}$. Graphs are scaled with
respect to the total time of Backward Euler, denoted as
$\mathrm{time}_{\mathrm{BE}}$. In (a), $n$ denotes the total amount of
iterations carried out by the nonlinear solver. In (d), the task Build
represents the computation of triplets ($i$, $j$, Ai,j), doublets ($j$, bj),
and matrix/vector assembly; Solve involves a pattern analysis of matrix
$\mathbf{A}$, an LU factorization of $\mathbf{A}$, and solving a linear
system; and Others consists of updating nodal state variables according to the
time integration scheme, and verification of convergence criteria.

(a) Time per iteration (b) Total time (c) Time in PFEM code (d) Time in FEA
code

##  5 Conclusions

This work incorporates the Generalized-$\alpha$ time integration scheme into
the Particle Finite Element Method (PFEM) to solve the incompressible Navier-
Stokes equations in a Lagrangian framework. The discretized system of
equations is derived for both velocity-pressure and displacement-pressure
formulations. The time integration scheme is compared with the Backward Euler
and Newmark, which are well established in the PFEM literature. Numerical
comparisons of several benchmark problems and a literature review of the
Generalized-$\alpha$ method for computational fluid mechanics lead to the
following conclusions:

  * •

The Generalized-$\alpha$ method in PFEM outperforms the classical Backward
Euler and Newmark schemes, since it exhibits less numerical damping at large