The Finite Mass Method.
The finite mass method, a new Lagrangian method for the numerical simulation of gas flows, is presented and analyzed. In contrast to the finite volume and the finite element method, the finite mass method is founded on a discretization of mass, not of space. Mass is subdivided into small mass packets of finite extension, each of which is equipped with finitely many internal degrees of freedom. These mass packets move under the influence of internal and external forces and the laws of thermodynamics and can undergo arbitrary linear deformations. The method is based on an approach recently developed by Yserentant and can attain a very high accuracy.
Introduction
. Fluid mechanics is usually stated in terms of conservation
laws that link the change of a quantity like mass or momentum inside a given volume
to a flux of this quantity across the boundary of the volume. The finite volume
method is directly based on this formulation. Space is subdivided into little cells, and
the balance laws for mass, momentum and energy are set up for each of these cells
separately. Similarly, also the finite element method is based on a discretization of
space and a choice for the trial functions on the resulting cells.
In contrast, the finite mass method is founded on a discretization of mass, an idea
which is at least as obvious and that can be traced back to the work of von Neumann
[10] or of Pasta and Ulam [9] in the late 1940's and the 1950's. Instead of dividing
space into elementary cells, we divide mass into a finite number of mass packets
of finite extension, each of which equipped with a given number of internal degrees
of freedom. These mass packets move under the influence of internal and external
forces and the laws of thermodynamics and can intersect and penetrate each other.
They can contract, expand, rotate, and even change their shape. Their internal mass
distribution is described by a fixed shape function, similarly as with finite elements.
Although the finite mass method is a purely Lagrangian approach, it has not much
to do with particle methods as used for Boltzmann-like transport equations; in some
way, it is much closer to finite element and finite volume schemes. The approximations
it produces are differentiable functions and not discrete measures. The method is
basically of second order, and in some experiments we observed even fourth order
convergence!
The Lagrangian form of description of fluid flows can have many advantages. For
example, there are no problems with free surfaces, and no convection terms arise. Such
features make numerical methods based on the Lagrangian view especially attractive
unbounded space. In fact, one of the most popular methods of this
type, Monaghan's smoothed particle hydrodynamics [8], has its origins in astrophysics.
Like the smoothed particle hydrodynamics, the finite mass method is a completely
grid-free approach, but it is not a method de facto imitating statistical mechanics and
possesses a sounder mathematical and physical foundation.
The finite mass method is a generalization and extension of the particle model of
compressible fluids that had been proposed by the last author in [11], [12] and [13] and
is based on the principles developed there. The compactness and convergence results
Mathematisches Institut der Universitat Tubingen, 72076 Tubingen, Germany
obtained in [11] and [13] concerning the transition to the continuum limit transfer to
the present situation. One of the essential differences to the approach in the articles
mentioned above is that the single mass packets can now undergo arbitrary linear
deformations and not only rotations and changes of size. Although the dimension of
the configuration manifold of the single mass packet increases because of that, this
strongly simplifies the equations of motion the mass packets are subject to because
the configuration manifold is now a linear space. The equations of motion take the
same form for all space dimensions. The main advantage, however, are the superior
approximation properties, due to the fact that the mass packets can now be deformed
by the flow.
The main issue with a method like ours is how the equations of motion for the
mass packets or particles, as we often prefer to say, are set up. We start from the basic
physical principles that finally lead to the Euler and Navier-Stokes equations, not from
these equations themselves and a least squares or Galerkin approach. In consequence,
the equations of motion also do not break down when particles completely cover
each other. In the most simple case of an adiabatic, inviscid flow, the equations
of motion for the particles are derived from a Lagrange-function with the internal
energy as potential energy. To damp the fluctuation part of the local kinetic energy
that necessarily arises with every such model, frictional forces vanishing in the limit
of particle sizes tending to zero are added to these potential forces. They break the
invariance to time reversal and make the method consistent with the second law of
thermodynamics.
The paper is organized as follows. In x2, the finite mass method is explained
and derived in detail, both for inviscid and for viscous fluids. It is shown how the
mass density, the velocity field and, as second independent thermodynamic quantity,
the entropy are discretized, and it is discussed what sort of approximation properties
can be expected. The equations of motion describing the local interaction of
the mass packets and the time evolution of the system are set up. A main feature
of the approach is that mass, momentum, angular momentum and energy are exactly
conserved. The conservation of mass follows immediately from the construction
and is discussed in x2. The conservation of energy, momentum and angular momentum
is studied in x3. The conservation of angular momentum is a remarkable fact
as in continuum mechanics the conservation of angular momentum is hidden in the
symmetry of the stress tensor and does not appear explicitly as a conservation law.
Consequently, most discretizations violate this principle.
The finite mass method as described in x2 is invariant to arbitrary translations
and rotations and, as it concerns the shape the mass packets can attain, even to
every linear transformation of space. These properties are reflected by the quadrature
formulas that are developed in x4 to evaluate the integrals which define the forces
acting upon the particles. Conservation of momentum, angular momentumand energy
are maintained. In x4, we also discuss how the forces and moments acting upon the
particles can be calculated on the computer.
In x5, a suitable time discretization is presented. The proposed scheme is an
exponential integrator in the spirit of the recent paper [4] by Hochbruck and Lubich.
In case of pure pressure forces, it transfers to the well-known Stormer-Verlet method
for second order equations and conserves momentum and angular momentum exactly.
Finally, in x6, some typical test calculations for flows in two space dimensions are
documented. These examples clearly demonstrate the potential and the high accuracy
of the method.
We restrict our attention in this article to free flows in vacuum. For the inviscid
case, the reflection laws needed for particles touching the walls of bounded volumes
have already been given in [11], [12] and [13]. Their numerical realization will be
discussed elsewhere.
For a certain background in mechanics and fluid dynamics, we refer to the text-books
[5], [6] by Landau and Lifschitz on one hand and [1] by Chorin and Marsden
on the other hand, and to the classical monograph [2] by Courant and Friedrichs.
2. The particle model of compressible fluids. The basic ingredient of the
finite mass method is a continuously differentiable shape function / : R d ! R, d the
space dimension, with compact support that attains only values  0. This function
describes the internal mass distribution inside the mass packets into which the fluid
is subdivided. We assume that
Z
Z
The second property states that the origin of the body coordinate system attached to
a particle is its center of mass. Further we suppose that
Z
with the y k the components of y.
For example, the function / may be built up from a piecewise polynomial function
e
in one space variable. If we let
d
Y
e
the conditions (2.1) are equivalent to
Z
e
Z
e
The second condition in (2.4) also implies that the integrals (2.2) vanish for k 6= l.
The constant J is given by
Z
e
and does not depend on the space dimension. A suitable choice for e
/, that we have
used in the numerical computations documented in x6, is the normalized third order
B-spline given by
e
for jj  1 and by e
1. It can be composed of smaller copies and
be used to build up a basis for the cubic spline functions on a uniform grid. This e
fulfills the conditions (2.4), and the constant (2.5) takes the value
An advantage of such tensor product like shape functions are their good approximation
properties.
The points y of the particle i move along the trajectories
The vector q i (t) determines the position of the particle and the matrix H i (t) its size,
shape and orientation in space. Correspondingly,
are the body coordinates of the point at position x in space at time t. Let
denote the mass of the particle i. The total mass density
then results from the superposition of the mass densities of the single particles.
The points y of the particle i have the velocity
Inserting the expression (2.9) for y, one gets the velocity field
of the particle i related to the space coordinates. The total mass flux density
again results from the superposition of the mass flux densities of the single particles.
With the local mass fractions
the velocity field v of the flow, that is defined by the relation
is the convex combination
of the velocity fields of the single particles.
To simplify notation, we introduce the abbreviation
where again
The gradient of / i , that is the derivative of / i with
respect to x, is
and can, as the derivatives
be expressed in terms of the internal variables y.
To derive the expression above for the derivative with respect to the matrix H i ,
we suppress the index i serving for the distinction of the particles and write /   instead
of / i . First we observe that, for H fixed and E tending to zero,
E)
and, because of
for A tending to zero, similarly
E)
where
k;l
denotes the inner product of two matrices A and B. One findsdet(H +E) /
det H (r/)(y)
This means that the linear mapping
det H /(y)
is the total derivative of
det H /
at given H. Inserting Ej l for the entries of E and using (2.17), (2.18),
@/
follows.
The continuity equation expressing the conservation of mass is automatically
satisfied with the given ansatz. Independently of the size and shape of the particles
and their distribution in space,
@ae
@t
This follows from the fact that such an equation holds for every single particle, as one
proves using (2.19) and div v
6 CHRISTOPH GAUGER, PETER LEINEN, HARRY YSERENTANT
To study the potential accuracy of the approach, we start from a given twice continuously
differentiable velocity field u. Then fixing the particle trajectories q i (t) and
the matrices H i (t) for t  t 0 , with given initial values, as solutions of the differential
equations
the velocity field (2.12) of the particle i reads
and is therefore a second order approximation of u in a neighborhood of
As the mass fractions (2.14) form a partition of unity, the resulting overall velocity
field
remains a second order approximation of u on the region occupied by mass, independent
of where the particles are located. The mass density is an exact solution of the
transport equation (2.20) with respect to this perturbed velocity field and therefore,
with corresponding initial values, a good approximation of the true density in the
sense of a backward error analysis. For a more careful and detailed analysis of this
type also covering the motion in an external force field, we refer to [14].
For a complete description of the thermodynamic state of a compressible fluid,
besides the mass density ae a second thermodynamic quantity like temperature is
needed. Most convenient for our purposes is the entropy density s that is given here
in the form
where the S i (t) denote the specific entropies of the single particles. In particular, the
specific entropy
is the convex combination
of the specific entropies of the single particles. As with the velocity, one recognizes
that this representation potentially leads to a first order approximation of the exact
specific entropy independent of the distribution of the particles, where the actual
accuracy can again be much higher.
The pressure, the absolute temperature and the internal energy per unit volume
are functions of the mass density and the entropy density. These functions are not
independent of each other but are connected by the Gibbs fundamental relation of
thermodynamics taking the form
@ae
@s
@s
in the present variables. We assume that the internal energy "(ae; s) is defined for all
ae ? 0 and all real s and is twice continuously differentiable on this set. We suppose
that
@"
@s
for all these ae and s and require that "=ae and the first order partial derivatives of "
can be extended by the value 0 at (ae; to functions that are continuous on
the sectors jsj  S   ae, S   ? 0. Barytropic ideal gases with
ae 0
s
ae
represent a simple example. The constant  0 ? 0 is a characteristic pressure, the constant
mass density, and the constant c v ? 0 the characteristic
heat for constant volume. The constant fl ? 1 is the ratio of the specific heats for
constant pressure and constant volume. A typical value is
The particle i has the kinetic energy
that attains the closed representation
denotes the Euclidean norm of a vector and the Frobenius norm of a matrix,
respectively. The total kinetic energy of the system is
From the kinetic energy (2.33) and the internal energy
Z
the Lagrange-function
of the system is formed. For the completely adiabatic case that is described by the
Euler-equations of gas dynamics, the equations of motion
d
dt
d
dt
are derived from this Lagrange-function. With the normalized forces
explicitly given by
@ae
@"
@s
@ae
@"
@s
these equations read
Note that the masses m i of the particles cancel out and appear only implicitly in
the mass density ae and the entropy density s. The derivatives of / i occurring in the
expressions above can be written as functions of the internal particle coordinates (2.9)
and are given by (2.19). The pressure forces (2.37) acting upon a group of particles
can be interpreted as a kind of surface force [11].
As long as the specific entropies S i of the single particles are assumed to be
constant, the time evolution of the system is completely determined by the equations
of motion (2.36). This system is time reversible, a fact that contradicts the behavior of
actual fluids where heat is generated in shock fronts. Therefore it has been proposed
in [13] to add the (normalized) frictional force
F (r)
Z
to the right-hand side of the first equation in (2.40). This force damps local velocity
fluctuations and couples the particles softly together. The difference between the
velocity of the given particle and the velocity field of the surrounding flow determines
the direction of the force and the scalar function R  0 its size. The quantity in the
equation for the H i corresponding to the force (2.41) is
In smooth flows, the forces (2.41) and (2.42) will be comparatively small and will
vanish with the second power of the local particle size, but they can dominate in
shocks.
For a better quantitative understanding of the frictional forces (2.41), (2.42), we
pause for a moment and consider a little model problem. We neglect the pressure
forces, keep R constant, and fix the velocity field The equations of motion for
a single particle then read
that is, within the time the velocity of the particle is reduced by the factor
1=e. This justifies calling the local relaxation time of the system. With
the trajectory of the particle is
Thus the particle is stopped in distance T jv 0 j from its position at the given time t 0 .
The friction among the particles generates heat. Therefore the specific entropies
are no longer constant and increase in time. They obey the differential equations
The quantity ' i defined by
is the mean value
Z
of the absolute temperature around the particle and
Z
the specific heat supply. The quantity
is the density of the fluctuation energy. It can be rewritten as
which is more convenient for computational purposes. The equation (2.43) is the
second law of thermodynamics. It ensures conservation of energy as will be shown in
x3. Note that the specific entropy of a particle never decreases.
If the q i and H i and the function R are kept fixed,
is a linear mapping. To study this mapping, we utilize the energy norm given by
and the corresponding energy inner product, respectively. Note that the square of
the energy norm of the vector consisting of the velocity components q 0
i is the
total kinetic energy (2.33) of the system. Denoting by b
v i the velocities (2.16) formed
with b
i instead of q 0
, one has
Z
with the symmetric bilinear form
in the velocities v i and b v i . This shows that the linear mapping (2.49) is negative
semidefinite and selfadjoint with respect to the energy inner product. For R constant,
its eigenvalues range in the interval [\Gamma1=T; 0]. The proof of (2.51) is based on the
same arguments as the proof of Theorem 1 below.
In viscous fluids, an additional force
Z
acts upon the mass contained in a volume W t transported with the flow, that is
the mass initially contained in a given volume inside the region occupied
by mass. The symmetric tensor T is the viscous part of the stress tensor and is a
function of the density, the temperature and of the symmetric part
of the gradient of the velocity field. In Newtonian fluids, it depends linearly on D
and v, respectively, and is given by
d
in d space dimensions. The first term on the right-hand side has trace zero and
describes the shear forces in the fluid. The second part comes from the bulk viscos-
ity. The viscosity coefficients j and i are nonnegative functions of ae and ', which
guarantees
a fundamental property needed to satisfy the second law of thermodynamics. For
most gases, can be assumed.
To model the viscous forces, one starts from the observation that
Z
Z
where the  ffi are continuously differentiable functions with values between 0 and 1
vanishing outside W t and tending pointwise to the characteristic function of W t for ffi
tending to 0. If I denotes the set of indices of the particles initially contained in the
region W , the mass fraction
i2I
serves as our discrete counterpart of the functions  ffi . This leads to the viscous force
Z
acting upon this group of particles, or to the normalized force
F (v)
Z
acting upon the particle i, where the integrals extend over the region occupied by
mass at the given time t. Note that the expression (2.59) correctly takes into account
only that part of the boundary of the moving volume that is touched by mass from its
exterior. Therefore the construction also covers the case not directly included above
when a part of the boundary of the moving volume belongs to the boundary of the
region occupied by mass. The quantities M i corresponding to the forces (2.60) are
Z
Z
and the additional specific heat supply to the particle i is
Z
Provided the given particle is located in the interior of the region occupied by mass,
the forces (2.60) and (2.61) formally vanish for stress tensors T of divergence zero, as
continuum mechanics requires. The proof is by partial integration, where the second
integral on the right-hand side of (2.61) cancels.
The problem with this approach is that, in more than one space dimension, the
derivatives of the mass fractions  i can have singularities and are not always square
integrable. Fortunately, this effect is often compensated by the behavior of T and the
properties of the shape function /. Assume that T is of the form (2.55) and that the
viscosity coefficients behave like
with ff ? 0. Let the gradient of / satisfy the pointwise inverse estimate
holds for the cubic splines (2.3), (2.6) with 3. Then there exist
constants C i depending on the number of the locally interacting particles and their
mass, size, and shape with
ae ff jr
Provided that ff \Gamma 2=k ? 0, this demonstrates that the integrands above are continuous
and tend to zero at the points at which ae tends to zero. Therefore, under the given
circumstances, one can refrain from replacing the  i by regularized mass fractions
such as proposed in [11] and [12]. However, our considerations would immediately
transfer to this case.
To incorporate heat conduction, the right-hand side of the entropy equation (2.43)
has to be supplemented to
Z
The heat flux vector k is mostly given by Fourier's law
where the coefficient function   0 is the heat conductivity. For the heat flux to
be well defined, k has to satisfy an estimate similar to the estimate above for the
viscosity coefficients.
For the rest of this section, we restrict ourselves to Newtonian fluids (2.55). For
given q i and H i , the mapping
is then linear. Denoting by b
D the tensor (2.54) formed with the velocities b q 0
instead of q 0
Z
Ddx
holds, with the integrals discretized as described above. Because
D \Gammad
(tr D)(tr b
this proves that also the linear mapping (2.68) is symmetric and negative semidefinite
with respect to the energy inner product given by (2.50). For the proof of (2.69), we
refer again to the proof of Theorem 1 below.
3. The conservation of energy, momentum, and angular momentum.
The conservation of energy, momentum, and angular momentum are basic physical
properties of any closed system and must therefore be reproduced by the finite mass
method. Moreover, energy estimates are basic for the mathematical examination of
the model and, in particular, are needed to transfer the compactness and convergence
results from [11] and [13] to the present situation of particles underlying arbitrary
linear deformations. The conservation of energy also prevents the determinants of
the H i from becoming arbitrarily small since, with equations of state like (2.30), this
would require too much energy. Our first result is:
Theorem 1. The total energy composed of the kinetic energy (2.33)
of the particles and the internal energy (2.34) is a constant of motion.
Proof. Utilizing the representation (2.32) of the kinetic energy of a single particle,
one first obtains
d
dt
Therefore the equations of motion yield
d
where the terms coming from the pressure forces (2.37) have canceled with the corresponding
derivatives of V ,
comes from the frictional forces (2.41), (2.42), and
from the viscous forces (2.60), (2.61).
Before we start calculating
3 , we state two simple algebraic relations
that are repeatedly used in this section, namely that
A
for all square matrices A and B and all vectors a and b and that
A
for all square matrices A, B and C.
For example, with I the first of these two relations yields
for arbitrary vectors f . Choosing
Z
follows. Because, by the definition (2.16) of v,2
with (2.48) this leads to the representation
Z
Rq dx
of the first of the two terms above. Inserting Tr i for f in (3.3), one finds
Z
for the second part coming from the viscous forces (2.60) and (2.61). As
by the symmetry of T and
by the definition of v i and (3.2), the second integral in the formula above cancels and
the sum attains the value
Z
14 CHRISTOPH GAUGER, PETER LEINEN, HARRY YSERENTANT
As by (2.43), (2.44), (2.46), (2.62) and (2.66)
Z
Z
this proves the proposition and demonstrates that exactly the right amount of kinetic
energy is converted to heat.
Surprisingly, for vanishing frictional and viscous forces, one has both conservation
of energy and entropy, which contradicts the usual conception of gas flows and would
not be possible in continuum mechanics, but which is explained by the fact that the
kinetic energy (2.33) considered here is composed of the kinetic energies of the single
particles and is not identical with the mean kinetic energy
E(t) =2
Z
known from continuum mechanics; it consists additionally of the local fluctuation
energy
Z
ae
In smooth flows, this fluctuation energy is a negligibly small part of the total kinetic
energy and will vanish with the fourth power of the particle size, but it can dominate
where the particles clash. The role of the frictional forces (2.41) and (2.42) is to
convert this kind of fluctuation energy into true internal energy.
The total momentum of the system is defined as
Z
Because of (2.1), it has the closed representation
Theorem 2. The total momentum (3.6) is a constant of motion.
Proof. The representation (3.7) and the equations of motion yield
d
dt
0with the three parts
coming from the pressure forces (2.38), the frictional forces (2.41), and the viscous
forces (2.60). With (2.19), the first part reads
@ae
@"
@s
By the definitions (2.10) of the mass density and (2.24) of the entropy density, this
means
Z
As ", under the given assumptions, is a continuously differentiable function with
compact support, this yields P 0
2 , with (2.41) one obtains
Z
As, by the definition (2.16) of the velocity field v,
As the  i form a partition of unity, finally also the part
Z
resulting from the viscous forces vanishes.
Last, we consider the scalar quantities
Z
ae(x;
with fixed skew-symmetric matrices W, that have, by the definition of v and because
of (2.1) and (2.2), the closed representation
In three space dimensions, the components of the total angular momentum
Z
are quantities of this form, and vice versa all quantities of the form (3.8) can be
composed of the components of the angular momentum (3.10). Therefore, for the
three-dimensional case, the following theorem states that the angular momentum of
the system is a constant of motion. In other space dimensions, one gets less or more
first integrals, corresponding to the dimension of the space of the skew-symmetric
matrices.
Theorem 3. For arbitrarily given skew-symmetric matrices W, the scalar quantity
L defined by (3.8) is a constant of motion.
Proof. As a \Delta vectors a and A \Delta
one gets
dL
dt
from the representation (3.9) of L. The equations of motion therefore yield
dL
0where the first part
comes from the pressure forces (2.38) and (2.39), the second part
from the frictional forces (2.41) and (2.42), and the third part
from the viscous forces (2.60) and (2.61).
We show that each of these three parts vanishes separately. From the representation
of the partial derivatives of / i , (3.1) and (3.2), and the skew-symmetry
of W, first
follows. Thus (2.38) and (2.39) yield
@ae
@"
@s
or, again taking into account the definitions of ae and s, summed up
Z
As " is a continuously differentiable function with compact support, partial integration
gives
and therefore L 0
because W is skew-symmetric. In the same way, (3.1) yields
and therefore, with (2.41) and (2.42),
Z
As, by the definition (2.16) of the velocity field v,
one obtains L 0
(3.2), the symmetry of the viscous stress tensor T and the
skew-symmetry of W give
and therefore as above
Z
As the  i form a partition of unity, finally also L 0
4. The discretization of the integrals. Our particle model of compressible
fluids is purely Lagrangian. It is invariant to arbitrary translations and rotations and,
as it concerns the shape the particles can attain, even to every linear transformation
of space. These properties must be reflected by the quadrature formula needed to
evaluate the integrals that define the forces acting upon the particles.
We start from the observation that the integral of a scalar or vector function f
weighted by the mass density (2.10) can be written as a sum
Z
Z
of contributions from the single mass packets. The integrals on the right-hand side of
(4.1) live on the reference domain on which the shape function / is strictly positive.
They are replaced by a fixed quadrature formula
Z
ff  g(a  )
with weights ff  ? 0 and nodes a  inside the support of the shape function /, which
is considered as the weight function here. This results in the composite rule
Z
Z
f d
with weights m j ff  and nodes q j +H j a  that is based on the given discretization of
mass and not on a subdivision of space into elementary cells. Therefore it has the
desired invariance properties.
Utilizing the quadrature rule (4.3), one replaces the potential energy (2.34) by
the fully discrete expression
Z
denotes the specific internal energy. With this discrete potential, the
normalized forces
(4.
and the temperatures
are formed. As both the specific internal energy itself and the quadrature points
depend on the q i and H i , the forces F i and M i correspondingly split into two parts
that are of different structure. Into
F (1)
@ae
@~"
@s
@ae
@~"
@s
d
all quadrature points q j +H j a  contained in the support of the given particle i enter
whereas
F (2)
and the corresponding term
M (2)
are completely determined by the function values and the first order derivatives of
the mass density and the entropy density at the quadrature points q i +H i a  assigned
to the particle i itself. The local temperature (4.6) around the particle i transfers to
Z @~"
@s
The forces (4.7) replace the forces (2.38) and (2.39) in the Lagrangian equations of
motion (2.36) and (2.40), respectively, and the local temperature (4.6) the local temperature
(2.44) in the entropy equation (2.43). Note that, through the assumptions
we made on "(ae; s) and because
ff
the quantities above behave numerically well. The frictional forces (2.41) and (2.42)
and the heat supply (2.46) are directly discretized using the quadrature rule (4.3).
For the viscous forces (2.60) and (2.61), the heat supply (2.62) and the heat flux in
(2.66) one can proceed correspondingly.
The method reacts sensitively on the choice of the weighted quadrature rule (4.2)
for the reference domain. Experience has shown that quite a few quadrature points
are needed to exploit the full accuracy of the approach. A too small number of
quadrature points leads to instabilities, in particular when the quadrature points are
not properly spaced; a high polynomial accuracy alone does not suffice. For the

Table
The basic one-dimensional quadrature rule for cubic B-splines
nodes 4131656631641tensor-product third order B-splines described at the beginning of x2, we had good
experience with the tensor-product counterpart of the one-dimensional quadrature
rule given by Table 1. This quadrature formula is exact for fifth order polynomials
and assigns 5 2 or quadrature points to each particle in two space dimensions.
The total number of quadrature points q j +H j a  entering into the computation for
a given particle i depends on the overlap of the B-splines and will generally be much
higher, not much less than in two space dimensions.
The conservation of energy, momentum, and angular momentum transfers to the
present case of discretized integrals. The proofs from the last section can be taken
over with minor changes concerning only the pressure terms. For the proof that the
total momentum (3.6) remains a constant of motion, one has only to utilize that
Z
Z
r~" d
such that the quantity
vanishes. In the proof that the quantities (3.8) associated with the angular momentum
are constants of motion, the relations
Z
Z
enter, giving L 0
The essential idea to compute the forces acting upon the particles is to separate
the quadrature points from the particles. This reflects the fact that particles do not
interact directly with each other but only with global fields like the mass density or
the velocity.
Two different data structures are used. The first data structure is associated with
the particles themselves. It contains the fixed particles masses m i , the values q i , H i ,
i and S i finally to be determined, storage for the local temperatures ' i , and for
the forces acting upon the particles and the heat supply to the particles, of course.
The second data structure is associated with the quadrature points. First, it contains
their positions and the weights and then all necessary
global information at xZ like the values ae Z of the mass density, s Z of the entropy
density, j Z of the mass flux density or vZ of the velocity, quantities associated with
the viscous stress tensor T, should the occasion arise, and the other field information
needed in the computation.
The procedure to compute the integrals involved in the differential equations is
then quite simple and consists of three phases. We will describe this procedure only
for the inviscid case; the case of additional viscous forces is similar.
In the first phase, the quadrature points are generated and the values ae Z of the
mass density (2.10), s Z of the entropy density (2.24), j Z of the mass flux density
of the intermediate quantity
as well as the gradients of ae and s at the quadrature points are assembled. This is
a loop on all particles. In this phase, information is transferred from the particles to
the quadrature points.
In the second phase, the values of the derivatives
@ae
@s
and, to incorporate the frictional forces (2.41) and (2.42), of
ae Z
are computed. These operations work only on the single quadrature points.
In the final third phase, the discrete integrals determining the forces, the heat
supplies and the local temperatures are computed using the results of phase 1 and
phase 2. In this phase, information is transferred from the quadrature points back to
the particles.
Phase 1 and phase 3 require an efficient access to the quadrature points contained
in the support of a given particle. Search trees can be used for this purpose, which
have to be set up after the quadrature points have been generated. The information
needed to access the quadrature points assigned to a given particle i has
to be stored separately.
5. A time-stepping procedure. The remaining big system of differential equations
for the particle positions q i , the deformation matrices H i and the entropies S i
has to be solved numerically. In this section, we present a simple, robust second order
method for that which is adapted to the structure of this system and the properties
of its solutions. The method has been proposed to us by Christian Lubich [7] and is
inspired by the recent work of Hochbruck and Lubich [4] on exponential integrators.
Its stability range is de facto independent of the strength of the frictional and viscous
forces but shrinks, as with all schemes for this type of equation, when approaching
incompressibility.
First, we combine the vectors q i and the matrices H i to a big vector
of dimension (d and the S i to a vector z of dimension N . The system of
differential equations (2.40) together with the additional frictional forces (2.41) and
(2.42) and the viscous forces (2.60) and (2.61) can then be written in the form
and the entropy equation reads
The first term on the right-hand side of (5.2) corresponds to
with the F i and M i the discrete counterparts (4.7) of the pressure forces (2.38) and
(2.39). The second term
corresponds to the discretized versions of the frictional forces (2.41) and (2.42) and
of the viscous forces (2.60) and (2.61) in Newtonian fluids (2.55). The notation (5.5)
reflects that these forces depend only linearly on q 0
. As the considerations in
x2 have shown, the matrices are symmetric negative semidefinite with
respect to the energy inner product determined by (2.50). Finally, the right-hand side
of (5.3) corresponds to
with the ffi Q i the discretized versions of (2.46) and (2.62), respectively, plus eventually
a term coming from heat conduction. This function depends nonlinearly on all its
arguments.
The proposed method for the numerical solution of the system (5.2), (5.3) is a
leap-frog scheme. To come from the approximations for y 0 at time t k \Gamma =2 and for y
and z at time t k to the new approximations at times t k +=2 and t k respectively,
one first computes
The new value for the derivative of y is then
and the new value for y itself is
where OE 0 is given by
\Theta OE 2

and OE 1 and OE 2 are the entire functions
The new value z k+1 finally is implicitly determined by the equation
z
22 CHRISTOPH GAUGER, PETER LEINEN, HARRY YSERENTANT
Using the new b k+1 from the next time step, the approximation
is obtained at little additional cost. This approximation does not enter into the further
computations. The computation starts with the approximation
for y 0 at t 0 +=2 and the approximation
for y itself at the quantities z 1 and y 0
are then computed as above.
The method is constructed such that it would yield the exact values of y and y 0
for f and A constant. For vanishing frictional and viscous forces, i.e. for
the scheme transfers to the time-reversible integrator
known as Verlet method in molecular dynamics, where (5.13) reduces to
With the starting step above, the Verlet method is equivalent to the one-step method
k+1=2 is now considered as the intermediate value
The Verlet method is symplectic and has therefore very favorable properties for Hamiltonian
systems as arising in the present case of pure pressure forces [3]. In particular,
it preserves momentum and angular momentum: Let
denote the total momentum (3.7) and the scalar quantity (3.9).
As\Omega is a skew-symmetric
matrix and
by Theorem 2 and Theorem 3 and the considerations in x4, respectively, both quantities
retain their value in the transition from one time level to the next.
The matrix-vector products OE( A)b can be computed by a Krylov-space method
based on the symmetric Lanczos process using the energy inner product. Started with
the vector b, after m steps the Lanczos algorithm delivers a matrix V consisting of
columns orthonormal with respect to this inner product, a symmetric tridiagonal
matrix H of dimension m \Theta m, and a rank one matrix R such that
The matrix-vector products are then approximated by
With a spectral decomposition
of the small symmetric tridiagonal matrix H, OE( H) is given by
Note that the computation of Ax corresponds to an evaluation of the discretized
frictional forces (2.41) and (2.42) and, if necessary, of the discretized viscous forces
(2.60) and (2.61) and can be performed without any explicit knowledge of A. The
approximation (5.23) converges very fast for entire functions OE, much faster than the
conjugate gradient method for the solution of a linear system with the coefficient
Usually, very few Krylov steps suffice.
In the examples presented in the next section, we solved the equation (5.12)
approximately by two steps of a simple fixed point iteration. More sophisticated
methods will be needed when heat conduction is present.
6. Examples. In this section, we have compiled some two-dimensional examples
that exhibit a typical behavior and compare the numerical with the known exact
solutions. As shape function / of the particles we have used the bicubic B-spline
(2.3), (2.6) together with the 25-point quadrature rule described in x4.
To find a good approximation of the given initial data, we assume that the region
occupied by mass at time contained inside an axiparallel rectangle and cover
this rectangle by a regular grid of gridsize h. The initial matrices are then H i
and the initial positions are selected from the gridpoints, that we denote by x k here.
To determine the particle masses, we first compute the coefficients a k  0 minimizing
an appropriately chosen distance between the linear combination
a k
and the given initial mass density ae, for example by the projected Gau-Seidel
method. At the points x k for which a k ? 0, particles with are located.
The remaining gridpoints are ignored in the sequel.
The initial velocities can be determined by solving the interpolation problem
at the particle positions q Setting
the discrete mass flux will virtually be a fourth order approximation to the continuous
mass flux aev at time sufficiently smooth. The specific entropies S i (0) are
computed in the same way interpolating the entropy density s. The alternative is to
choose the initial velocities correspondingly to (2.21) as
which works better with dominating external forces and has the advantage that velocity
fields linearly depending on x like rigid body motions are exactly reproduced.
For reference, we recall the Navier-Stokes equations
@ae
@t
ae
@t
@s
@t
that govern the time evolution of a smooth compressible flow in the absence of heat
conduction and that transfer to the Euler equations for the inviscid case
These equations have to be completed by material laws like (2.30) and (2.55) and the
basic thermodynamic relation (2.28). Often, the equations are stated in a different,
but mathematically equivalent form, with (6.6) also written as conservation law and
replaced by the conservation law for the energy. For the special case of inviscid
ideal gases (2.30), the entropy equation (6.7) can be replaced by the equation
@
@t
coupling pressure and velocity.
6.1. Gas clouds. The first example serves, more or less, as consistency check
and shows how accurate the method can be. We reproduce the self-similar solutions
of finite mass of the Euler equations. The functions (6.9), (6.10) satisfy the continuity
equation (6.5) regardless of the choice of the density profile \Psi and the matrix function
H(t), and of course also the entropy equation (6.7). Starting from the internal energy
of a barytropic ideal gas (2.30), the momentum equation (6.6) reads
ae 0
has been set and H is an abbreviation for the determinant of H.
Up to normalization,
ae
is the only function with compact support that satisfies the relation (6.11) for an
appropriately chosen function H(t). For 1! 2, the function (6.12) is continuously
differentiable at 1. The corresponding matrix functions H(t) are the solutions of
the differential equation
ae 0
det H(t)
with initial values det H(0) ? 0. This differential equation coincides with the differential
equation from (2.40) one would obtain for the motion of a single
with the shape function (6.12).
If one lets the particles move according to the differential equations (2.21) in the
velocity field (6.10), the whole configuration would simply undergo a series of linear
transformations. So it is no wonder that a high accuracy can be reached in this
example. In the computation presented here, we started from the material constants
and the initial values
yielding a rotating, first contracting and then again expanding gas ball. The initial
positions q i (0), the masses m i and the matrices H i (0) have been determined as described
in the introduction to this section. To reproduce the initial velocity field (6.10)
exactly, the initial velocities q 0
fixed by (6.4), and not via
the linear system (6.2) as in the other examples. The initial entropies are S i
We set R = 250 in this example.
We solved the problem over the period 0  t  100 with the stepsize
in the time-stepping procedure from x4, starting from a grid of gridsize
21 \Theta 21 B-splines finally yielding 325 particles. Fig. 1 shows the initial configuration
of the particles at time and the final configuration at time 100, the latter
rescaled from the radius determined by (6.13) and (6.15) to the initial radius. The
Fig. 1. The particles in the gas cloud example at times
initial configuration of the particles is almost retained over this long period, although
the gas ball first shrinks to the radius 0:515 and then again expands to the final radius
516:9, that is by more than the factor thousand. The exact and the approximate
solution cannot be distinguished with the naked eye; the maximum norm of the error
is about one per thousand of the maximum norm of the solution. Fig. 2 shows a cross
section of the mass density together with the differences to the approximate mass
26 CHRISTOPH GAUGER, PETER LEINEN, HARRY YSERENTANT
densities along the line x
are again rescaled to the same radius and the errors have been multiplied by the
factor thousand in order to compare them with the exact solution. Astonishingly, the
Fig. 2. Mass density and absolute error in the gas cloud example at times
relative size of the error does not increase. Probably one could run this example to
infinity.
There is a viscous counterpart of the solutions above, with the same density profile
and the viscosity coefficients j and i constant multiples of 'ae. Also these solutions are
well reproduced, even if with a slightly reduced accuracy. In examples like these, the
shear viscosity must counterbalance the bulk viscosity; otherwise physical instabilities
occur.
6.2. Shock fronts. The second example shows how the particle method behaves
when shocks develop in an inviscid fluid and demonstrates how the frictional forces
(2.41), (2.42) work. In this example, the transformation of kinetic energy into internal
energy plays a dominant role.
We consider a spherically symmetric shock front ct in an ideal gas (2.30)
and make for the velocity, the density, and the pressure the ansatz
r
x
r
r
r
for ct and
with constant values ae 1 and  1 for jxj ! ct.
To determine the values c, ae 1 and  1 and the functions V (), M () and P () for
1=c from the initial state
and
we first observe that the differential equations (6.5), (6.6) and (6.8) transfer to
where d is again the space dimension. For d ? 1, these differential equations can be
used to compute V , M and P numerically. For
The next step is to express c, ae 1 and  1 as functions of the values
fixed for the time being. Let u c be the radial components of
the velocity on both sides of the shock front relative to the velocity of the front itself.
As continuum mechanics teaches (see [1], for example), the mass flux, the momentum
flux and the energy flux across the shock front are continuous. This is equivalent to
the relations
With help of (6.24) and (6.25), the equation (6.26) can be replaced by the Hugoniot
equation
that equivalently expresses the conservation of energy.
For an ideal gas (2.30), the pressure and the internal energy are coupled by the
equation of state
The conservation of mass (6.24) leads to
c
and with (6.28), the Hugoniot equation (6.27) takes the form
Inserting (6.29) and (6.30) into (6.25), one finally obtains a quadratic equation for
the velocity of the shock front. Only the positive solution
of this equation is of interest; the negative solution leads to a rarefraction shock that
is not compatible with the second law of thermodynamics.
Taking into account the functional dependence (6.23), (6.31) represents an equation
for the unknown velocity c and the corresponding value
This equation can be solved numerically. With (6.29) and (6.30), one then yields also
the values ae 1 and  1 for the density and the pressure inside the shock. As an example,
with
28 CHRISTOPH GAUGER, PETER LEINEN, HARRY YSERENTANT
one obtains
in two space dimensions. The values (6.32) mean that, at time the fluid collides
in the origin with nearly the double speed of sound, which takes the value p
here.
We set R = 250ae in this example. The idea behind this form of R is that friction
should grow when the density of matter increases, as physical intuition suggests. In
our actual computation we placed 251 \Theta 251 particles on a regular grid covering the
square [\Gamma5; 5] 2 and followed their motion in the time interval 0  t  1. Approximately
10000 of these particles are contained in the region [\Gamma2; 2] 2 of interest here
at time slightly more than 20000 at time 1. With the time stepsize
the exponential integrator needed 955 Krylov steps, at most 5 for each of
the 200 time steps, a typical behavior also observed in other examples.
Fig. 3 shows a cross section of the exact and the approximate solution along the
on the interval \Gamma2  x 1  2 and Fig. 4 the corresponding contour lines of the
approximate density at time 1. Both the position and the height of the shock and
the solution profile outside the shock are very well reproduced and the approximate
solution retains its spherical symmetry. The shock is practically as well resolved as
possible with particles of the given size, and as a closer look at the results of the
single time steps shows, the shock does not smear during the computation. Except
for the shock region itself and a small neighborhood of the origin, the error of density
and pressure is less than one per thousand. The same holds for the velocity field
outside the shock. When passing the shock, the particles almost loose all their kinetic
energy such that the velocity is practically zero inside the shock. The error in the
mass density around the origin is probably due to the discontinuous initial data; the
region over which this error extends seems to shrink proportionally to the particle
size. An interesting observation is that behind the shock the particles again rearrange
to a rectangular grid, as can be seen in Fig. 5.
Fig. 3. Density, pressure, and the velocity in the shock front example at time
Fig. 4. Contour lines of the density in the shock front example at time
Fig. 5. Particle contours in a neighborhood of the shock.
6.3. Shear flows. Our last example serves as a test for the modeling of the
shear viscosity. We keep the mass density and the pressure
constant, neglect the heat production, that is set the right-hand side of (6.7) to zero,
and consider velocity fields of the form
where n is a given unit vector fixing the direction of the flow and
component of x in direction of a unit vector e orthogonal to n. Again, the continuity
equation is fulfilled independent of the choice of the functions OE and c. For a
Newtonian fluid (2.55) with constant shear viscosity j, the momentum equation (6.6)
reads
This means that either OE() is linear and c(t) constant, that is
up to a translation in direction of e, or correspondingly
with
The solution corresponding to the velocity field (6.37) is remarkable in so far as
it is reproduced exactly provided that the initial data are exact and the integrals
defining the forces are evaluated exactly. The reason is that all forces acting upon the
particles then vanish such that all particles are distorted in the same way and mass
density and pressure remain constant.
The solution with the velocity field (6.38) represents a more severe test. If one
lets the particles move according to the differential equations (2.21) in the velocity
field (6.38), the particle positions q i and the matrices H i at time t would be
with the function ff(t) given by
Thus, in a little while, the particles have practically lost all their kinetic energy and
get stuck.
In our experiments, we used the data
fixing the physical properties of the fluid, set v as in the
first example. To exclude influences coming from an alignment of the particles with
the flow and to avoid problems with the infinite extension of the flow region, we set
p' 3'
p:
The problem then becomes periodic with period interval [0; 3] \Theta [0; 2]. All computations
were for the time interval 0  t  1 with time stepsize
The particles in Fig. 6 stem from such a computation. They have initially been
located on an axiparallel grid of sidelength Fig. 7 shows how the velocity
Fig. 6. The deformation of particles in the shear flow example.
field and the errors corresponding to the initial gridsizes
evolve in time, where the function values have been sampled on a 500 \Theta 500
grid and the associated l 2 -norm has been taken as distance measure. For better
comparison, the errors have been multiplied by the factors 25,
respectively. The picture demonstrates that the viscous time scale is
Fig. 7. The time evolution of the velocity and the rescaled errors in the shear flow example.
perfectly reproduced. In the transition from 30 \Theta 20 to 60 \Theta 40 and from 60 \Theta 40 to
120 \Theta 80 particles, the error decreases approximately by the factor 16, a clear fourth
order convergence, and a fine confirmation of the finite mass method.



--R

A mathematical introduction to fluid mechanics
New York
Solving ordinary differential equations
On Krylov subspace approximations to the matrix exponential operator
Lehrbuch der Theoretischen Physik
Lehrbuch der Theoretischen Physik


Heuristic numerical work in some problems of hydrodynamics
Proposal and analysis of a new numerical method in the treatment of hydrodynamical shock problems
A particle model of compressible fluids
Particles of variable size
Entropy generation and shock resolution in the particle model of compressible fluids
A convergence analysis for the finite mass method for flows in external force and velocity fields
--TR
finite mass method;gridless discretizations;compressible fluids