1 Legato group, Faculty of Science, Technology and Medicine, University of
Luxembourg, Luxembourg 2Fraunhofer ITWM, Germany.

E-mail<EMAIL_ADDRESS><EMAIL_ADDRESS>
# Volume and Mass Conservation in Lagrangian Meshfree Methods

Pratik Suchde1, Christian Leithäuser2, Jörg Kuhnert2, Stéphane P.A. Bordas1

###### Abstract

Meshfree Lagrangian frameworks for free surface flow simulations do not
conserve fluid volume. Meshfree particle methods like SPH are not mimetic, in
the sense that discrete mass conservation does not imply discrete volume
conservation. On the other hand, meshfree collocation methods typically do not
use any notion of mass. As a result, they are neither mass conservative nor
volume conservative at the discrete level. In this paper, we give an overview
of various sources of conservation errors across different meshfree methods.
The present work focuses on one specific issue: unreliable volume and mass
definitions. We introduce the concept of representative masses and densities,
which are essential for accurate post-processing especially in meshfree
collocation methods. Using these, we introduce an artificial compression or
expansion in the fluid to rectify errors in volume conservation. Numerical
experiments show that the introduced frameworks significantly improve volume
conservation behaviour, even for complex industrial test cases such as
automotive water crossing.

Conservation; Meshfree; Lagrangian framework; CFD; Collocation; SPH

##  1 Introduction

Over the last few decades, meshfree or meshless methods have become a popular
alternative to conventional mesh-based solution procedures. Their biggest
advantage lie in avoiding the generation of a mesh to discretize the
computational domain. Meshfree domain discretization procedures that produce a
point cloud or a particle cloud are faster, and can be automated to a higher
degree, even for complex domains [57]. Another major advantage of meshfree
methods, especially when used in a Lagrangian framework, is the ease of
capturing large deformations and displacements in the computational domain. In
this context, the meshfree equivalent of remeshing is much cheaper and faster
[59]. These advantages of meshfree methods, naturally, come with their own
challenges. One of the biggest drawbacks of meshfree methods is a lack of
conservation.

This lack of conservation can appear in two forms. The first is the lack of
conservative properties inherent in derivative approximation procedures. This
has been a topic of significant study in recent years [31, 60, 65], though
several open questions still remain. The second issue of conservation, which
has been largely ignored in meshfree literature, concerns the absence of mass
and volume conservation. This becomes especially relevant in flow problems
with free surfaces which are captured using Lagrangian frameworks. While the
latter point is the topic of investigation of the present work, we explain the
details of both notions of conservation in this paper.

For the purpose of the present work, we consider meshfree methods for fluid
flow divided into two categories. We refer the reader to review papers [5, 39,
8] for an overview and classification of meshfree methods for structural
simulations, which are not considered here. The first category of meshfree
methods for fluid flow is particle-based meshfree methods such as Smoothed
Particle Hydrodynamics (SPH) [19], where each node of the domain
discretization is a mass particle. The second category is meshfree collocation
methods (for example, [24, 23]), where a node is only an approximation
location, and does not inherently carry mass. The former mass-based nodes have
the advantage of a simple notion of mass conservation. Since the mass of a
particle is constant in time, and particles are neither added nor deleted
during a simulation (barring inflows and outflows), the total mass contained
in the system remains constant. As a result, these methods are typically said
to be mass conservative. From a physical or continuum perspective, for
incompressible flow, mass conservation and volume conservation are equivalent.
However, this is not the case for SPH-type meshfree methods, where the fluid
volume is not conserved, not even for incompressible flows with a constant
density. In incompressible free-surface flow simulations, the geometric volume
occupied by the computational domain often decreases as the simulation
progresses, sometimes significantly so. An example of this is shown for a
sloshing simulation in Figure 1. On the other hand, for approximation node
based meshfree collocation methods, since nodes do not carry mass, mass can
only be defined as an auxiliary property. As a result, there is no notion of
mass or volume conservation in these methods.

Figure 1: Lack of volume conservation: Sloshing of an incompressible fluid
with a particle method. The gravity vector is rotated to simulate sloshing,
after which the fluid is brought back to rest. Counter-clockwise from top
left. The top row shows the initial state (left) and the final state (right)
of the fluid, while the bottom row shows intermediate states. The red dashed
line indicates the level of the fluid at the initial state. The figures shows
that the volume of the final state is lower than that present initially. Since
the number of particles are constant, the total numerical mass is conserved.
However, the fluid volume is not conserved.

We note that this notion of volume conservation defects is also an issue for
mesh-based methods for capturing free surface flow. Volume errors are observed
in many mesh-based volume of fluid (VOF) approaches [1], as well as hybrid
Eulerian mesh/Lagrangian particle methods [29].

In mesh-based methods, the definition of an element or control volume
prescribes a discrete volume packet. In contrast, for meshfree methods, there
is no obvious notion of a discrete volume attached to a node. As a result,
meshfree literature uses various different and unequivalent ways to define the
volume of a point or particle. This compounds the issue of volume conservation
errors, and also introduces errors in post-processing simulation results. In
meshfree collocation type methods, a further issue in accurate post-processing
is the lack of native mass definitions, which makes it cumbersome to track,
for example, the transportation and accumulation of mass. A detailed example
of this is presented later in Section 7.5.

The main contributions of this paper are four-fold: (1) We present a first of
its kind survey of different conservation issues in meshfree methods for fluid
flow. (2) We give an overview of different notions of defining discrete
volumes for meshfree methods, and present their advantages and disadvantages.
(3) We introduce the novel notion of representative masses and densities for
meshfree methods, which simplify post-processing in collocation type methods.
(4) Using the introduced representative masses and densities, we introduce an
artificial point displacement to significantly reduce volume conservation
errors. This is done based on the difference between the representative
density and the actual (physical) density of the fluid. We note that the
notion of artificial displacement has been widely adopted in the SPH community
(for example, [32, 33, 45]). SPH literature uses artificial displacement to
avoid particle clustering, which leads to numerical instabilities. In
contrast, in the present work, we use a similar notion, but towards the goal
of fluid volume conservation. While this work focuses on the mass and volume
conservation in the context of collocation type meshfree methods, many of the
ideas are also directly applicable to mass particle based meshfree methods.

The paper is organized as follows. In Section 2 we introduce some further
preliminaries about different types of meshfree methods, the notation used,
and the governing equations of flow used in the present work. We then give a
brief overview of different conservation issues across all meshfree methods in
Section 3. In Section 4, we present different methods for defining discrete
volumes in meshfree methods. Section 5 then introduces the notion of
representative masses, followed by the notion of representative densities in
Section 6, which also introduces the volume correction algorithm. A series of
numerical test cases are presented in Section 7 followed by a summary and
conclusion in Section 8.

##  2 Preliminaries

In this section we introduce some preliminaries regarding the notation and
governing equations used in the present work, along with a description of
different types of meshfree methods.

###  2.1 Types of meshfree methods

As explained above, in the context of the current work, we divide Lagrangian
meshfree methods for fluid flow into two broad classes, each with their own
advantages and disadvantages.

  1. 1.

Mass-particle based meshfree approaches. Here, each node is a mass particle.
These types of methods include the popular Smoothed Particle Hydrodynamics
(SPH) [19] and derived methods, such as Moving Particle Semi-explicit (MPS)
[28], Moving Particle Explicit (MPE [63] or MPS [68]) and Discrete Droplet
Method (DDM) [7], and also some hybrid mesh-meshfree methods like Material
Point Methods (MPM) [40].

  2. 2.

Approximation-point based meshfree methods. Here, a node does not carry mass,
and is only a collocation point where the governing PDEs are solved. While
Radial Basis Function (RBF) based collocation methods have been widely used
(for example, [11, 46, 52]), Generalized Finite Difference Methods (GFDM) are
more commonly used in Lagrangian approaches for fluid flow (for example, [13,
56]).

The main advantage of meshfree particle methods is that particles having mass
directly implies a notion of mass conservation, which is not present in
meshfree collocation approaches. On the other hand, the primary advantage of
approximation point based meshfree methods is the ease of adding and
merging/deleting points. This is very helpful to fix distortion as a result of
movement. Points can easily be deleted or merged when they come too close, and
points can be added when artificial “holes” are created in the point cloud.
This advantages also means that adaptive refinement is very straight forward
in these methods. While SPH type methods are developing multi-resolution
frameworks [27, 70] with two or three discrete levels of refinement in the
domain, collocation type meshfree methods have for long had the ability to
achieve continuous resolution adaptivity [6, 13, 18, 59].

Without the ability to freely add and delete particles, particle distortion
can not be fixed as easily in SPH type methods. Two common approaches are used
for this: (i) redistributing particles based on a background mesh [41, 42,
43], or the more widely adopted (ii) artificially moving the particle cloud.
For this, a non-physical or artificial movement of particles is performed to
avoid particle clustering [32, 33, 45]. These approaches have been referred to
as particle shifting. In the present work, we introduce the use of similar
ideas to provide improved volume conservative behaviour.

