1. Introduction

Nowadays, comparative verification of numerical methods is getting more and more in-demand in science and technology. Usually, the comparative verification is conducted on a class of problems that have an exact or numerical solution, acknowledged as reference, or experimental data.

This paper constitutes the part of investigations on comparative verification of numerical methods, held with the help of building up a generalized computational experiment. Generalized computational experiment is computational technology, built on synthesis of mathematical modeling approaches, parallel computations and visual analysis of multi-dimensional data. Methods of making the generalized computational experiment are extensively described in works [1, 2]. Implementation of the generalized computational experiment helps to obtain and process the numerical solutions not only for one particular problem, but also for the class of problems, which is set in the area of space of determining parameters. This turned out to be quite a useful feature in organizing the comparative verification of numerical methods. The investigations dedicated to comparative verification of numerical methods with the help of the generalized computational experiment, are presented in works [3–6].

Presently, the methods of numerical modeling find expanding applications as addition to the experiment and as research techniques replacing the experiment in different academic fields [7].

In solving evolutionary problems the three classes of numerical methods of modeling are widely used: finite-difference methods, finite element methods and particle methods [8].

Particle methods are actively used in modeling the problems of gas dynamics. The model of continuous media is replaced with discrete model, assembly of particles. Each particle has a set of attributes, such as mass, speed, position in space. The condition of the physical system is defined by the set of attributes of a finite number of particles, and evolution of the system is defined by the laws of particles’ interaction.

There are three basic types of particle methods distinguished in scientific literature [9]: particle-particle (PP), particle-mesh (PM) and particle-particle-particle-mesh (P3M).

PP methods use the Lagrangian approach, where the particles move together with the medium. The the PM methods use Eulerian-Lagrangian approach. The computational domain is divided with a fixed mesh (Euler approach), but the particles which move through Eulerian mesh are also taken into consideration (Lagrangian approach). Particles serve for defining liquid’s parameters (mass, energy, speed), and Eulerian mesh is used for defining field parameters (pressure, density, temperature).

In PP methods the force influencing each single particle is calculated by means of summing up the forces from the side of all other particles, in the PM methods the force is a field’s value and is approximated on a mesh. The P3M methods is a hybrid of PP and PM, where for nearby particles (till some preset distance) the force is defined the same way as in PP methods, and for more distant ones, the same way as in PM methods.

In PP methods the condition of the physical system is described by the set of positions and speeds of particles. With transition to a new time layer these values recounted with the use of interaction forces and motion equations. Such methods are computationally expensive, in the context of PM and P3M.

In PM methods the field values, which fill the entire space of the physical system, are approximately presented as values in regularly positioned mesh nodes. As a result, the force is calculated more economically and quickly, but way less precisely rather than with the usage of PP method.

The smoothed particle hydrodynamics (SPH) [10], particle finite element method (PFEM) [11] and discontinuous particle method [12, 13] are referred to as the first type. The second type includes the particle-in-cell method (PIC) [14], large-particle method [15], grid-characteristic methods [16, 17].

As it was shown in articles [12, 18], the discontinuous particle method enables the calculation of discontinuities with high precision, which is extremely important for nonlinear equations in partial derivatives of hyperbolic type.

2. Fundamentals of the Discontinuous Particle Method

Suppose there are material points that are at the initial moment in time in coordinates and moving at speeds . This verbal formulation corresponds to the Cauchy problem:

(1)

The article [10] shows the transition from (1) to the transfer equation in differential form:

(2)

That is, if the coordinates of the points change according to the system of equations (1), then the density is a generalized solution to the Cauchy’s problem for the transport equation (2).

Let's describe the modification of the discontinuous method with a new variant of density correction, previously presented in [19]. The particles selected for correction will be called interacting, and the correction process will be called interaction. Let's introduce a uniform grid in time with a step . We consider the system as a set macroparticles. To describe the particles, we introduce the following notation: is the coordinate of the center of the i-th particle at the k-th moment of time, is the speed of the particle, is the height (density) of the particle. Also each particle has a time invariant mass, which indicates that the method is conservative. If in the past we used the width of the particle (Fig. 1), then the new algorithm is based on the conservation of the mass between the particles.

 

