A discontinuous shapeless particle method for the quasi-linear transport

This paper considers a new version of the discontinuous particle method, whose higher accuracy is based on the “predictor-corrector” scheme. The peculiarity of this version is a new criterion of rearranging particles at the “corrector” stage. In contrast to the previously used version with the analysis of overlapping particles, which required an assumption about their form, we use another key characteristic of particles, namely, their mass, more precisely, the assumption that in the nonlinear elastic transport not only particle masses are conserved but also the mass located between the centers of these particles. This requirement leads to the fact that changing a distance between particles in the process of their movement and conservation of mass in the space between them, lead to a change in the density of one of the particles. A new version arose in the solution of the two-dimensional transport problems. We emphasize that the discontinuity is smeared into a single particle, which indicates to a high accuracy of the method. The construction of the method for a simple nonlinear transport problem is a necessary step to simulate the complex gas dynamics problems.


Introduction
The main feature of gas dynamics is the emergence of discontinuities, or more precisely, regions of strong gradients in the flow field due to the nonlinear transport. The quality of computational methods is evaluated, first of all, by their ability to transfer such a behavior of the solution as adequately as possible. The discontinuous particle method [1][2][3] allows one to cope with these difficulties better than the alternative, conventionally widespread classes of difference and finite element methods. This is achieved because the particle method is based on a microscopic view of a continuous medium, or the Lagrangian approach. In addition, the particle methods have a constructive propensity for paralleling, are efficient in terms of multidimensionality, and are ideologically organic with respect to hierarchical transitions between micro -macro models of the phenomena under consideration. The main ideas of constructing a new (discontinuous) version of the method and examples of its application will be presented.
The quasilinear transport equation, the Burgers equation [4,5], viscous and non-viscous, are simple models, which exhibit the key property for the whole gas dynamics: the existence of zones of sharp gradients of the solution. Having correctly numerically simulated such zones, one can expect building effective computational methods, not only and not so much for macro-models, but also for mesostages of the hierarchy in the representation of the Kolmogorov-Fokker-Planck equations [1,6,7]. Hockney and Eastwood [8] divide the whole set of particle methods into 3 subclasses: particle-grid (PM) methods, particle-particle method (PP) and hybrid particle-particle -particle-grid methods (PPPM or P3M). In this classification, the discontinuous particle method can be classified as a particle-particle method. The smoothed particle hydrodynamics method and its numerous modifications [9] or the material point method [10] belong to the same subclass. The particle method is used to solve problems in the plasma physics and gas dynamics [11][12][13]. Despite the extensive literature, numerous applications and the obvious physical origin as a division of the whole mass into separate parts followed by tracking their behavior, in our opinion, the particle method requires a clearer mathematical formalization.

The one-dimensional linear transport
Consider the Cauchy problem for a system of differential equations describing the motion of N material points with the given velocities ( , ) where ( ) is the trajectory of motion of the i-th material point, (0) is its initial position.
Such problems are typical of many natural processes, the diffusion of gases or liquids, or the movement of planets and galaxies. If the number of particles is relatively small, the problem can be solved analytically or by some numerical method. However, some phenomena, such as already mentioned diffusion of gases and liquids, require huge numbers of particles for their accurate description. Typical values of the number of molecules in a gas volume are 10 19 -10 23 . Such quantities cannot be simulated even by the most powerful supercomputers.
The solution is to replace some sets of microparticles with a single macroparticle, thus significantly reducing the number of entities in the model and enabling numerical calculations of the problem.
Let us determine the density of the distribution ( , ): or (as is common in physics): Let us differentiate both parts of (2) with respect to time, use the formula for differentiating a complex function, and, as is shown in [2], move from problem (1) to the generalized transport equation: and then to the transport equation in the differential form: That is, the particle density function is a generalized solution of the Cauchy problem for the onedimensional transport equation.
The basic idea of the numerical particle method is that by solving system (1) for an acceptable (relatively small) number of particles, we obtain an approximate solution (4).

The one-dimensional quasilinear transport
Consider the following ODE system: If in the right-hand side of (1) we take the velocities, which are the density values, we obtain a nonlinear Cauchy problem which has no direct physical sense, but leads, by the above described scheme, to a quasi-linear transport equation. Thereby there appears a relation between trajectories of particles and their density of distribution, significantly complicating the solution of the problem:

The two-dimensional quasilinear transport
Consider the following ODE system:  By analogy with the transition from system (5) to equation (6), we carry out the transition for system (7):