######  Remark.

The topic of volume conservation has been widely discussed in recent SPH
literature [25, 36, 44, 62]. However, those discussions are in the context of
particle shifting. Classical particle shifting methods introduced non-physical
volume, while the volume conservation work in SPH aims to conserve volume
during the particle shifting process. In contrast, the present work deals with
volume conservation defects from other sources of the numerical scheme, and
aims to perform an artificial point shifting to correct that. Another approach
for volume conservation considered in the context of meshfree collocation
methods is the so-called Position-Based Dynamics [2, 3, 4], which aims to keep
the volume of each point constant as points move in a Lagrangian fashion. In
contrast, here we aim to keep the global fluid volume constant.

Another significant distinction between these two types of meshfree methods is
the ease of enforcing boundary conditions. While a wide range of boundary
conditions can be easily incorporated in meshfree collocation type methods, it
is not straight forward to do the same in particle based meshfree methods
[34].

###  2.2 Notation

The bulk of this work uses an approximation-point based collocation meshfree
method. We consider a computational domain $\Omega=\Omega(t)$, with boundary
$\partial\Omega=\partial\Omega(t)$. The domain is discretized with a point
cloud composed of $N=N(t)$ points, consisting of points both in the interior
and on the boundary of the domain.

For a point $i=1,2,...N$, all approximations are done on a set of nearby nodes
referred to as its neighbourhood or support $S_{i}$ consisting of $n_{i}$
nearby points, within a distance of $h=h(\vec{x},t)$. $h$ denotes the so-
called smoothing length or interaction radius, which governs not just the
support size, but also the resolution of the point cloud. This is because
inter-point distances are controlled to fixed multiples of $h$. The minimum
distance between two points is controlled to be greater than $r_{min}h$, and
maximum hole size in the point cloud is given by $r_{max}h$. This is enforced
by merging particles that come too close, and adding them in holes where no
points are present. For more details of these addition / deletion procedures,
we refer to [15, 51, 57, 59].

The point clouds used in the present work are irregularly spaced, without any
background grid or mesh. All points clouds are generated using an advancing
front point generation method [57, 54].

###  2.3 Governing equations / model problem

The governing equations considered here are a unified framework for fluids and
solids, consisting of the conservation of mass and momentum equations, used in
a Lagrangian framework. D →xDt &= →v ,
D ρDt = - ρ∇⋅→v ,
D →vDt = 1ρ ∇⋅S \- 1ρ ∇p + →g , for position $\vec{x}$, velocity $\vec{v}$,
pressure $p$, density $\rho$, gravity and other external forces $\vec{g}$, and
a stress tensor $\mathbf{S}$ consisting of both viscous and solid terms.

$\mathbf{S}=\mathbf{S}_{\text{visc}}+\mathbf{S}_{\text{solid}}\,,$ (1)

with

$\mathbf{S}_{\text{visc}}=\eta\left(\left(\nabla\vec{v}\right)+\left(\nabla\vec{v}\right)^{T}-\frac{2}{3}\nabla\cdot\vec{v}\,\mathbf{I}\right)\,.$
(2)

The exact formulation of $\mathbf{S}_{\text{solid}}$ or of other material
properties, such as $\rho$ or the viscosity $\eta$ is not very important in
the present context. The ideas introduced in this paper could be used with a
variety of material modelling frameworks. For simplicity, we use
incompressible Newtonian flow with a constant density and viscosity, with
$\mathbf{S}_{\text{solid}}=0$ throughout this paper. If required, an energy
conservation equation could also be used. While the ideas of this work can be
easily applied to any governing equations, here we restrict the use to fluid
flow.

##  3 Conservation and Meshfree Methods

While we focus on one aspect of conservation, it is important to understand
the other reasons of conservation defects to correctly interpret numerical
results. In this section, we give a brief overview of different conservation
related errors in meshfree methods for fluid flow. Each of these apply to both
SPH-type and collocation meshfree methods. We note that these issues only
cover Lagrangian meshfree methods for fluid flow. Eulerian weak form meshfree
methods with stationary background meshes used in structural modelling are not
covered here.

###  3.1 Conservation of fluxes

Consider a conservation law

$\frac{\partial\phi}{\partial t}+\nabla\cdot\mathbf{J}=0\,,$ (3)

with $\mathbf{J}=\mathbf{J}(\phi)$. For sufficiently smooth $\phi$ and
$\mathbf{J}$, integrating Eq. Eq:ConservationLaw over the entire domain, and
an application of the divergence theorem leads to the integral form of the
conservation law,

$\frac{d}{dt}\int_{\Omega}\phi\,dV=-\int_{\partial\Omega}\vec{n}\cdot\mathbf{J}\,dA\,.$
(4)

Due the local nature of meshfree methods and the absence of a global mesh
discretizing the domain, most meshfree methods do not posses a discrete
divergence theorem. As a result, at the discrete level, there is no connection
between the differential law Eq. Eq:ConservationLaw and the integral form Eq.
Eq:EnergyConservation. No matter how accurately the differential form is being
solved at each node, there is no discrete equivalent of the global (integral)
conservation law.

Due to the absence of a mesh, a definition of flux is not natural to meshfree
discretizations. To introduce global conservation in meshfree methods, [12,
31, 65] introduce different algebraic notions of fluxes to obtain a discrete
divergence theorem. However, each of these methods require a global
computation either of the differential operators [31] or of the fluxes [65].
As a result, they are only used as a pre-processing step for meshfree
simulations on fixed point clouds, and can be quite expensive for moving point
clouds. To overcome this, [60] introduced a local definition of fluxes that is
much quicker to compute, but only provides approximate conservation. For
moving Lagrangian meshfree frameworks, achieving global conservation without
expensive global computations remains an open problem.

###  3.2 Volume conservation

As discussed earlier and illustrated in Figure 1, the lack of volume
conservation of the discrete domain is another source of numerical inaccuracy.
This topic is the main focus of the present work.

####  3.2.1 Lagrangian motion

One of the biggest source of volume defects in Lagrangian methods is
inaccurate advection. Lagrangian advection does not introduce diffusion, as is
typically the case in Eulerian advection. It is thus often a preferred choice
for capturing free surface flow. However, inaccurate advection for flows with
free surfaces introduces volume conservation errors.

Most Lagrangian meshfree methods consider a Lagrangian movement step by
updating node positions using a first order accurate method, by treating the
fluid velocity as constant between time steps. This has been shown to be
extremely inaccurate, and higher order Lagrangian movement has been shown to
significantly reduce volume conservation defects [58]. In our present work, we
follow a higher order movement approach introduced in [58], to reduce the
influence of inaccurate Lagrangian motion on volume conservation errors. It is
important to note here that this issue is relevant in all Lagrangian, semi-
Lagrangian and ALE frameworks, covering meshfree, mesh-based, and hybrid
methods. However, this aspect has been widely ignored in literature.

The topic of volume conservation as a result of Lagrangian motion is usually
mentioned in the context of the distortion of point or particle clouds. As
described in the remark of Section 2.1, volume conservation in SPH generally
refers to preventing volume changes in the artificial particle shifting step,
and not overall. While in meshfree collocation methods, local volume
conservation is enforced to conserve individual particle volumes [4].

####  3.2.2 Mass and volume definitions

An important factor in the lack of volume conservation in Lagrangian meshfree
methods is that a volume definition is not native to meshfree methods. This
topic will be discussed in more detail in Section 4.

####  3.2.3 Other modelling and numerical inaccuracies

Volume conservation errors could also arise due to various other numerical
inaccuracies. These could be due to the design of the complete numerical
scheme (for instance, using a weakly compressible scheme to approximate
incompressible equations), the tolerance of the iterative solvers, the
(derivative) approximation error, or some other numerical source of error.

Our present work combines issues relating too all these sources of volume
errors. We introduce a volume correction scheme to address volume inaccuracies
that could have arisen due to any of the above mentioned reasons.

##  4 Volume definition

Unlike the case of meshes, volume definitions are not native to meshfree
methods. Before taking a look at the different approaches used in literature,
we first define the concept of a consistent volume definition for a meshfree
node cloud. A node here could be either a mass-carrying particle or an
approximation point. For a volume definition to be consistent, the sum of
volumes of all meshfree nodes in a domain should represent the physical volume
occupied by the computational domain. For example, if a cube of side $l$ is
discretized by meshfree nodes, the sum of volumes of all nodes should
approximate the volume of the cube, $l^{3}$. Such a consistent volume
definition is needed for accurate post-processing of simulation results.

######  Remark.

In the Lagrangian frameworks considered here, the node positions change in
each time step. Thus, the volume of each node must be recomputed at every time
step.

(a) (a) Spherical volume

(b) (b) Global Voronoi decomposition

(c) (c) Local Voronoi decomposition (with partial overlapping volumes)