Figure 1 – Approximation of the function by a set of rectangular particles

 

The mass between the coordinates of the particles is equal to the half-sum of the particle masses, and, in the absence of diffusion, it should also remain constant. Let’s introduce a notation is the trace of the mass located between the (i-1)-th and i-th particles. We calculate the trace of mass as the area of trapezoids:

(3)

 

Figure 2 – Introduction of the mass trace

 

Remember the values at the initial point of time . Let's perform the procedure of initializing the parameters of the particles at the initial point in time. Let be given the initial density . The coordinates of the particles can be evenly arranged on the design region, where .

(4)

As shown in [10], the coordinates of the particles in solving the Hopf’s equation must satisfy the system of equations:

(5)

Recall that the algorithm of the particle method is built as a predictor-corrector. First, we solve the system of ordinary differential equations by Euler's explicit method:

(6)

After the particles shift, the distances between them change, which leads to a change in the area of the trapezoid. Therefore, at the stage of the corrector, it is necessary to change the heights of the particles so that the mass between the particles remains constant. Consider the possible cases of particle interaction (Fig. 3):

 

Figure 3 – Particle interaction options

 

1. A particle with a higher density runs over a particle with a lower density, which leads to a decrease in the trapezoidal area between the particles. In this case, in order to preserve the area of the trapezoid, we will increase the height of the particle with a lower density.

2. A particle with a higher density moves away from a particle with a lower density, which results in an increase in the trapezoidal area between the particles. In this case, to preserve the area of the trapezoid, we will reduce the height of the particle with a higher density.

Using these rules of particle rearrangement and selection criteria and as a result of the interaction of the particle with one of the neighbors, the area of the trapezoid with another neighbor for which correction has already been made can change, which indicates the error of the algorithm. The particle interactions that arise in this way are not taken into account.

The corrector changes the height of the i-th particle so that the trapezoidal area between the particles remains constant:

(7)

Therefore, the height of the i-th particle in the new step in time is defined as:

(8)

3. Particle method for gas dynamics equations

The equations of gas dynamics are expressions of the general laws of mass, momentum, and energy. Following [20, 21] let's write down a system of equations for the two-dimensional case in Euler’s variables:

(9)

Perfect gas, . , , , , is density, and -components of velocity, pressure and total energy.

The algorithm for solving a two-dimensional problem is similar to a one-dimensional case. The heights of the particles are found from the initial condition, functions , , , calculated at the centers of particles bases (the points , ).

First, as in the one-dimensional case, we solve the systems of ordinary differential equations for the coordinates of 4 types of particles by the Euler’s method.

In the two-dimensional case, the partner for interaction is chosen by minimizing the "impact" parameter — the angle between the vector of relative velocity and the vector connecting the centers of the particles (algorithmically maximizing the cosine of this angle).

Let be a vector connecting the centers of particles i and j, and be the velocity vector of the particles. Then (Fig. 4):

(10)

 

Figure 4 – Finding the impact parameter

 

Selecting the j-th particle for interaction, we proceed to the one-dimensional problem (Fig.5).

 

Figure 5 – Transition to a one-dimensional problem

 

Next, with the help of a corrector, we change the height of the i-th particle similarly (7) so that the trapezoidal area between the particles remains constant:

(11)

From here we find the preliminary height (so far without taking into account the pressure forces):

(12)

The next step of the algorithm is to take into account the forces of pressure. The difference in pressures to the left and right of the particle leads to a change in the momentum and energy of the particle, that is, to an increase in the volume of the corresponding particles. Similarly [10], we obtain calculation formulas:

(13)

The density, momentum, and energy values calculated in the previous step make it possible to determine the pressure at the center of the particle. To do this, you need to use the equation of state.

