In
recent years, particle methods have become a valuable tool for numerically
solving partial differential equations and have been successfully applied to a
wide range of problems in astrophysics, plasma physics, solid-state physics,
medical physics, and hydrodynamics [1–4]. In these methods, the solution is
sought as a linear combination of
δ-functions,
where the positions and
coefficients represent the locations and weights of particles, respectively.
The solution is then obtained by tracking the temporal evolution of particle
positions and weights according to the system of ordinary differential
equations derived from considering the weak formulation of the problem. To
recover pointwise values of the computed solution at some time t > 0,
regularization of the particle solution is necessary. Therefore, the efficiency
of the particle method depends on the quality of regularization procedures,
allowing the reconstruction of an approximate solution based on the particle
distribution. Typically, regularization of the particle solution involves
convolution with a so-called smoothing function, which serves as a smooth
approximation of the
δ-function
and,
after appropriate scaling, accounts for the tightness of particle
discretization.
Particle
methods offer numerous advantages over finite-difference methods. The numerical
viscosity introduced by the discretization of convective terms in most
finite-difference methods can significantly degrade the accuracy of the
computational method, especially when a coarse grid is used. Lagrangian-type
methods, on the other hand, can alleviate many issues associated with numerical
viscosity since particles provide a non-dissipative approximation of
convection. Additionally, in scientific applications like kinetic theory,
finite-difference methods may be impractical for realistic scenarios due to
problem dimensionality [5]. In contrast, particle schemes concentrate particles
in the relevant phase space region, optimizing computer memory usage. Being
meshless, particle methods are highly flexible and advantageous for problems
with complex geometry and/or moving boundaries.
However,
it's essential to consider that the self-adaptation of particle positions to
the local flow map comes at the expense of particle distribution regularity.
Distances between particles can vary over time, and particles may cluster near
discontinuities while excessively spreading apart near non-smooth fronts. This
can lead not only to poor resolution of the computed solution but also to
extremely low method efficiency. The latter is related to the fact that the
time step of the ODE solver used to evolve the particle system in time
generally depends on the distance between particles. Therefore, the success of
various particle methods relies not only on the accuracy of reconstruction
procedures, allowing the recovery of pointwise values of the numerical solution
from its particle distribution but also on precise and efficient redistribution
algorithms providing adequate resolution in different regions of the
computational domain.
The
study of applying the discontinuous particle method to viscous problems was
conducted as part of a comprehensive research effort comparing the relative
accuracy of numerical methods on benchmark solutions. This research was
previously carried out for inviscid gas dynamics problems with benchmark
solutions presented in [6–8]. Currently, research is ongoing to comparatively
assess the accuracy of numerical methods for viscous problems using benchmark
solutions. Particularly interesting is the application of the discontinuous
particle method, given its fundamentally different nature compared to common
numerical methods.
Let
there be
material
points, located at the starting moment in coordinates
and moving
with the speed
.
Such verbal formulation corresponds to the
initial value problem (Cauchy problem).
|
|
(1)
|
In
articles [9, 10], the transition from equation (1) to the transport equation in
a differential form is shown:
|
|
(2)
|
As
to say, if the coordinates of points change according to the system of
equations (1), then the density
is a
generalized solution to the Cauchy problem for the transport equation (2).
Let's
describe a modification of the discontinuous method with a different density
correction approach. The particles chosen for correction will be called
interacting, and the correction process will be called interaction. Introduce a
uniform time grid with a step. Consider the system as a set of
macroparticles.
To describe the particles, we introduce the following notations:
— the
coordinate of the center of the
i-th particle at the
k-th moment
in time,
— the velocity
of the particle,
— the height
(density) of the particle.
Also,
each particle has a time-invariant mass, indicating the conservativeness of the
method. The new algorithm is based on the conservation of mass between
particles. The mass between particle coordinates is equal to half the sum of the
masses of the particles and, in the absence of diffusion, should also remain
constant. Let
denote the
mass between the (i-1)-th particle and the i particle. We calculate the mass
as the area of
trapezoids (Figure 1):
Figure 1
– Particles forming a trapezoid
With
the absence of mass diffusion, the mass between particles remains constant over
time. We'll record the values of
at the initial
time moment
.
Let's proceed
with the initialization procedure of particle parameters at the initial time.
Let’s suppose we have an initial density
.
The
coordinates of particles
can be
uniformly distributed over the computational domain, where
.
Coordinates
of particles when solving the Hopf equation must be correct for the system of
equations [11]:
Let's not forget that the particle method algorithm
is constructed as a predictor-corrector. First, we solve the system of ordinary
differential equations using the explicit Euler method:
After
the particle shift, the distance between them changes, leading to changes in
the trapezoid areas. Therefore, at the correction stage, it is necessary to
adjust the heights of the particles so that the mass between the particles
remains constant. Let's consider possible cases of particle interactions
(Figure 2).
Figure 2
– Examples of
particle interactions
A
particle with a higher density collides with a particle with a lower density,
resulting in a reduction of the trapezoid area between the particles. In this
case, to maintain the trapezoid area, we will increase the height of the
particle with lower density.
A
particle with higher density moves away from a particle with lower density,
leading to an increase in the trapezoid area between the particles. In this
case, to maintain the trapezoid area, we will decrease the height of the
particle with higher density.
By
applying these rules for particle rearrangement and selection criteria, the
interaction result with one of the neighbors may change the trapezoid area with
another neighbor, for which correction has already been made, indicating the
algorithmic error. Interactions between particles arising in this way are not
considered.
The
corrector adjusts the height of the
i-th
(
particle in
such a way that the trapezoid area between the particles remains constant:
|
|
(3)
|
Consequently, the height of the
i
particle on the new moment of time
is calculated
this way:
Let’s
analyze Burgers’ equation:
where
— diffusion coefficient.
This equation describes quasi-linear advection with diffusion. In the initial
stage, we solve the advection equation without considering diffusion:
where
— mass,
located between the (i-1)-th and
i
particles on the
temporal
layer,
is the
intermediate value of the particle height. The particle coordinates form a
non-uniform grid. Let's introduce the following notations:
Let's
write out the finite difference approximation of the second derivative on a
non-uniform grid and obtain the value of
on the new
temporal layer:
|
|
(4)
|
Also,
due to diffusion, the mass between particles will change, so it's necessary to
find the new trapezoid areas. Let's consider the
i
trapezoid between the
(i-1)-th and
i-th particles. We'll calculate the flux density at
the trapezoid boundaries according to Fick's law:
The
flux density at the right boundary of the
i-th trapezoid is determined
by the expression:
In
the same way, for the flux density at the left boundary of the
i
trapezoid, what we have is:
The
mass that has flowed from the
of the (i+1)
trapezoid to the
i
is equal to
over the time
step. The mass that has flowed from the (i-1) to the
i
is equal
to
.
The new mass of the
trapezoid is now:
|
|
(5)
|
Among
Burgers’ equation solutions, there is a smoothed wave:
where
,
.
Let it be so that
,
,
,
.
Let’s take
as the initial
condition The comparison between the numerical and analytical solutions is
depicted in Figure 3. The animation starts at time
and ends at
.
Figure
3
–
Applying
particle metho
d to Burgers’ equation
The black dots are the centers of the particles,
the red line is the exact solution. For better visual perception of the
particle density, vertical lines are drawn from the particle centers. It can be
seen that the particle method allows solving such problems with a given accuracy.
It is also worth noting the densification of the particles, due to which the
time step is reduced. In problems with strong gradient densification helps to
better calculate regions with strong gradient, however, in the Burgers’
equation such a high density of particles reduces the space step too much.
Thus, the authors see the necessity of using the "birth-death" of
particles scheme previously applied in [11].
Let us consider the following microscopic Cauchy problem for a
2D case:
which solves the inviscid Burgers’
equation in a weak sense.
We have introduced two sets of particles with the coordinates
,
and
,
for the
functions
and
,
respectively.
The coordinates of the particles
,
and
,
will match each
other if the same number of particles is selected for representing of
u
and
v:
.
Thus, we will
solve the system:
In the two-dimensional case, for selecting an interacting neighbor
one needs to use the aiming parameter, which is equal to the cosine of the
angle between the velocity vector of the particle and the vector connecting the
centers of the particles.
For example,
,
are the coordinates
of the first interacting particle,
,
are the coordinates
of the second particle,
,
are the coordinates
of the first particle’s velocity vector. If we use the definition of an inner
product, then:
Let’s
analyze the test problem for the two-dimensional linear advection equation
[12]:
where
,
.
With such a velocity problem, the initial
profile
should rotate
around the point
without
changing its shape. Let's compare two numerical solutions obtained by the
particle method we described, in the predictor stage, we use the first-order Euler
approximation and the modified second-order Euler method.
Let’s
use the modified second-order Euler method. First, take a time step of
In
the second step, we obtain the coordinates of the particle center
,
at the new
time step using the values
,
:
Let's
take the initial function
in the form of
a right circular cone with a radius of 5 and a height of 1, whose center is
located at the point with coordinates (25, 10):
Figure 4 shows the initial condition.
Figures 5 and 6 show solutions by two methods after three circles.
On the left is a three-dimensional image, and on the right is a top-down
view, more suitable for visual analytics. It is evident that the solution
obtained with the modified method better preserves the shape of the initial
profile.
Figure 4
– Initial conditions
Figure 5
– Numerical solution with particle method using Euler method
Figure 6
– Numerical solution with particle method using
modified Euler method
The
equations of gas dynamics express the general laws of mass, momentum, and
energy. Following [13, 14], let's write the system of equations for the
two-dimensional case in Eulerian variables:
Ideal
gas,
.
,
,
,
,
,
— density,
and
components of
velocity, pressure, total energy, and stress tensor.
The
algorithm for solving the two-dimensional problem is similar to the
one-dimensional case. The computational domain is divided into a finite number
of regions serving as the bases of particles. Particle heights are determined
from the initial condition – functions
,
,
,
computed at
the particle base centers (points
,
).
The radius
of particle bases of all types is assumed to be the same.
First
of all, just like in the one-dimensional case, we solve systems of ordinary
differential equations for the coordinates of the four types of particles using
the Euler method. In the two-dimensional case, a partner for interaction is
chosen by minimizing the "aiming" parameter—the angle between the
relative velocity vector and the vector connecting the centers of the particles
(algorithmically, maximizing the cosine of this angle). Once we've selected the
j-th particle for interaction, we move on to a one-dimensional problem.
Then, using a corrector, we adjust the height of the i particle similarly to
(3) to keep the area of the trapezoid between the particles constant.
From
here, we can find the preliminary height (without considering pressure forces
yet):
The
next step in the algorithm involves accounting for pressure forces. The
pressure difference on the left and right sides of a particle leads to changes
in its momentum and energy, increasing the volume of the respective particles.
Similarly to [8], we arrive at the following computational formulas:
The
computed values of density, momentum, and energy from the previous step allow
for determining the pressure at the particle's center. To
achieve
this,
the
equation
of
state
is
applied.
To
determine the pressure and velocity values at the particle boundary, a pressure
calculation scheme based on particle "interaction" is employed to
determine the pressure and velocity values at the particle boundary. If an
interaction occurred at the particle boundary during the time step (according
to the criteria described above for the one-dimensional configuration), then
the pressure and velocity at that boundary are set equal to the pressure and
velocity of the particle that caused the rearrangement. If no interaction
occurs, the pressure at the boundary is set equal to the pressure at the
particle's center. Consequently, the volumes of particles
,
,
and
are
additionally increased.
Next,
we need to consider the diffusive part. It is needed to take a particle
and select all
neighboring particles
within a
certain interaction radius R. For each pair
(
),
simplify the
problem to make it one-dimensional and use formulas (4) and (5). It's worth
stating that the current algorithm requires finding neighbors at each time
step, significantly impacting the method's performance and eliminating one of
the advantages compared to traditional numerical methods where grid nodes are
stationary, and neighbors are always known.
As
a test case to develop the capabilities of the particle method, the problem of
a compressible gas flow around a plate is considered. The boundary conditions
on the plate are standard no-slip conditions (u=v=0). A significant difference
from the Blasius problem is that a singularity forms at the plate's nose. For
the boundary condition, special particles with zero velocity are used; their
height does not change.
For
the laminar boundary layer, the layer thickness
for the
incompressible case is determined by the relation [15]:
With
this, it is possible to calculate the thickness of the boundary layer for
compressible gas
:
Figure
7 shows the balanced velocity field. The Mach number
,
free stream temperature
,
free stream
density
,
thickness for compressible gas
.
Figure 7
– Velocity
distribution in the Blasius problem
As
a result of the establishment, a qualitative flow pattern has been obtained.
However, there are visible white spots in the figure, indicating the absence of
particles at the spout. This result can be explained by the fact that at the
nose of the plate the particles experience sharp braking and fly up, thereby
forming a narrow strip without particles. This can also be eliminated by the
mechanism of “birth-death” of particles.
With
the use of the discontinuous particle method without shape consideration two diffusion
problems were solved: the one-dimensional problem for Burgers’ equation and
Blasius’ two-dimensional problem. It was also discovered that the use of a
high-precision circuit at the corrector stage increases the final accuracy of
the particle method. However, in problems with viscosity, regions arise with an
increased number of particles, which is why the time step has to be reduced,
and regions in which there are no particles, that is, the numerical result is
not defined. It is proposed to solve such problems by the mechanism of
“birth-death” of particles, which will be done in subsequent works. The
developed approach ensures the participation of the discontinuous particle
method in the general comparative analysis of the accuracy of numerical methods
on reference solutions for viscous problems.
Calculations
were performed on the hybrid supercomputer K-100 installed in the Supercomputer
Centre of Collective Usage of KIAM RAS.
[1]
Harlow F.H. The
Particle-in-Cell Computing Method for Fluid Dynamics // Methods in
Computational Physics. Vol. 3, ed. B Alder, S Fernbach, M Rotenberg. New York:
Academic Press. 1964. pp.
319–343.
[2]
Hockney R.W., Eastwood
J.W. Computer simulation using particles, New York: Taylor&Francis, 1988.
doi:10.1201/9780367806934
[3]
Choquin J.P., Huberson S. Particle
simulation of viscous flow // Comput. & Fluids. 1989 Vol. 17 no.2. pp. 397–410.
doi:10.1016/0045-7930(89)90049-2
[4]
Hermeline F. A
deterministic particle method for transport diffusion equations: Application to
the Fokker-Planck equation // J. Comput. Phys. 1989. Vol. 82. no.1. pp. 122–146.
doi:10.1016/0021-9991(89)90038-7
[5]
Dimarco G., Pareschi L. Numerical
methods for kinetic equations // Acta Numerica. 2014. Vol. 23. pp. 369–520.
doi:10.1017/S0962492914000063
[6]
Alekseev A.K., Bondarev
A.E., Kuvshinnikov A.E. Comparative analysis of the accuracy of OpenFOAM
solvers for the oblique shock wave problem // Mathematica Montisnigri. 2019.
Vol. XLV. pp. 95–105. doi:10.20948/mathmontis-2019-45-8
[7]
Alekseev A.K.,
Bondarev A.E., Kuvshinnikov A.E. On uncertainty quantification via the ensemble
of independent numerical solutions // Journal of Computational Science 2020. Vol. 42.
101114. doi:10.1016/j.jocs.2020.101114
[8]
Bondarev A.E, Kuvshinnikov
A.E. Analysis and Visualization of the Computational Experiments Results on the
Comparative Assessment of OpenFOAM Solvers Accuracy for a Rarefaction Wave
Problem // Scientific Visualization. 2021. Vol. 13. no. 3. pp. 34–46.
doi:10.26583/sv.13.3.04
[9]
Bogomolov S.V.,
Kuvshinnikov A.E. Discontinuous particles method on gas dynamic examples // Math.
Models Comput. Simul
.
2019 Vol.11 no.5. pp. 768–777. doi:10.1134/S2070048219050053
[10]
Bogomolov
S.V, Bondarev A.E., Kuvshinnikov A.E. Comparative Verification of Numerical
Methods Involving the Discontinuous Shapeless Particle Method// Scientific
Visualization. 2022. Vol. 14. no.4. pp. 97–109. doi:10.26583/sv.14.4.09
[11]
Bogomolov
S.V., Zvenkov D.S. Explicit particle method, non-smoothing gas-dynamic discontinuities.
Matem. Mod. 2007. Vol. 19. no.3. pp. 74–86. (in Russian)
[12]
Crowley
W.P. Numerical advection experiments // Monthly Weather Rev. 1968. Vol. 96.
no.1. pp. 1–11. doi:10.1175/1520-0493(1968)096<0001:NAE>2.0.CO;2
[13]
Rozdestvenskii
B.L., Janenko N.N. Systems of Quasilinear Equations and Their Applications to
Gas Dynamics,
American Mathematical Society, 1983. doi:
10.1090/mmono/055
[14]
Landau
L.D., Lifshitz E.M. Fluid Mechanics, Pergamon Press, Oxford, 1987. doi:10.1016/C2013-0-03799-1
[15]
Schlichting
H., Gersten K. Boundary layer theory, 9th ed. Berlin Heidelberg, Springer-Verlag,
2016. doi: 10.1007/978-3-662-52919-5