Figure 2: Different volume definition of point cloud based meshfree methods.

In meshfree literature, several different approaches have been considered to
prescribe a notion of a volume to each meshfree point or particle. Both
geometric and algebraic notions of volume have been proposed. A few commonly
used examples are considered below:

  * •

The simplest approach is to associate a node with the volume of a sphere of a
specified diameter. This could be based on the neighbourhood size or kernel
width, or any other prescribed diameter. This method is not consistent due to
either large overlapping volumes, or unaccounted regions in the middle of the
domain, as is evident from Figure 2.

  * •

Another way to prescribe volumes to each point or particle is to construct a
tessellation on the point cloud. Occasionally, for example [22], meshfree
literature has used the notion of a globally constructed Voronoi tessellation
of a point cloud. While this method is consistent, it can prove to be very
expensive, especially when it needs to be reconstructed at every time step.

  * •

Another possible way to define volumes is through local tessellations. For
each point $i$, its volume $V_{i}$ is given by the volume of the Voronoi cell
locally defined on the support $S_{i}$. The locally defined Voronoi cells need
not stitch together to form a global tessellation of the domain, as can be
seen in Figure 2. However, the amount of overlap is quite small [60] and thus
this method is consistent. It is also much faster than computing a global mesh
on the set of points at every time step. We note here that actual local
Voronoi tessellation need not be computed, as the shape of the elements and
their faces are not needed. Figure 2 only shows the shapes for visualization.
In reality [55], only the nearest neighbours in a Delaunay sense need to be
identified locally, with the resultant simplex volumes giving the volume for
the point in question.

  * •

In most mass-particle meshfree methods, each particle carries a mass and a
density. This results in the obvious notion of volume as
$V_{i}=\frac{m_{i}}{\rho_{i}}$, for mass $m$ and density $\rho$ [34, 53].
During the initial particle generation process, masses are distributed to
particles in a manner so as to make this method consistent. However, this
method does not remain consistent once particles start moving. Other notions
of volume include different algebraic manipulations of the mass and density
[17, 48].

  * •

Several other algebraic notions of volume have been proposed. These include
virtual notions of volume to enforce a discrete divergence theorem, as
explained in 3.1. These methods typically solve global constrained
optimisation problem to enforce an algebraic constraint, while maintaining a
consistent volume definition in the whole domain [31]. However, due to the
optimisation process, the volumes are no longer consistent on sub-domains.

  * •

Another algebraic approach is to use partition of unity functions. A
differential volume at any arbitrary location is partitioned to the nearest
nodes using a weighted kernel function with a normalization to ensure a
partition of unity. This results in a Voronoi tessellation with smoothed edges
[21].

######  Remark.

Analogous to the definition of volume of each meshfree node is the question of
prescribing an area for each boundary node. Similar approaches to the volume
definition can be used here, though this topic has not been widely discussed.

An important special case of prescribing volumes is for nodes with very few
neighbours, which often occurs in violent free surface flows. Section 7.5 is
an example of such a flow pattern. When a node does not have sufficient
neighbours, several of the consistent methods for defining a volume element
explained above can not be used. In the present work, we use local Voronoi
tessellations for defining volumes of each collocation point, as described in
[55, 60]. If the number of neighbours are insufficient to compute a local
Voronoi tessellation, we use a spherical volume such that the volume of the
point closely matches its volume before it became “isolated”.

##  5 Representative Masses

We now introduce the notion of representative mass for points in a meshfree
collocation method. As a stand alone method, this does not affect the
numerical scheme in any way, and is only a very useful post-processing tool.
The need for the representative masses for accurate post-processing is
highlighted in the numerical example in Section 7.5. This framework will also
be the base on which the volume correction algorithm (see Section 6 ) will be
built.

Consider a point $i$ with neighbourhood $S_{i}$, density $\rho_{i}$, and
notional volumes $V_{i}$ computed as explained in Section 4. The
representative mass of point $i$, denoted by $\hat{m}_{i}$, at the first time
step is given simply by the physical density times the volume

$\hat{m}_{i}(t=0)=\rho_{i}V_{i}\,.$ (5)

We note that unlike the case of mass particles in SPH-type methods, these
representative masses can be shifted from one point to another when needed.
This becomes essential when points are added or deleted for adaptive
refinement, or to maintain a quasi-regularity of the point cloud. The
procedures for the addition and deletion of points has been well established
in meshfree collocation literature (for example, [13, 59]), and will not be
covered in the present work. The following subsections introduce a procedure
to modify the representative masses as points are advected, added or deleted.

###  5.1 Mass adaptation

We first establish a general adaptation procedure for redistribution of
representative masses when the point cloud changes. Special cases of this
procedure will be used when points are added or deleted. The same procedure
will also be used for a “smoothing” of representative masses. If this
smoothing procedure is not used, we empirically observe a high mass collection
for a few points near where a lot of point adaption is needed, for example,
near moving wall boundaries.

For a point $i$ with representative mass $\hat{m}_{i}$ computed at the
previous time step, the change in representative mass in the next time step is
denoted by $\Delta\hat{m}_{i}$. We use the convention where a positive sign
$\Delta\hat{m}_{i}>0$ means that point $i$ gains mass. For the mass adaptation
procedure, we first establish a set of conditions or rules that need to be
fulfilled.

  1. 1.

_Total mass is constant_. Barring inflows and outflows (treated separately in
Section 5.5 ), no mass should be created or destroyed during the mass
adaptation process. Thus, we have

$\sum_{i=1}^{N}\Delta\hat{m}_{i}=0\,,$ (6)

where $N=N(t)$ is the total number of points in the domain. Thus, any mass
being added to a point must be removed from another, and vice-versa. For this,
we introduce the notation $\Delta\hat{m}_{j\rightarrow i}$ to denote mass
being transferred from point $j$ to point $i$. Thus, the change of
representative mass of point $i$ can be given by

$\Delta\hat{m}_{i}=\sum_{\substack{j=1\\\ j\neq
i}}^{N}\Delta\hat{m}_{j\rightarrow i}\,.$ (7)

This uses the same convention as before, $\Delta\hat{m}_{j\rightarrow i}>0$
implies point $j$ is losing mass, and point $i$ is gaining mass.

  2. 2.

_Mass adaptation must be local_. This procedure has to be entirely local, in
the sense that a mass redistribution from points $j$ to $i$ should only happen
if $j$ is in the neighbourhood of $i$. Thus,

$\Delta\hat{m}_{j\rightarrow i}=0\,,\text{if}j\notin S_{i}\,.$ (8)

Thus, Eq. Eq:Mass_GlobalSum reduces to

$\Delta\hat{m}_{i}=\sum_{\substack{j\in S_{i}\\\ j\neq
i}}\Delta\hat{m}_{j\rightarrow i}\,.$ (9)

  3. 3.

_Distance dependence of mass adaptation_. Within the support domain of a point
$i$, the mass transfer from a neighbouring point $j$ to $i$ should decrease as
the distance between the points increases.

  4. 4.

_Neighbour mass dependence of mass adaptation_. In the support domain of a
point $i$, for equidistant neighbours $j_{i}$ and $j_{2}$,

$\hat{m}_{j_{1}}>\hat{m}_{j_{2}}\implies\Delta\hat{m}_{j_{1}\rightarrow
i}\geq\Delta\hat{m}_{j_{2}\rightarrow i}\,.$ (10)

More mass should be transferred from a neighbour with more mass, all other
conditions being equal. This prevents mass clumping in certain points.

  5. 5.

_Computations must be local_. The computation of the redistribution should be
a local procedure. Since points move in a Lagrangian framework in each time
step, and addition/deletion of points also has to be done in each time step,
the mass adaptation procedure is run in each time step. Avoiding a global
calculation of the representative mass redistribution can significantly speed
up the redistribution algorithm.

Under theses rules defined above, the updated representative mass of each
point should be modified to satisfy, if possible,

$\hat{m}_{i}+\sum_{\substack{j\in S_{i}\\\ j\neq
i}}\Delta\hat{m}_{j\rightarrow i}=\hat{m}_{i}^{\text{target}}\,,$ (11)

where $\hat{m}_{i}^{\text{target}}$ is the target mass desired. For most
cases, we have

$\hat{m}_{i}^{\text{target}}=\rho_{i}V_{i}\,,$ (12)

which relates the updated representative mass to the physical density and
volume, as done during the initialization step, Eq. Eq:RepMassInitial. Outflow
boundaries are an exception to this target mass prescription, and their
treatment will be described in Section 5.5. For the introduced representative
mass redistribution task, the following observations can be made

  * •

There is no unique solution. To enforce uniqueness, an additional optimality
constraint would need to be added, for example, to minimize the representative
mass changes, $\text{min}\sum_{i=1}^{N}|\Delta\hat{m}_{i}|$.

  * •