To determine the values of the pressure and velocity at the boundary of the particle, a pressure account scheme based on the "interaction" of the particles is used. If at step in time at one of the boundaries of the particle an interaction occurred (in accordance with the criteria described above for a one-dimensional configuration), then the value of the pressure and velocity at this boundary was assumed to be equal to the pressure and velocity of the particle that caused the rearrangement. If no interaction occurred, the pressure at the boundary was supposed to be the same pressure as the center of the particle.

Thus, the volume of , and further increases.

4. Comparison of numerical methods

We will compare the discontinuous particle method (Particles) with the solvers of the open source software package OpenFOAM [22]: rhoCentrlFoam (rCF), pisoCentralFoam (pCF), QGDFoam (QGDF). To do this, by all methods we will solve the classical two-dimensional non-viscous problems in modeling the oblique shock with all methods.

The general flow scheme is presented in Figure 6. A supersonic gas flow with Mach number Ì falls on a flat plate at an angle of β. Before the start of the plate, an oblique shock S occurs [21]. This problem is considered within the system of the Euler equations and has an analytical solution [23].

 

Figure 6 – Flow diagram

 

At the input boundary, the parameters of the unperturbed incoming flow are specified. On the part of the lower boundary corresponding to the flat plate, a non-flow condition is specified. At the output boundary, the boundary conditions are set equal to zero of the derivatives of gas-dynamic functions at the normal to the boundary. At the upper boundary for the velocity components, the boundary conditions are set similarly to the conditions for the input boundary. For the remaining gas-dynamic functions of the upper limit, the boundary conditions are set to zero of the derivatives of the gas-dynamic functions at the normal to the boundary. At the upper boundary for the velocity components, the boundary conditions are set similarly to the conditions for the input boundary. For the other gas-dynamic functions at the upper boundary, the conditions are set similarly to the conditions for the output boundary.

Fig. 7 shows the established solution for the pressure field. Angle of incidence of the incoming flow β = 10°, Mach number M = 2. As a result of the steadying, a qualitative picture of the flow was obtained, corresponding to the analytical solution. It should be noted that the particles are presented in the form of circles for ease of perception.

 

Figure 7 – Pressure distribution

 

Figure 8 shows an animation of how the particle method works. It can be seen that in the region of large gradients, the particle size decreases, which corresponds to the thickening of the grid for grid methods.

 

Figure 8 – Animation of the solution of the problem of supersonic flow around a wedge by the particle method

 

Carefully consider the behavior of gas-dynamic functions in the vicinity of the oblique shock. Figure 9 shows a comparison of all solvers in the form of pressure distribution along the horizontal line AA1 crossing the design area at a distance from the lower boundary of y = 0.15 (see Fig. 6). The exact solution is shown by the dotted line. Numerical methods are indicated by different colors shown in the corresponding table in Fig. 9.

 

Figure 9 – Pressure distribution along the horizontal line

 

The presented figure allows you to judge the degree of smearing of the shock wave front for all the methods considered. The particle method smears the gap into the smallest number of cells, but the shock wave front is shifted towards the area with increased density. Of the remaining methods, the best result is given by the solver rhoCentralFoam. The solver QGDFoam smears the front of the shock wave, and the resulting oscillations are also noticeable at the top of the shock wave.

Also, to assess the deviation of the obtained numerical results from the known exact solution in the entire calculated area, we use an analogue of the norms of L₂:

(14)

Here y_m is the pressure p of the particle m, S_m is the area of the particle. All calculations were carried out when setting the following parameters: flow angle β is 6°, 10°, 15°, 20°, Mach number M_∞ varies from 2 to 3 in increments of 0.5. Thus, the solution was implemented in the region of the space of the determining parameters. In the incoming flow, the following gas-dynamic parameters were set: pressure P_∞ = 101325 Pa, temperature T_∞ = 300 K. The set of calculations with variable parameters is part of the generalized computational experiment described in the works [2–6], where the results of comparative verification of numerical methods implemented in the solvers of the open source software package OpenFOAM [21] were presented. Tables 1–3 show the result of calculating the error rate for pressure.