The numerical method (a one-dimensional case)
Let us introduce the following notations: is the coordinate of the center of the i-th particle at the kth moment of time, is the velocity of the particle, ℎ is the height (density) of the particle, is the area (mass) of the particle, is the width of the particle. The particle method algorithm is constructed as a predictor-corrector. First, at the predictor stage, we solve the system of differential equations using the explicit Euler scheme: After occuring the shift, there may appear intersections or gaps among the particles, which represent the approximation error. Let us choose for each particle a partner with the strongest influence on it, or let us leave the particle unchanged. Let us consider possible cases of interaction of particles: 1. The higher-density particle flies away from the lower-density particle. In this case, it is necessary to expand the particle with a higher density, thereby eliminating the space among the particles so that they touch each other. As a particle expands, the density of the particle decreases. From a physical point of view, a denser particle should have a lower density.
2. The particle with a higher density collides with the particle with a lower density. In this case, it is necessary to eliminate the resulting intersection and to reduce the width of the particle with a lower density so that it only touches the neighboring one. Reducing the width of the particle increases its density at the expense of mass conservation. From the physical point of view, the interaction of particles with different densities increases the density of a less dense particle.
Thus, at the corrector stage, we will change the widths and heights of the particles so that the intersections and gaps be minimized. We will keep the particles symmetric with respect to the center.
In the first case, we check that the left edge of the right particle is to the right of the right edge of the left particle and that the height of the right particle is greater than the height of the left one. In the second case, the left edge of the right particle is to the left of the right edge of the left particle and the height of the left particle is greater than the height of the right one. In both cases, the same expression can be applied to determine a relative position of the neighboring particles: The positions of particles after the shift do not change, the width is adjusted to the width of the neighboring particle, and the height is calculated from the requirement that the particle area (mass) is constant. As a result we have: By looping through all particles at each time layer and recalculating their parameters according to the specified rules, we will obtain an implementation of the one-dimensional discontinuous particle method.

Features of the particle method algorithm in the two-dimensional case
The general scheme of the discontinuous particle method in the two-dimensional case repeats the one for the one-dimensional case: first, the prediction by the explicit Euler scheme takes place, then the parameters of the particles are corrected, which move apart or vice versa overlapping each other. However, in the one-dimensional case we had to trace only interactions of the neighboring particles, whereas in the two-dimensional case, the number of the neighbors of each particle can be an arbitrary number, and with each neighbor the particle can cross or move away in different degree.
As studies have shown, taking into account all interactions leads to a low accuracy of the simulated solution, so in this study we will apply the method of selecting a neighbor particle, and in the multidimensional case it will be done by the aiming parameter: for each particle i, such a particle j will be selected that the angle between the vector connecting their centers and the velocity vector be minimal.
By choosing a particle in this way, we have reduced the consideration of the interaction to a pair of particles, i.e., to the one-dimensional case. Let us list the basic configurations in which the correction will take place.
1. A particle with a higher density collides with the particle with a lower density: ℎ > ℎ . In this case, we constrict the denser particle so as to eliminate the intersection.
2. The particle with a higher density flies away from the particle with a lower density, i.e. ℎ > ℎ . In this case we expand the particle to eliminate the gap. The volume of each particle is invariant, the width of particles does not also change, which allows unambiguous calculation of the new height. Thus, if one of these two conditions is fulfilled, it is necessary to change the parameters of the i-th particle according to the following formulae: ̂= .
7. An alternate approach in the "corrector" phase Until now, in our previous papers [1][2][3]14], we based on the fact that for the recalculation of a height of particles, which is an approximation of the sought for function, it was necessary to know the base area and based on this fact and the fact of constancy of the volume of particles, which requires the conservativity of the method, the necessary recalculation could be made. At the same, time the shape of particles determines their interactions with each other, which directly affects the algorithm chosen by researchers. In order to get rid of the constant attachment to the shape and size of the particles, it is necessary to introduce another invariant, which ideally would depend only on the heights of the particles and on their position in space.
Consider the two neighboring particles i and j. Let us construct a trapezium whose bases are the heights of the particles, and the segment connecting their centers is the sideline. Assume the area of this trapezium , to be the new invariant (figure 1). ℎ and ℎ is the heights of particles i and j, d is the length of the segment connecting the centers of the particles. Then it is no longer necessary to store the length and width of the base of the particle, and the formula for calculating the changed height takes the following form (figure 2): Thus, we have rid of the binding to the form of the particle. The new invariant can be interpreted as the conservation of mass in the gap between a pair of particles, which is justified by physical considerations.

Numerical examples
First, let us solve the Cauchy problem for the one-dimensional quasi-linear transport equation (5) with the initial conditions of the "triangle with background" kind:  There is an exact solution for such an initial condition: Figure 4. Numerical solution of the one-dimensional quasilinear transport equation by particle method after discontinuity formation.  (5) with initial conditions (9). The red line shows the exact solution. One can see that the rupture velocity is correctly calculated, the rupture is smeared on one particle. After "hitting" a high particle with a higher velocity on a low particle with a lower velocity, the low particle grows to the height of the high particle. Thus, their speeds are compared, the low particle stops growing. In its turn, it "collides" with the next particle.
Then we numerically solve the Cauchy problem for the two-dimensional quasilinear transport equation (7) with the initial conditions of the form: The discontinuity propagation in the two-dimensional problem is shown in figure 5. For the ease of perception, the particles in figure 5 are depicted as circles. It should be understood that the shape of the particles is not used anywhere in the algorithm itself. The exact solution is indicated as a thick black line. You can see in figure 5 that particles with a higher velocity collide with particles with a lower velocity. Those taper, and a one-particle wide gap is formed. This shows the low approximation viscosity of our method. The particles then continue to collide with each other, so that the gap moves at a velocity consistent with the analytical solution. A visual representation of the displacement process of the discontinuity in the computational domain allows us to evaluate the properties of the numerical method used.

Conclusion
We have shown the capabilities of the new version of the discontinuous particle method on an example of the numerical solution of a system of two-dimensional quasilinear transport equations, confirming its high accuracy, which is especially evident on the examples with discontinuous solutions. In our opinion, the clear formalization of our method is also important.