The mass of a point $i$ can only change when the mass of at least one of its
neighbours $j$ is also changing. Thus, determining any solution to this
problem, with or without the added optimality constraint, would require a
global procedure to determine the mass changes. However, this violates rule 5.

To prevent the need for global computations, we modify the redistribution
procedure to only approximate Eq. Eq:Mv3. In a first step, for each point $i$,
local mass changes $\Delta\hat{m}_{j\rightarrow i}$ are computed to satisfy
Eq. Eq:Mv3 exactly. These mass changes are only computed and stored in this
step, without determining the resultant new masses. Once the local mass
changes for each point have been determined, they are communicated to
neighbouring points across MPI processes. Finally, the sum of changes is
applied to each point. Thus, a point $i$ will have a change in representative
mass due the redistribution Eq. Eq:Mv3 computed at point $i$ itself, and also
for Eq. Eq:Mv3 computed at each neighbouring point of $i$. For example, the
net change of mass of point $i$ due to its neighbour $j$ is given by

$\Delta\hat{m}_{ij}=\Delta\hat{m}_{j\rightarrow i}-\Delta\hat{m}_{i\rightarrow
j}\,,$ (13)

where $\Delta\hat{m}_{j\rightarrow i}$ is computed in the local computation of
point $i$, and $\Delta\hat{m}_{i\rightarrow j}$ is computed in the local
computation of point $j$. The entire procedure is also explained in the
algorithm in Section 5.6.

For the local redistribution process at each point $i$, we need to prescribe
the form of $\Delta\hat{m}_{j\rightarrow i}$. For this, we introduce the
following ansatz

$\Delta\hat{m}_{j\rightarrow i}=\zeta_{i}K_{ij}\hat{m}_{j}\,,$ (14)

with transfer coefficients $K_{ij}$ and a mass transfer function $\zeta_{i}$.
By definition, rule 4 is satisfied. To satisfy the distance dependence rule 3,
the transfer coefficients are based on the distance between the nearby points

$K_{ij}=\exp(-\alpha_{K}r_{ij}^{2})\,,$ (15)

with a kernel width $\alpha_{K}=2.0$, and where $r_{ij}$ is the relative
distance between points $i$ and $j$

$r_{ij}^{2}=\frac{\|\vec{x}_{j}-\vec{x}_{i}\|^{2}}{\frac{h_{i}^{2}+h_{j}^{2}}{2}}\,.$
(16)

Now, from Eqs. Eq:Mv3, Eq:mTarget and Eq:DeltaM, we get the mass transfer
function

$\zeta_{i}=\frac{\rho_{i}V_{i}-\hat{m}_{i}}{\sum_{\substack{j\in S_{i}\\\
j\neq i}}K_{ij}\hat{m}_{j}}\,.$ (17)

First, $\zeta_{i}$ is computed and stored for each point. After that, the
representative masses are adapted for each point based on its own mass
transfer condition as well as its neighbours’ mass transfer conditions.

###  5.2 Addition of points

As explained above, collocation points are added to the domain to prevent
artificial hole formation due to the Lagrangian movement, and when adaptive
refinement is desired. A representative mass is assigned to a newly added
point $i$ by running the mass adaptation procedure explained in Section 5.1
with $\hat{m}_{i}=0$, indicating that the point had no mass before the
redistribution procedure. The mass transfer function is then reduced to

$\zeta_{i}=\frac{\rho_{i}V_{i}}{\sum_{\substack{j\in S_{i}\\\ j\neq
i}}K_{ij}\hat{m}_{j}}\,.$ (18)

where $S_{i}$ is the new support of the newly added particle. In addition, we
also need to set the transfer coefficient $K_{ij}=0$ if point $j$ is also a
newly added point in that time step.

###  5.3 Deletion / Merging of points

When two points come too close to each other, one of them is deleted, and the
other is moved to the mid-point of the original two points. To adjust the
representative masses, first all the relevant points are flagged for deletion.
After that, the masses are redistributed, followed by the actual deletion of
flagged points. If point $i$ is being deleted, or is being made inactive such
that it no longer takes part in the approximations, we carry out the mass
adaptation procedure laid out in Section 5.1 with $V_{i}=0$, indicating that
the deleted point should have no volume. Thus, the mass transfer function
given in Eq. Eq:Zeta becomes

$\zeta_{i}=\frac{\hat{m}_{i}}{\sum_{\substack{j\in S_{i}\\\ j\neq
i}}K_{ij}\hat{m}_{j}}\,.$ (19)

Additionally, if point $j$ is also flagged for deletion, we set the transfer
coefficient $K_{ij}=0$, indicating that no mass is transferred between the two
points.

###  5.4 Smoothing procedure

To prevent an artificial accumulation of representative masses in a region,
and to ensure that the representative masses are never negative, we run the
general adaptation procedure established in Section 5.1 for each point in the
domain as a sort of smoothing of the masses.

###  5.5 Inflow and outflow boundaries

To maintain a notion of quasi-regularity at inflow and outflow boundaries,
points on inflow and outflow boundaries are kept fixed in an Eulerian
framework. The Lagrangian framework is only used for interior points, free
surface points, slip boundary points, and interface points between two phases.
Inflow points inject new points into the domain along the inward pointing
normal. The new representative mass of these points is determined in the same
as that for other newly added points, see Section 5.2.

For a point $i$ on an outflow boundary or a porous medium boundary, the target
representative mass in Eq. Eq:Mv3 is computed by considering the mass flux
across the boundary

$\frac{d\hat{m}_{i}}{dt}\approx\frac{\hat{m}_{i}^{\text{target}}-\hat{m}_{i}}{\Delta
t}=-\rho_{i}\left(\vec{v}_{i}\cdot\vec{n}_{i}\right)A_{i}\,,$ (20)

where $\vec{n}_{i}$ is the outward facing unit normal, and $A_{i}$ is the area
occupied by the boundary point, computed with local tessellations in a similar
manner to that done for volumes, as explained in Section 4.

###  5.6 Representative mass summary

The procedure for using the introduced representative mass as a post-
processing tool in a meshfree collocation flow solver is summarized in
Algorithm 1, where the italicized steps indicate those not done in a typical
flow solver.

Algorithm 1 Representative masses for post-processing

1:_Initialize representative masses_ (Eq. Eq:RepMassInitial)

2:while Time-stepping loop do

3: Move point cloud in a Lagrangian sense [58]

4: Update neighbour tree

5: Add new points, flag points to be deleted

6: Volume computation (Section 4)

7:  _Update representative masses for addition, deletion_ (Sections 5.2, 5.3)

8: Delete points flagged for deletion

9:  _Smooth representative masses_ (Section 5.1)

10: Flow solver: velocity, pressure solve, stress tensor update if applicable

11: Post-processing calculations

12:end while

##  6 Volume Correction Algorithm

Using the notion of representative mass introduced in Section 5, we now
introduce the volume correction mechanism.

###  6.1 Representative Density

As a preliminary step for volume correction, we introduce a notion of
representative density based on the updated representative masses as

$\hat{\rho}_{i}=\frac{\sum_{j\in S_{i}}\hat{m}_{j}}{\sum_{j\in
S_{i}}V_{j}}\,.$ (21)

Note here that $i\in S_{i}$. We emphasize the notation difference: $\rho$
denotes the physical density, while $\hat{\rho}$ denotes the representative
density.

###  6.2 Volume correction

Now we introduce the novel volume correction mechanism based on the relative
difference between the representative density and the actual physical density.
For each point $i$, a local correction is introduced which resembles an
artificial velocity divergence.

$\left(\nabla\cdot\vec{v}\right)^{\text{artificial}}_{i}=\frac{1}{\Delta
t}\text{min}\left\\{\frac{\hat{\rho}_{i}-\rho_{i}}{\rho_{i}},\gamma\right\\}\,,$
(22)

where $\left(\nabla\cdot\vec{v}\right)^{\text{artificial}}_{i}$ is the
artificial or correction velocity divergence at point $i$, and $\gamma$ is a
user defined parameter signifying the maximum allowed correction per time-
step. Eq. Eq:DivVcorr prescribes and additional compression or expansion rate
of the point cloud to artificially generate or destroy volume.

The artificial velocity divergence introduced in Eq. Eq:DivVcorr can be used
to correct the fluid volume in one of two ways. The first is to include it
within the main numerical scheme. This can be done by adding the term in the
right hand side of the pressure Poisson equation, which is straight forward
irrespective of whether a segregated or coupled approach is being used to
solve for the velocity and pressure fields. The advantage of this approach is
that it is fast, as there no extra global steps. However, the significant
disadvantage is that the actual pressure and velocity are affected.
Empirically, we observe that this approach leads to fluctuations in the
pressure field, and often also leads to unstable simulations.

An alternative approach to use the artificial velocity divergence introduced
in Eq. Eq:DivVcorr, which will be employed in the present work, is to solve
for an artificial displacement independent of the main numerical scheme for
velocity and pressure. For this, we solve a pressure Poisson equation similar
to that obtained in classical projection schemes [9, 10]