Table 1 – Deviation from the exact solution, U=2M

Angle	Particles	rCF	pCF	QGDF
6	0.012982	0.013287	0.013744	0.014393
10	0.019326	0.020839	0.021740	0.021850
15	0.029315	0.029893	0.031227	0.028868
20	0.037731	0.038307	0.040417	0.032726

 

Table 2 – Deviation from the exact solution, U=2.5M

Angle	Particles	rCF	pCF	QGDF
6	0.014283	0.015357	0.016346	0.018502
10	0.023192	0.025023	0.026259	0.029376
15	0.035997	0.036192	0.037452	0.042189
20	0.046724	0.045692	0.047210	0.051230

 

Table 3 – Deviation from the exact solution, U=3M

Angle	Particles	rCF	pCF	QGDF
6	0.017582	0.017717	0.018736	0.022639
10	0.030211	0.029721	0.030812	0.037448
15	0.046274	0.043788	0.045160	0.055111
20	0.057836	0.055751	0.057216	0.068286

 

Figure 10 shows the error surfaces in the L₂ norm for all four methods involved in the comparison. Calculations for the QGDFoam solver were made with a selected value of α = 0.1.

It can be seen that the error surfaces of the solvers rhoCentralFoam and pisoCentralFoam are very close, that is, the methods have similar characteristics. At β = 10° and M = 3, the particle method gives less accuracy than the solver rhoCentralFoam. At β = 15°, the particle method gives less accuracy than the solver QGDFoam at M = 2 and less accuracy than the solvers rhoCentralFoam and pisoCentralFoam at M = 3. At β=20°, the particle method gives less accuracy than the solvers rhoCentralFoam and pisoCentralFoam at M = 3. In other cases, the discontinuous particle method is more accurate than the other compared methods.

 

Figure 10 – Error surfaces

 

Fig. 11 shows a close-up of Fig. 10. The error surface of the QGDFoam solver is hidden, as it is very different from other surfaces, which is noticeable in the figure without magnification. As can be seen, the particle method for a large Mach number and angle of attack is only slightly worse in accuracy than the rhoCentralFoam solver.

Thus, the generalized computational experiment suggests that the discontinuous particle method is suitable for solving problems with a strong gradient. Not as high accuracy as in solving one-dimensional problems is due to the fact that although the particle method smears the shock wave front much less (as can be seen from Fig.9), the front is shifted, which negatively affects the accuracy of the method. However, there is a great potential for using the discontinuous particle method to solve practical problems of mathematical modeling.

 

Figure 11 – Close-up of Figure 10

5. Conclusion

With the help of the discontinuous shapeless particle method, the problem of gas dynamics on the formation of an oblique shock is solved. The results obtained were compared with the exact solution in the L₂ norm. Also, the results were processed using scientific visualization tools. It can be seen that the gap is smeared into 2-3 cells relative to the front of the oblique shock, which is less than the smearing of the gap with a solver rhoCentralFoam on comparable grids. Together, this suggests that the discontinuous particle method solves two-dimensional problems of gas dynamics, is especially well suited for solving two-dimensional problems with emerging shock waves. The obtained result is very interesting from the point of view of implementing comparative verification of numerical methods on a reference solution using a generalized computational experiment [1–6].

Historically, the particle method has been used in the problems of finding the interfaces of media, gas-dynamic problems of flowing around bodies, and problems of dynamics of multiphase media. The developed method should increase the computational efficiency of the application of particle methods in these traditional areas, for which research is to be done on comparative verification and evaluation of the effectiveness of the method.

6. Acknowledgements

Calculations were carried out using the hybrid supercomputer K100, installed in the Supercomputer Center for Collective Use of the Keldysh Institute of Applied Mathematics RAS.

References

[1] Bondarev A.E. On the Construction of the Generalized Numerical Experiment in Fluid Dynamics // Mathematica Montisnigri. 2018. Vol. XLII. pp. 52–64.