$-\frac{1}{\rho}\frac{\partial\rho}{\partial p}\frac{1}{\Delta
t}p^{\text{artificial}}+\nabla\cdot\left(\Delta t\frac{1}{\rho}\nabla
p^{\text{artificial}}\right)=\left(\nabla\cdot\vec{v}\right)^{\text{artificial}}\,,$
(23)

for the artificial pressure $p^{\text{artificial}}$. This is subsequently used
to create the artificial displacement, which can be computed locally at each
point

$\Delta\vec{x}^{\,\text{artificial}}=-\Delta t\left(\frac{\Delta
t}{\rho}\nabla p^{\text{artificial}}\right)\,.$ (24)

Since this artificial pressure is only used to create an artificial point
displacement, the actual velocity and pressure are not affected. As a result,
this approach does not introduce any spurious effects into the numerical
solution.

###  6.3 Volume correction summary

The resultant procedure for a volume corrected meshfree collocation flow
solver is summarized in Algorithm 2, where the italicized steps indicate those
not done in a typical flow solver.

Algorithm 2 The volume correction mechanism

1:_Initialize representative masses_ (Eq. Eq:RepMassInitial)

2:while Time-stepping loop do

3: Move point cloud in a Lagrangian sense [58]

4: Update neighbour tree

5: Add new points, flag points to be deleted

6: Volume computation (Section 4)

7:  _Update representative masses for addition, deletion_ (Sections 5.2, 5.3)

8: Delete points flagged for deletion

9:  _Smooth representative masses_ (Section 5.1)

10:  _Compute Representative densities_ (Eq. Eq:RepDens)

11: Flow solver: velocity, pressure solve, stress tensor update if applicable

12:  _Compute artificial divergence_ (Eq. Eq:DivVcorr)

13:  _Compute artificial pressure_ (Eq. Eq:VolCorr)

14:  _Artificial point movement_ (Eq. Eq:PointDispArt)

15: Post-processing calculations

16:end while

##  7 Numerical Examples

We now present the application of the introduced representative mass and
volume correction frameworks to a series of numerical tests of varying
complexities, starting from simple academic test cases to actual industrial
test cases. The test cases used and the relevance of each test case is
summarized in Table 1.

To quantify the loss in volume, we use an error defined as the relative
deviation of the numerical volume of the fluid being simulation from the
expected value. Thus, the error is given by

$\epsilon_{V}=\frac{|V_{\text{analytical}}-V|}{V_{\text{analytical}}}\,,$ (25)

where $V=\sum_{i=1}^{N}V_{i}$ is the sum of volumes of each point in the
domain, and $V_{\text{analytical}}$ is the analytical volume expectation.

Table 1: Numerical experiments used to highlight the effectiveness of the representative mass and volume correction frameworks. $d$ shows the dimensionality of the problem, with the computational domain $\Omega\subset\mathbb{R}^{d}$. For the three dimensional test cases, the Reynolds number is indicated by $Re$.  |  Test case |  Comments
---|---|---
1 |  Adaptive refinement |  $d=2$
Simplified test for stationary fluid
Illustration of representative mass adaptation
2 |  Dam break |  $d=2$
Common meshfree benchmark test case
3 |  Drops falling and collection |  $d=3$
Convergence with spatial resolution
$Re\approx 15$
4 |  Partially filled screw conveyor |  $d=3$
Industrially motivated test case
Fluid disappears without volume correction
$Re\approx 0.2$
5 |  Automotive water crossing |  $d=3$
Industrial test case
Violent free surface flows
$Re\approx 10^{7}$

### Numerical scheme

The representative mass and volume correction frameworks introduced in the
present work can be used with any choice of time integration scheme, and any
meshfree collocation method for derivative discretization. We briefly
introduce the methods used in the numerical examples for the sake of
completeness.

We use a Generalized Finite Difference Method (GFDM) approach to compute
numerical derivatives. The GFDM is a strong form meshfree collocation approach
that has been widely used [14, 35, 67, 69], and shown to be a robust framework
for a variety of flow applications [38, 64] including flow through porous
media [47], non-Newtonian flow [50, 66], and even for soil mechanics [37], and
elasticity problems [49]. The GFDM approach used here is second order
accurate, with polynomial basis functions.

The overall time integration procedure is explained in Algorithm 2. The basic
scheme starts with Lagrangian movement with a second order approach [58]. This
is followed by a monolithic scheme where the velocity and pressure are solved
together in a single system, using an implicit first-order accurate time
integration scheme. A verification and validation of the basic scheme used
(without volume correction), and the implementation, can be found in our
earlier work [13, 60]. For more details of the scheme used, we refer to [26,
30, 61]. The numerical solvers used and the novel frameworks introduced in the
present work are part of the in-house developed software suite MESHFREE [16].

###  7.1 Adaptive Point Cloud Refinement

For the first test case, we consider a stationary fluid with a free surface in
a square container. In the middle of the domain, the point cloud is refined
(see Figure 3) till the total number of points in the domain approximately
doubles. The point cloud is then coarsened towards a uniform resolution, with
the simulations terminated at a fixed time $t=2$. This example illustrates one
of the main advantages of meshfree collocation methods over meshfree particle
methods: the ease of performing adaptive refinement (see Section 2.1). Since
no external force is applied, the fluid should remain at rest, which is also
observed in the simulations.

The simulation domain is taken as a unit square $[0,1]^{2}$, with the free
surface located at $0.7$ in the height direction. Thus, the analytical volume
occupied by the domain is $0.7$. Figure 4 shows that the total numerical
volume $V=\sum_{i=1}^{N}V_{i}$ remains constant during both refinement and
coarsening, and closely matches the analytical volume. Quantitatively, the
total numerical volume deviates from the analytical one by less than $0.5\%$.
This shows that using the locally defined tessellations for determining
volumes (see Section 4) produces very good and consistent results. It is
important to note here that the presence of sharp boundaries in the present
case, and curved boundaries in a general scenario, will always result in
imperfections in the prescription of discrete volumes.

Figure 4 also shows that the total representative mass
$\hat{\mathbb{M}}=\sum_{i=1}^{N}\hat{m}_{i}$ in the domain remains conserved
during the adaptive refinement process. The fluid density is taken as
$\rho=10$, thus giving the analytical total mass as $7$ units. The numerical
mass closely matches this value throughout the simulation, with the error
bounded by the same amount as the deviation in the total volume.

The examples illustrates that the representative masses adaptation works well
while adding and deleting points. This also serves as a preliminary test to
show that no spurious effects are introduced by the volume correction
mechanism.

Figure 3: Adaptive refinement: Gradually refining the point cloud. The
computational domain is marked in black. The stationary fluid only fills a
part of the domain. The colour indicates the numerical volume of each point.
The evolution of the numerical volume and the representative mass are shown in
Figure 4.

Figure 4: Adaptive refinement test case: Evolution of the total mass
($\sum_{i=1}^{N}\hat{m}_{i}$) and volume ($\sum_{i=1}^{N}V_{i}$) in the
computational domain as the point cloud is refined (left), and the evolution
of the number of the points ($N$) in the domain (right).

###  7.2 Dam Breaking

We now consider a two dimensional dam break problem, which is a classical
benchmark test case for free surface flows. The dam break simulation is run
till the fluid comes to a rest. To prevent the influence of the numerical
definition of volumes, the final volume of the fluid is determined by the
height of the fluid at rest multiplied by the width of the domain. This is
then compared to the analytical initial volume.

The evolution of the numerical volume for the standard collocation approach
and the volume conservative method is shown in Figure 5. The standard approach
has a steady drop in numerical volume till the fluid starts coming to a rest,
with the biggest volume drop occurring when the fluid reaches the wall. On the
other hand, the volume conservative simulations only show minor fluctuations
of volume around the analytical volume.

Figure 5: Dam breaking test case: Evolution of numerical volume for the
collocation simulations, with and without the volume conservation algorithm.

###  7.3 Droplets Falling

While the earlier two test cases had two dimensional spatial domains,
henceforth, all examples considered are in $\mathbb{R}^{3}$. We now consider a
more complex case of rain-like droplets falling onto an inclined plate and
collecting in a bowl below the plate, as shown in Figure 6. Upon impinging on
the inclined plate, the droplets slide down and then collect in a cylindrical
container below. Each droplet considered is volumetrically resolved. The
inflow of droplets of a constant diameter occurs above the inclined plate at a
fixed rate. The height of the droplets above the inclined plate, and viscosity
of the fluid considered is chosen such that the droplet does not splash upon
impact. A test case with splashing and spraying fluid is considered later in
Section 7.5.

(a) (a)

(b) (b)

(c) (d)

(d) (c)

Figure 6: Droplets falling test case: Volumetrically resolved fluid droplets
falling onto an inclined plate and then collecting in a cylindrical container.
Image sequence is clockwise from top-left.

Each fluid droplet is spherical with diameter $0.15\text{m}$. A newly
generated droplet has no initial velocity, and falls due to gravity. The fluid
considered has density $\rho=1\text{kg/m}^{3}$, and dynamic viscosity
$\eta=0.051\text{Pas}$. This results in a Reynolds number of about $Re\approx
15$. Surface tension and contact angle hysteresis effects are not considered.
Homogeneous Neumann pressure boundary conditions are applied on all wall
boundaries. Slip boundary conditions are applied for the velocity at both the
inclined plate and the container.

Since fluid enters the domain at a fixed rate at the inflow, the total
geometrical volume of the fluid domain considered is known analytically at
every time step, $V_{\text{analytical}}=V_{\text{drop}}\sum_{\text{drops}}1$,
where $V_{\text{drop}}$ is the volume of single droplet, and the summation is
over the number of droplets injected into the domain until the current time.
Thus, the analytical volume expectation has a jump every time a new droplet is
added to the simulation domain.

The total numerical volume in the domain $V$ as compared to the analytical
expectation is plotted in Figure 7. The figure shows the evolution of the
numerical volume for both the standard simulations and the conservative one,
for two different resolutions. In both cases, the conservative simulation
shows a numerical volume very close to the expected one. On the other hand,
the standard simulation deviates significantly from the analytical
expectation. The biggest drop in volume occurs each time a droplet impinges on
the inclined plate, with a further decline in numerical volume after it
reaches the container.

(a) (a) $h=0.02$. Points per droplet $=3\,807$

(b) (b) $h=0.04$. Points per droplet $=29\,353$

Figure 7: Droplets falling test case: Numerical volume compared to analytical
volume. In both figures, the actual analytical volume (dashed green line) and
the numerical volume in the conservative simulations (dotted blue line) are
almost overlapping (see Table 2 for a quantified comparison), while the
numerical volume in the standard simulations (solid red line) deviates
significantly.

A quantification of errors in volumes (see Eq. Eq:Verr) for both the standard
and conservative simulations is tabulated in Table 2 for multiple resolutions.
The volume error for the conservative simulations is almost constant as the
resolution is made smaller. For the coarsest simulation considered, $h=0.04$
which corresponds to $619$ points per droplet, the volume error for the
standard case is $29\%$, while that for the conservative case is three orders
of magnitude smaller at $0.047\%$. Whereas for the finest simulation
considered, $h=0.01$ which corresponds to $29353$ points per droplet,the error
in the conservative case is two orders of magnitude smaller than that in the
standard case. The number of points needed to resolve each droplet for the
different resolutions used is listed in Table 2, along with the total number
of points in the simulation domain at the end of the simulation at $t=3$. Note
that the number of point at the end of the simulation varies slightly between
the standard and conservative simulations due to the addition and deletion
algorithms (see Sections 2.1 and 2.2).

Table 2: Droplets falling test case: Error in volume, $\epsilon_{V}$ (see Eq. Eq:Verr ) at $t=3$, number of points per droplet, and total number of points in the simulation domain at the end of the simulation ($t=3$) for the different resolution sizes considered. Resolution | Points per droplet | Points at $t=3$ |  $\epsilon_{V}$ standard |  $\epsilon_{V}$ conservative
---|---|---|---|---
$h=0.04$ | $\phantom{00}\,619$ | $3.1\times 10^{4}$ | $2.9\times 10^{-1}$ | $4.7\times 10^{-4}$
$h=0.02$ | $\phantom{0}3\,807$ | $1.8\times 10^{5}$ | $7.3\times 10^{-2}$ | $5.1\times 10^{-4}$
$h=0.01$ | $29\,353$ | $6.8\times 10^{5}$ | $2.6\times 10^{-2}$ | $4.8\times 10^{-4}$

###  7.4 Partially Filled Screw Conveyor

Screw conveyors are used for material processing, for example for polymer
extrusion. They use sophisticated systems with interacting screws and often
complex rheology is involved. Some systems are only partially filled with
material thus giving rise to a free surface between the material and
surrounding atmosphere. Meshfree methods come with a great advantage when
simulating such partially filled systems. However, volume loss can be a severe
issue especially at machine parts which are moving fast against each other
like a screw and its casing. This is a problem often seen in low resolution
meshfree simulations. This problem could be circumvented with very small
spatial and temporal resolutions. However, that would effectively make
simulations of larger systems impossible. Using the proposed volume correction
algorithm enables us to simulate large systems with moderate resolution
settings.

Motivated by such an industrial polymer extrusion problem, we demonstrate the
effects of volume correction using a simplified test case of a screw conveyor
as shown in Figure 8. The screw is rotating clockwise around the x axis
(pointing from left to right in the figure). The screw is tightly enclosed in
a non-rotating casing. Polymer melt is injected on the left. The inlet is
modeled as a small circular surface which rotates together with the screw. The
material is then transported towards the right side by the moving screw. For
simplicity, the material is modeled as a Newtonian fluid with constant
viscosity. For a real setting, more complex rheological models are needed. The
fluid used has a desity $\rho=10^{3}\text{kg/m}^{3}$, and dynamic viscosity
$\eta=200\text{Pas}$. The diameter of the screw is $0.04\text{m}$, which
results in a Reynolds number of about $Re\approx 0.2$.

Figure 8: Screw conveyor test case: Illustration of the domain. The fluid is
marked in green, and the fluid inlet is marked in blue on the left of the
domain. The screw rotates clockwise around its axis, while the casing is
fixed. As the screw rotates, the fluid is transported from the left to the
right side in the figure.

Figure 9: Screw conveyor test case: Evolution of extrusion without the volume
correction mechanism. The fluid is marked in green. This highlights a common
problem in meshfree methods where the fluid volume disappears when the
resolution used is not fine enough. In contrast, this problem is no longer
present in the volume corrected simulations, as shown in Figure 10.

Figure 10: Screw conveyor test case: Evolution of extrusion without the volume
correction mechanism. The fluid is marked in green. The issue of mass loss
present in the standard simulations (Figure 9) is no longer present here.

The evolution of the extrusion process is shown in Figure 9 without the volume
correction mechanism. Mass is lost quickly, mainly at the edge between screw
and casing which are moving against each other. The loss happens so fast that
after two or three revolutions no material is left. The material never reaches
the right side of the device. In order to have material reach the right side
for the setting without volume correction the resolution would have to be
chosen extremely small. The same setting with volume correction enabled is
shown in Figure 10. Both simulations use the same resolution. No mass loss is
observed for this setting. The material is transported by the screw and it
reaches the right side after some time.

###  7.5 Industrial Application: Automotive Water Crossing

We now consider an actual industrial application of a car moving through a
shallow pool of water, as shown in Figure 11. Such simulations have been
referred to as water crossing or water wading simulations. The primary goal in
these simulations is typically to check how much water reaches different parts
of the car, especially in the under-body of the car. An example of this is to
determine how much water, if any, enters the air intake where no water should
go. This test cases represents a much more complex geometry and a non-smooth
flow pattern than that considered in previous sections.

Figure 11: Automotive water crossing: A car crossing a shallow pool of water a
constant velocity. The car is shown with a slight transparency to improve the
visualization of the fluid.

This test case is a perfect example of the need for meshfree methods over
mesh-based ones. Meshing the complex under-body of the car can be quite
challenging. Moreover, since the car is moving through the water, the mesh
needs to be updated at every time step. The domain occupied by the fluid
spraying is much larger than the initial domain of the pool, as can be seen in
Figure 11. Rather than meshing the entire domain, a meshfree approach only
tracks points or particles where the fluid is actually present, with no
discretization away from the fluid. Furthermore, the fluid sprays to regions
where very small amount of fluid is present, which would require a very fine
mesh to capture accurately.

The ease of adaptive refinement near the car geometry in meshfree collocation
approaches give them an advantage over particle-type meshfree methods (see
Section 2.1). However, due the absence of a native notion of mass, it is not
possible to determine how much water truly goes into different parts of the
under-body. The introduction of the notion of representative masses in the
present work solves this problem.

The car geometry used is the open-source DrivAer car model [20] 111The estate
back configuration with the detailed under body is used.. The majority of the
geometry is treated as a rigid body moving at a fixed velocity. The wheels and
tyres rotate at a fixed angular velocity around their respective axle such
that it matches with the linear velocity of the car. Since the application is
intended towards cars crossing rain water, the fluid simulated is water with
$\rho=10^{3}\text{kg/m}^{3}$ and $\eta=10^{-3}\text{Pas}$. This results in a
Reynolds number of about $Re\approx 2\times 10^{7}$.

The evolution of the total fluid volume as the car crosses the water pool, for
both the standard and conservative simulations, is shown in Figure 12. Similar
to the previous test cases, the numerical volume in the standard simulation
drops as the simulation progresses, with the volume starting to increase once
the water spray starts increasing. On the other hand, the volume in the
conservative simulation remain approximately constant, with the fluctuations
increasing once the fluid spray increases. About a $2\%$ deviation in volume
is observed in the standard simulations, despite the use of fine resolution in
which the initial water pool contains about $2\times 10^{5}$ points.