[2] Bondarev A.E. Multidimensional data analysis in CFD problems / Scientific Visualization. 2014. V. 6, no. 5, pp. 59–66. (in Russian)

[3] Bondarev A.E. On the Estimation of the Accuracy of Numerical Solutions in CFD Problems // ICCS 2019, Lecture Notes in Computer Science (LNCS). 2019. V. 11540, pp. 325–333. doi:10.1007/978-3-030-22750-0_26

[4] On Applying of Generalized Computational Experiment to Numerical Methods Verification / Alekseev A., Bondarev A., Galaktionov V., Kuvshinnikov A., Shapiro L.  // CEUR Workshop Proceedings, 2020, V. 2744, paper 19, Proceedings of the 30th International Conference on Computer Graphics and Machine Vision (GraphiCon 2020), Saint Petersburg, Russia, September 22-25, 2020.

[5] Bondarev A.E. Processing of Visual Results of a Generalized Computational Experiment for the Problem of Supersonic Flow Around a Cone at an Angle of Attack // Scientific Visualization. 2021. V. 13. pp. 104–116. doi:10.26583/sv.13.2.08

[6] Bondarev A.E., Kuvshinnikov A.E. Parametric Study of the Accuracy of OpenFOAM Solvers for the Oblique Shock Wave Problem // IEEE The Proceedings of the 2019 Ivannikov ISPRAS Open Conference (ISPRAS-2019). 2019. pp. 108-112. doi:10.1109/ISPRAS47671.2019.00023

[7] Samarskii A.A., Mikhailov A.P. Principles of Mathematical Modeling. Ideas, Methods, Examples, Taylor and Francis, London and New York, 2002, 349 p. doi:10.1201/9781482288131

[8] Bogomolov S.V. Particle method. Incompressible fluid // Matem. mod. 2003 V. 15. no 1. pp. 46–58. (in Russian)

[9] Hockney R.W., Eastwood J.W. Computer simulation using particles, New York: Taylor&Francis, 1988. doi:10.1201/9780367806934

[10] Liu G.R., Liu M. B. Smoothed Particle Hydrodynamics: A Meshfree Particle Method. World Scientific Publishing, 2003. 472 p. doi:10.1142/5340

[11] The particle finite element method. An overview. / E. Oñate, S.R. Idelsohn, F. Del Pin, R. Aubry // Int. J. Comput. Methods. 2004. V. 1, pp. 267–307. doi:10.1142/S0219876204000204

[12] 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)

[13] 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

[14] 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.

[15] Belotserkovskii O.M., Davydov Yu.M. Large-Particle Method in Gas Dynamics, Nauka, Moscow, 1982. (in Russian)

[16] Magomedov K.M., Kholodov A.S. Grid-characteristic numerical methods, Nauka, Moscow, 1988. (in Russian)

[17] Petrov I.B., Favorskaya A.V., Sannikov A.V., Kvasov I.E. Grid-characteristic method using high-order interpolation on tetrahedral hierarchical meshes with a multiple time step // Math. Models Comput. Simul. 2013. V.5. no.5. pp.409–415. doi:10.1134/S2070048213050104

[18] Bogomolov S.V., Kuznetsov K.V., Particle method for system of gas dynamics equations // Matem. Mod. 1998. V. 10 no.7. pp. 93–100.

[19] Bogomolov. S.V., Kuvshinnikov A.E. A discontinuous shapeless particle method for the quasi-linear transport // Journal of Physics: Conference Series. 2021. Ò. 2099. ¹012009. doi:10.1088/1742-6596/2099/1/012009

[20] 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

[21] Landau L.D., Lifshitz E.M. Fluid Mechanics, Pergamon Press, Oxford, 1987. doi:10.1016/C2013-0-03799-1

[22] OpenFOAM, website: http://www.openfoam.org (accessed May 7, 2022).

[23] Ames Research Staff. Equations, Tables, and Charts for Compressible Flow. NACA TR-II35. 1953.