Figure 12: Automotive water crossing: The total numerical volume of all points
in the simulation for the conservative simulations (dotted blue line), and the
standard simulations (solid red line). Since there is no inlet or outlet of
fluid in the simulations, the volume should remain constant. In the
conservative simulations, the volume is approximately constant, while it drops
and then increases in the standard simulation.

A visual comparison of the spray of water in both cases in shown in Figure 13.
The figure illustrates that the free surface spray patterns are similar in
both cases, with slight differences only towards the tail of the spray. This
illustrates that the actual difference made in the fluid profile is not
significant.

(a) Side view

(b) Bottom view with the water in the pool hidden

Figure 13: Automotive water crossing: Overlap of the standard simulation
(white) and the conservative simulations (blue). A slight transparency is
applied to both simulation results in both figures to enhance visualization.
The figure shows that the standard and conservative simulations show very
similar flow profiles.

##  8 Conclusion

In this work, we introduced a notion of mass, called representative mass, for
points in meshfree collocation type methods. This was shown to be essential
for proper post-processing, since a concept of mass is not inherently present
in these methods. This work also presented a first overview of different kinds
of conservation errors across different meshfree methods. We focused on the
issue of volume conservation, for which we highlighted the role of the
discrete definition of volume itself, while also giving an overview of
different volume definitions in meshfree methods.

Using the introduced representative masses, we computed a representative
density. Based on the difference between the representative and actual
densities, a novel volume correction mechanism was presented. Volume
conservation was acieved by introducing an artificial velocity divergence like
term to prescribe an artificial expansion or compression of the numerical
domain. Our numerical results across multiple test cases of varying complexity
show that the use of the volume conservation algorithm significantly
improvements conservative behaviour, without affecting the flow patterns in
any other way. We also showed that the introduced algorithm is very relevant
for actual industrial flow problems as well.

While the method presented here was in the context of collocation type
meshfree methods, it could be easily extended to be applied for particle type
meshfree methods like SPH. Here, the artificial movement was done based on the
relative difference between the representative density and the actual one. For
particle type meshfree methods, the volume correction could be prescribed
either in the same way, or based on the relative difference between the
representative mass and actual mass.

## Acknowledgements

Pratik Suchde would like to acknowledge support from the European Union’s
Horizon 2020 research and innovation programme under the Marie Skłodowska-
Curie Actions grant agreement No. 892761. Pratik Suchde and Stéphane P.A.
Bordas would like to acknowledge funding from the Institute of Advanced
Studies, University of Luxembourg, under the AUDACITY programme. Pratik Suchde
would also to thank Daniel Louis Louw for helpful discussions on FVM VoF, and
Mohsen Abdolazadeh for helpful discussions on SPH.

## References

  * [1] T. Arrufat, M. Crialesi-Esposito, D. Fuster, Y. Ling, L. Malan, S. Pal, R. Scardovelli, G. Tryggvason, and S. Zaleski.  A mass-momentum consistent, volume-of-fluid method for incompressible flow on staggered grids.  Computers & Fluids, 215:104785, 2021.
  * [2] M. Bašić, B. Blagojević, C. Peng, and J. Bašić.  Lagrangian differencing dynamics for time-independent non-newtonian materials.  Materials, 14(20):6210, 2021.
  * [3] J. Bašić, N. Degiuli, B. Blagojević, and D. Ban.  Lagrangian differencing dynamics for incompressible flows.  Journal of Computational Physics, 462:111198, 2022.
  * [4] J. Bašić, N. Degiuli, Šime Malenica, and D. Ban.  Lagrangian finite-difference method for predicting green water loadings.  Ocean Engineering, 209:107533, 2020.
  * [5] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, and P. Krysl.  Meshless methods: An overview and recent developments.  Computer Methods in Applied Mechanics and Engineering, 139(1):3–47, 1996.
  * [6] J. Benito, F. Ureña, L. Gavete, and R. Alvarez.  An h-adaptive method in the generalized finite differences.  Computer Methods in Applied Mechanics and Engineering, 192(5):735–759, 2003.
  * [7] A. S. Bharadwaj, J. Kuhnert, S. P. Bordas, and P. Suchde.  A discrete droplet method for modelling thin film flows.  Applied Mathematical Modelling, 112:486–504, 2022.
  * [8] S. Bordas, A. Menk, and S. Natarajan.  Partition of unity methods, 2022.
  * [9] D. L. Brown, R. Cortez, and M. L. Minion.  Accurate projection methods for the incompressible Navier–Stokes equations.  Journal of Computational Physics, 168(2):464–499, 2001.
  * [10] A. J. Chorin.  Numerical solution of the Navier–Stokes equations.  Mathematics of computation, 22(104):745–762, 1968.
  * [11] M. Depolli, J. Slak, and G. Kosec.  Parallel domain discretization algorithm for RBF-FD and other meshless numerical methods for solving pdes.  Computers & Structures, 264:106773, 2022.
  * [12] O. Diyankov.  Uncertain grid method for numerical solution of PDEs.  Technical report, NeurOK Software, 2008.
  * [13] C. Drumm, S. Tiwari, J. Kuhnert, and H.-J. Bart.  Finite pointset method for simulation of the liquid - liquid flow field in an extractor.  Computers & Chemical Engineering, 32(12):2946 – 2957, 2008.
  * [14] C.-M. Fan, C.-N. Chu, B. Šarler, and T.-H. Li.  Numerical solutions of waves-current interactions by generalized finite difference method.  Engineering Analysis with Boundary Elements, 2018.
  * [15] Y. Farjoun and B. Seibold.  Solving one dimensional scalar conservation laws by particle management.  In M. Griebel and M. A. Schweitzer, editors, Meshfree Methods for Partial Differential Equations IV, pages 95–109, Berlin, Heidelberg, 2008\. Springer Berlin Heidelberg.
  * [16] Fraunhofer Society.  MESHFREE.  https://www.meshfree.eu.
  * [17] D. García-Senz, R. M. Cabezón, and J. A. Escartín.  Conservative, density-based smoothed particle hydrodynamics with improved partition of the unity and better estimation of gradients.  Astronomy & Astrophysics, 659:A175, 2022.
  * [18] L. Gavete, M. Gavete, and J. Benito.  Improvements of generalized finite difference method and comparison with other meshless method.  Applied Mathematical Modelling, 27(10):831–847, 2003.
  * [19] R. A. Gingold and J. J. Monaghan.  Smoothed particle hydrodynamics: theory and application to non-spherical stars.  Monthly notices of the royal astronomical society, 181(3):375–389, 1977.
  * [20] A. I. Heft, T. Indinger, and N. A. Adams.  Introduction of a new realistic generic car model for aerodynamic investigations.  Technical report, SAE Technical Paper, 2012.
  * [21] P. F. Hopkins.  A new class of accurate, mesh-free hydrodynamic simulation methods.  Monthly Notices of the Royal Astronomical Society, 450(1):53–110, 2015.
  * [22] A. Iske.  On the construction of mass conservative and meshless adaptive particle advection methods.  In Advances in Meshfree Techniques, pages 169–186. Springer, 2007\.
  * [23] T. Jacquemin, P. Suchde, and S. P. Bordas.  Smart cloud collocation: geometry-aware adaptivity directly from cad.  Computer-Aided Design, 154:103409, 2023.
  * [24] T. Jacquemin, S. Tomar, K. Agathos, S. Mohseni-Mofidi, and S. P. Bordas.  Taylor-series expansion based numerical methods: a primer, performance benchmarking and new approaches for problems with non-smooth solutions.  Archives of Computational Methods in Engineering, 27:1465–1513, 2020.
  * [25] M. Jandaghian, A. Krimi, A. R. Zarrati, and A. Shakibaeinia.  Enhanced weakly-compressible MPS method for violent free-surface flows: Role of particle regularization techniques.  Journal of Computational Physics, 434:110202, 2021.
  * [26] A. Jefferies, J. Kuhnert, L. Aschenbrenner, and U. Giffhorn.  Finite pointset method for the simulation of a vehicle travelling through a body of water.  In M. Griebel and A. M. Schweitzer, editors, Meshfree Methods for Partial Differential Equations VII, pages 205–221, Cham, 2015. Springer International Publishing.
  * [27] Z. Ji, L. Fu, X. Y. Hu, and N. A. Adams.  A new multi-resolution parallel framework for SPH.  Computer Methods in Applied Mechanics and Engineering, 346:1156–1178, 2019.
  * [28] S. Koshizuka and Y. Oka.  Moving-particle semi-implicit method for fragmentation of incompressible fluid.  Nuclear science and engineering, 123(3):421–434, 1996.
  * [29] T. Kugelstadt, A. Longva, N. Thuerey, and J. Bender.  Implicit density projection for volume conserving liquids.  IEEE Transactions on Visualization and Computer Graphics, 27(4):2385–2395, 2019.
  * [30] J. Kuhnert.  Meshfree numerical scheme for time dependent problems in fluid and continuum mechanics.  In S. Sundar, editor, Advances in PDE Modeling and Computation, pages 119–136, New Delhi, 2014. Anne Books.
  * [31] E. Kwan-yu Chiu, Q. Wang, R. Hu, and A. Jameson.  A conservative mesh-free scheme and generalized framework for conservation laws.  SIAM Journal on Scientific Computing, 34(6):A2896–A2916, 2012.
  * [32] B.-H. Lee, J.-C. Park, M.-H. Kim, and S.-C. Hwang.  Step-by-step improvement of MPS method in simulating violent free-surface motions and impact-loads.  Computer Methods in Applied Mechanics and Engineering, 200(9):1113–1125, 2011.
  * [33] S. Lind, R. Xu, P. Stansby, and B. Rogers.  Incompressible smoothed particle hydrodynamics for free-surface flows: A generalised diffusion-based algorithm for stability and validations for impulsive flows and propagating waves.  Journal of Computational Physics, 231(4):1499–1523, 2012.
  * [34] S. J. Lind, B. D. Rogers, and P. K. Stansby.  Review of smoothed particle hydrodynamics: towards converged lagrangian flow modelling.  Proceedings of the royal society A, 476(2241):20190801, 2020.
  * [35] M. Luo, C. G. Koh, W. Bai, and M. Gao.  A particle method for two-phase flows with compressible air pocket.  International Journal for Numerical Methods in Engineering, 108:695–721, Nov. 2016.
  * [36] H.-G. Lyu and P.-N. Sun.  Further enhancement of the particle shifting technique: Towards better volume conservation and particle distribution in SPH simulations of violent free-surface flows.  Applied Mathematical Modelling, 101:214–238, 2022.
  * [37] I. Michel, T. Seifarth, J. Kuhnert, and P. Suchde.  A meshfree generalized finite difference method for solution mining processes.  Computational Particle Mechanics, 8(3):561–574, 2021.
  * [38] A. Möller and J. Kuhnert.  Simulation of the glass flow inside a floating process / Simulation de l’écoulement du verre dans le procédé float.  Revue Verre, 13(5):28–30, 2007.
  * [39] V. P. Nguyen, T. Rabczuk, S. Bordas, and M. Duflot.  Meshless methods: a review and computer implementation aspects.  Mathematics and computers in simulation, 79(3):763–813, 2008.
  * [40] S. B. Nguyen Vinh Phu, Alban de Vaucorbeil.  The Material Point Method: Theory, Implementations and Applications.  Springer, 2023.
  * [41] A. Obeidat, T. Andreas, S. P. Bordas, and A. Zilian.  Simulation of gas-dynamic, pressure surges and adiabatic compression phenomena in geometrically complex respirator oxygen valves.  Thermal Science and Engineering Progress, 24:100906, 2021.
  * [42] A. Obeidat and S. Bordas.  Three-dimensional remeshed smoothed particle hydrodynamics for the simulation of isotropic turbulence.  International Journal for Numerical Methods in Fluids, 2017.
  * [43] A. Obeidat and S. P. Bordas.  An implicit boundary approach for viscous compressible high reynolds flows using a hybrid remeshed particle hydrodynamics method.  Journal of Computational Physics, 391:347–364, 2019.
  * [44] G. Pahar and A. Dhar.  A robust volume conservative divergence-free ISPH framework for free-surface flow problems.  Advances in Water Resources, 96:423 – 437, 2016.
  * [45] G. Pahar and A. Dhar.  Robust boundary treatment for open-channel flows in divergence-free incompressible SPH.  Journal of Hydrology, 546:464 – 475, 2017.
  * [46] A. Petras, L. Ling, and S. J. Ruuth.  An rbf-fd closest point method for solving PDEs on surfaces.  Journal of Computational Physics, 370:43–57, 2018.
  * [47] X. Rao.  An upwind generalized finite difference method (GFDM) for meshless analysis of heat and mass transfer in porous media.  Computational Particle Mechanics, pages 1–22, 2022.
  * [48] T. R. Saitoh and J. Makino.  A density-independent formulation of smoothed particle hydrodynamics.  The Astrophysical Journal, 768(1):44, 2013.
  * [49] F. R. Saucedo-Zendejo.  A novel meshfree approach based on the finite pointset method for linear elasticity problems.  Engineering Analysis with Boundary Elements, 136:172–185, 2022\.
  * [50] F. R. Saucedo-Zendejo and J. M. Nóbrega.  A novel approach to model the flow of generalized newtonian fluids with the finite pointset method.  Computational Particle Mechanics, pages 1–11, 2021.
  * [51] B. Seibold.  M-Matrices in Meshless Finite Difference Methods.  PhD thesis, Kaiserslautern University, 2006.
  * [52] V. Shankar, G. B. Wright, R. M. Kirby, and A. L. Fogelson.  A radial basis function (RBF)-finite difference (FD) method for diffusion and reaction–diffusion equations on surfaces.  Journal of scientific computing, 63:745–768, 2015.
  * [53] L. D. G. Sigalotti, J. Klapp, and M. G. Gesteira.  The mathematics of smoothed particle hydrodynamics (SPH) consistency.  Frontiers in Applied Mathematics and Statistics, 7:797455, 2021\.
  * [54] J. Slak and G. Kosec.  On generation of node distributions for meshless PDE discretizations.  SIAM journal on scientific computing, 41(5):A3202–A3229, 2019.
  * [55] P. Suchde.  Conservation and Accuracy in Meshfree Generalized Finite Difference Methods.  PhD thesis, University of Kaiserslautern, Kaiserslautern, Germany, 2018\.
  * [56] P. Suchde.  A meshfree lagrangian method for flow on manifolds.  International Journal for Numerical Methods in Fluids, 93(6):1871–1894, 2021.
  * [57] P. Suchde, T. Jacquemin, and O. Davydov.  Point cloud generation for meshfree methods: An overview.  Archives of Computational Methods in Engineering, pages 1–27, 2022\.
  * [58] P. Suchde and J. Kuhnert.  Point cloud movement for fully lagrangian meshfree methods.  Journal of Computational and Applied Mathematics, 340:89–100, 2018\.
  * [59] P. Suchde and J. Kuhnert.  A fully lagrangian meshfree framework for PDEs on evolving surfaces.  Journal of Computational Physics, 395:38 – 59, 2019.
  * [60] P. Suchde, J. Kuhnert, S. Schröder, and A. Klar.  A flux conserving meshfree method for conservation laws.  International Journal for Numerical Methods in Engineering, 112(3):238–256, 2017.
  * [61] P. Suchde, J. Kuhnert, and S. Tiwari.  On meshfree GFDM solvers for the incompressible Navier–Stokes equations.  Computers & Fluids, 165:1 – 12, 2018.
  * [62] P. Sun, A. Colagrossi, S. Marrone, M. Antuono, and A.-M. Zhang.  A consistent approach to particle shifting in the $\delta$-plus-SPH model.  Computer Methods in Applied Mechanics and Engineering, 348:912–934, 2019.
  * [63] A. Tayebi and Y.-c. Jin.  Development of moving particle explicit (MPE) method for incompressible flows.  Computers & Fluids, 117:1–10, 2015.
  * [64] A. Tramecon and J. Kuhnert.  Simulation of advanced folded airbags with VPS-PAMCRASH/FPM: Development and validation of turbulent flow numerical simulation techniques applied to curtain bag deployments.  In SAE Technical Paper, Warrendale, PA, USA, 2013. SAE International.
  * [65] N. Trask, P. Bochev, and M. Perego.  A conservative, consistent, and scalable meshfree mimetic method.  Journal of Computational Physics, 409:109187, 2020.
  * [66] L. Veltmaat, F. Mehrens, H.-J. Endres, J. Kuhnert, and P. Suchde.  Mesh-free simulations of injection molding processes.  Physics of Fluids, 34(3):033102, 2022.
  * [67] H. Xia and Y. Gu.  Generalized finite difference method for electroelastic analysis of three-dimensional piezoelectric structures.  Applied Mathematics Letters, 117:107084, 2021.
  * [68] H. Xiao and Y.-C. Jin.  Development of explicit moving particle simulation method with applications.  Computers & Fluids, 235:105270, 2022.
  * [69] W. Zhan, X. Rao, H. Zhao, H. Zhang, S. Hu, and W. Dai.  Generalized finite difference method (GFDM) based analysis for subsurface flow problems in anisotropic formation.  Engineering Analysis with Boundary Elements, 140:48–58, 2022.
  * [70] C. Zhang, M. Rezavand, and X. Hu.  A multi-resolution SPH method for fluid-structure interactions.  Journal of Computational Physics, 429:110028, 2021.