**10. Method testing**

#### **10.1 Test problem 1. A moving cruciform density discontinuity**

Domain 0 < x < 12, 0 < у < 12 is divided into two subdomains (0 and 1). In subdomain 0: ρ<sup>0</sup> ¼ 1, e0 =0, ux ¼ 1, uy ¼ 1, no particles are specified; in subdomain 1: ρ<sup>0</sup> ¼ 10, e0 =0, ux ¼ 1, uy ¼ 1, each cell contains one particle. Р = 0 all over the domain, so the problem involves virtually no gas dynamics, only convective flow. The calculations were performed on a fixed grid of 60x60 cells.

Results calculated by the SP method and EGAK (in what follows, SP is the particle method with flux correction described above, PIC is the correction-free particle method similar to the PIC method, and EGAK is used to denote the particlefree code) are shown in the form of density distributions at t = 7.08 (**Figure 7**). The figure shows that the result produced by the SP method is exact.

Calculations were performed by the SP, PIC and EGAK methods. In the SP and PIC calculations, each cell initially contained four particles. In the SP problem calculation, the number of particles in the region with gas increased because of its strong expansion. The calculations were done on a fixed regular grid of 100х100 cells. The SP calculation results are presented in **Figure 9** as a density distribution at t = 100. **Figure 10** shows plots of density as a function of radius for all the cells in the domain occupied by the shock wave at a given instant. They demonstrate how efficient the methods are in preserving the flow's spherical symmetry. **Figure 11**

shows plots of density as a function of radius for section x = y.

*Test problem 3. Density across the cells as a function of radius, t = 100: a) EGAK; b) SP; c) PIC.*

*Test problem 3. Density as a function of radius in section x = y, t = 100.*

**Figure 9.**

**Figure 10.**

**Figure 11.**

**55**

*Test problem 3. Density distribution at t = 100.*

*A Monotonic Method of Split Particles DOI: http://dx.doi.org/10.5772/intechopen.97044*

### **10.2 Test problem 2. A one-dimensional steady shock wave**

We consider the following 1D problem statement. Domain 0 < x < 50, 0 < y < 4 is occupied by an ideal gas with ρ = 1, Р = 0, u = 0, γ = 3. A plane shock wave propagates in the material from left to right. Its parameters behind the shock front are ρ = 2, e = 2, Р = 8, u = 2. The calculations were performed on a fixed grid of 100x4 cells. In the PIC and SP calculations, there were four particles in each cell.

**Figure 8** shows the plots of density as a function of coordinate at t = 10 calculated by the SP, PIC and EGAK methods. One can see that the SP method gives nearly the same result as EGAK and a better result compared to PIC. The exact shock position is х = 40.

#### **10.3 Test problem 3. A point explosion**

Domain 0 < x < 20, 0 < y < 20 contains two materials: a circle of radius 0.1 with its center at the origin is occupied by an ideal gas with ρ = 1, е = 1, Р = 0, γ = 1.4; the remaining part is occupied by an ideal gas with ρ = 1, е = 0, Р = 0, γ = 1.4.

**Figure 7.** *Test problem 1. Density distributions: a) t = 0; b) t = 7.08, method SP; c) t = 7.08, method EGAK.*

**Figure 8.** *Test problem 2. Density as a function of shock position, t = 10.*

### *A Monotonic Method of Split Particles DOI: http://dx.doi.org/10.5772/intechopen.97044*

Results calculated by the SP method and EGAK (in what follows, SP is the particle method with flux correction described above, PIC is the correction-free particle method similar to the PIC method, and EGAK is used to denote the particlefree code) are shown in the form of density distributions at t = 7.08 (**Figure 7**). The

We consider the following 1D problem statement. Domain 0 < x < 50, 0 < y < 4 is occupied by an ideal gas with ρ = 1, Р = 0, u = 0, γ = 3. A plane shock wave propagates in the material from left to right. Its parameters behind the shock front are ρ = 2, e = 2, Р = 8, u = 2. The calculations were performed on a fixed grid of 100x4 cells. In the PIC and SP calculations, there were four particles in each cell. **Figure 8** shows the plots of density as a function of coordinate at t = 10 calculated by the SP, PIC and EGAK methods. One can see that the SP method gives nearly the same result as EGAK and a better result compared to PIC. The exact

Domain 0 < x < 20, 0 < y < 20 contains two materials: a circle of radius 0.1 with its center at the origin is occupied by an ideal gas with ρ = 1, е = 1, Р = 0, γ = 1.4; the remaining part is occupied by an ideal gas with ρ = 1, е = 0, Р = 0, γ = 1.4.

*Test problem 1. Density distributions: a) t = 0; b) t = 7.08, method SP; c) t = 7.08, method EGAK.*

figure shows that the result produced by the SP method is exact.

**10.2 Test problem 2. A one-dimensional steady shock wave**

shock position is х = 40.

**Figure 7.**

**Figure 8.**

**54**

*Test problem 2. Density as a function of shock position, t = 10.*

**10.3 Test problem 3. A point explosion**

*Recent Advances in Numerical Simulations*

Calculations were performed by the SP, PIC and EGAK methods. In the SP and PIC calculations, each cell initially contained four particles. In the SP problem calculation, the number of particles in the region with gas increased because of its strong expansion. The calculations were done on a fixed regular grid of 100х100 cells.

The SP calculation results are presented in **Figure 9** as a density distribution at t = 100. **Figure 10** shows plots of density as a function of radius for all the cells in the domain occupied by the shock wave at a given instant. They demonstrate how efficient the methods are in preserving the flow's spherical symmetry. **Figure 11** shows plots of density as a function of radius for section x = y.

**Figure 9.** *Test problem 3. Density distribution at t = 100.*

**Figure 10.**

*Test problem 3. Density across the cells as a function of radius, t = 100: a) EGAK; b) SP; c) PIC.*

**Figure 11.** *Test problem 3. Density as a function of radius in section x = y, t = 100.*

These figures demonstrate that the SP result is again nearly as good at the EGAK result and noticeably more accurate than the PIC result.

## **10.4 Test problem 4. A spherically converging shell**

The initial problem geometry is borrowed from paper [16] and shown in **Figure 12**. Domain 1: R1 = 0.8, ρ<sup>0</sup> = 0.01, e0 = 0, u0 = 0, ideal gas with γ ¼ 5*=*3. Domain 2: R2 = 1, ρ<sup>0</sup> = 10, e0 = 0, U<sup>R</sup> <sup>0</sup> = � 1, equation of state of Mie-Grueneisen type with constants ρ<sup>0</sup> = 10, c0 *=* 4, n *=* 5, γ *=* 2. Boundary condition at R2: pressure P = 0. Domain 3: vacuum with P = 0. Units of measurement: ρ - [g/cm<sup>3</sup> ], t - [10 μs], L - [cm]. The calculations were done on a fixed regular grid of 110x110 cells. The SP calculations involved particles in both gas and shell domains (one particle per cell, with a limitation of no more than four particles). Their results are compared with the EGAK results. An interesting feature of this problem is that the number of particles in the computation constantly decreases, because both materials are compressed.

The main target result in this problem is maximum gas compression. Note that the reference density obtained in convergence calculations by the 1D method [17] is ≈25. The maximum average gas density and respective time for SP and EGAK are 16.03 at t = 0.368 and 16.49 at t = 0.369, respectively. As an illustration, **Figure 13a** shows a fragment of the domain with particles at t = 0.368. The figure also shows density values across the cells of the compressed-gas domain from the calculations by SP (**Figure 13b**) and EGAK (**Figure 13c**).

The maximum gas compression ratios and their respective time obtained by EGAK and SP are close, but the maximum compression ratios are much lower than the reference solution, which is explained by a small number of cells in these calculations (the solutions converge to the reference solution with mesh size). A smaller compression ratio in the SP calculation compared to EGAK is attributed to the presence of gas spots "split-off" from the main domain in the SP calculation. As for the flow symmetry preservation, SP is nearly as good as EGAK.

**11. Conclusion**

**Figure 13.**

small number of particles.

the EGAK code.

**57**

**Acknowledgements**

The paper describes a monotonic split-particle method. The method has been developed to simulate multi-material gas dynamic flows using a combination of grid methods implemented in EGAK and the SP method for some layers. The calculations demonstrated that the SP method is close to EGAK in the accuracy of shockcapturing simulations and is much more accurate as applied to convective flow simulations, like the PIC method. At the same time, the SP method is free of the major drawback of the PIC method, namely the severe nonmonotonicity of its solution due to the discrete mass transfer. In addition, this method uses a relatively

*Test problem 4. a) Particle positions in the central domain; b) Density value across the gas cells calculated by*

*SP, and c) Density value across the gas cells calculated by EGAK.*

*A Monotonic Method of Split Particles DOI: http://dx.doi.org/10.5772/intechopen.97044*

Further prospects of the SP method are related to its application to the problems that require "remembering" the process history at Lagrangian points, like detonation and combustion of explosives, elastoplastic behavior and fracture of materials, etc. In particular, implemented to date have been the kinetic model of explosive HE transformation by Morozov and Karpenko [18] and the model of materials fracture by Kanel et al. [19]. This method has also been implemented in the 3D extension of

We would like to thank Tatiana Zezyulina for the translation of this paper.

**Figure 12.** *Geometry of test problem 4.*

*A Monotonic Method of Split Particles DOI: http://dx.doi.org/10.5772/intechopen.97044*

**Figure 13.**

These figures demonstrate that the SP result is again nearly as good at the EGAK

<sup>0</sup> = � 1, equation of state of Mie-Grueneisen type

], t - [10 μs], L -

The initial problem geometry is borrowed from paper [16] and shown in **Figure 12**. Domain 1: R1 = 0.8, ρ<sup>0</sup> = 0.01, e0 = 0, u0 = 0, ideal gas with γ ¼ 5*=*3.

with constants ρ<sup>0</sup> = 10, c0 *=* 4, n *=* 5, γ *=* 2. Boundary condition at R2: pressure P = 0.

The main target result in this problem is maximum gas compression. Note that the reference density obtained in convergence calculations by the 1D method [17] is ≈25. The maximum average gas density and respective time for SP and EGAK are 16.03 at t = 0.368 and 16.49 at t = 0.369, respectively. As an illustration, **Figure 13a** shows a fragment of the domain with particles at t = 0.368. The figure also shows density values across the cells of the compressed-gas domain from the calculations

The maximum gas compression ratios and their respective time obtained by EGAK and SP are close, but the maximum compression ratios are much lower than the reference solution, which is explained by a small number of cells in these calculations (the solutions converge to the reference solution with mesh size). A smaller compression ratio in the SP calculation compared to EGAK is attributed to the presence of gas spots "split-off" from the main domain in the SP calculation. As for the flow symmetry preservation, SP is nearly as good as

[cm]. The calculations were done on a fixed regular grid of 110x110 cells. The SP calculations involved particles in both gas and shell domains (one particle per cell, with a limitation of no more than four particles). Their results are compared with the EGAK results. An interesting feature of this problem is that the number of particles in the computation constantly decreases, because both materials are

result and noticeably more accurate than the PIC result.

**10.4 Test problem 4. A spherically converging shell**

Domain 3: vacuum with P = 0. Units of measurement: ρ - [g/cm<sup>3</sup>

Domain 2: R2 = 1, ρ<sup>0</sup> = 10, e0 = 0, U<sup>R</sup>

*Recent Advances in Numerical Simulations*

by SP (**Figure 13b**) and EGAK (**Figure 13c**).

compressed.

EGAK.

**Figure 12.**

**56**

*Geometry of test problem 4.*

*Test problem 4. a) Particle positions in the central domain; b) Density value across the gas cells calculated by SP, and c) Density value across the gas cells calculated by EGAK.*

### **11. Conclusion**

The paper describes a monotonic split-particle method. The method has been developed to simulate multi-material gas dynamic flows using a combination of grid methods implemented in EGAK and the SP method for some layers. The calculations demonstrated that the SP method is close to EGAK in the accuracy of shockcapturing simulations and is much more accurate as applied to convective flow simulations, like the PIC method. At the same time, the SP method is free of the major drawback of the PIC method, namely the severe nonmonotonicity of its solution due to the discrete mass transfer. In addition, this method uses a relatively small number of particles.

Further prospects of the SP method are related to its application to the problems that require "remembering" the process history at Lagrangian points, like detonation and combustion of explosives, elastoplastic behavior and fracture of materials, etc. In particular, implemented to date have been the kinetic model of explosive HE transformation by Morozov and Karpenko [18] and the model of materials fracture by Kanel et al. [19]. This method has also been implemented in the 3D extension of the EGAK code.

### **Acknowledgements**

We would like to thank Tatiana Zezyulina for the translation of this paper.

*Recent Advances in Numerical Simulations*

**References**

[1] Hirt CW. Nicols BD. Volume of Fluid (VOF) Method for the Dynamics of Free Boundaries. Journal of Computational

[11] Lapenta G, Brackbill JU. Dynamic and selective control of the number of particles in kinetic plasma simulations. Journal of Computional Physics. 1994;

[12] Welch DR, Genoni TC, Clark RE, Rose DV. Adaptive particle managment in

[13] Sapojnikov GA. A combine method of fluid fluxes and particle-in-cell for calculations of the gasdynamical flows. In "Voprocy razrabotki I ekspluatacii paketov prikladnyx program".

Novosibirsk: ITPM SO AN SSSR. 1981;

multidimensional multi-material flow

Transactions. Research Publication, Sarov: RFNC-VNIIEF. 2008; **12**: 54–65.

[15] Yanilkin YuV, Goncharov EA, Kolobyanin VYu, Sadchikov VV, Kamm JR, Shashkov MJ, Rider WJ. Multi-material pressure relaxation

hydrodynamics. Computers & Fluids.

[16] Yanilkin YuV, Toporova OO. Twodimensional scalar artificial viscosity of the EGAK code in spherical systems. VANT, Series MMFP. 2010; 3: 46–54.

[17] Bondarenko YuA. The order of approximation, the order of numerical convergence and cost efficiency of Eulerian multi-dimensional gas dynamics computations illustrated by "blast waves" problem simulations. VANT, Series MMFP. 2004; 4, 51–61.

[14] Yanilkin YuV, Belyaev SP, Bondarenko YuA, et al. EGAK and

TREK: Eulerian codes for

methods for Lagrangian

2013; 83: 137–143.

(In Russian).

(In Russian).

simulations. RFNC-VNIIEF

89–97. (In Russian).

(In Russian).

a Particle-in-Cell Code. Journal of Computional Physics. 2007; 227: 143–155.

115: 213–217.

Physics. 1981; 39: 201–225.

569. (In Russian).

[2] Bakhrakh SM, Glagoleva YuP, Samigulin MS, Frolov VD, Yanenko NN,

*A Monotonic Method of Split Particles DOI: http://dx.doi.org/10.5772/intechopen.97044*

Yanilkin YuV. Gas dynamic flow simulations by the concentration method. DAS USSR. 1981; **257 (3)**: 566–

[3] Harlow FH. The particle-in-cell computing methods for fluid dynamics. Meth. Comput. Phys. 1964; 3: 319–343.

Schneider R, Taccogna F. The Particle-In-Cell Method. Contributions to Plasma Physics. 2007; 47(8–9): 563–594.

[5] Shalaby M, Broderick AE, Chang P, Pfrommer C, Lamberts A, Puchwein E. SHARP: A Spatially Higher-order, Relativistic Particle-in-Cell Code. The Astrophysical Journal. 2017; 841(1): 52.

[6] Jiang C, Schroeder C, Selle A, et al. The affine particle-incell method. ACM

[7] Bogomolov SV, Zvenkov DS. Explicit particle method, nonsmoothing gas-

Matematicheskoe modelirovanie. 2007;

[8] Jiang C, Schroeder C, Teran J. An angular momentum conserving affineparticle-incell method. Journal of Computational Physics. 2017; 338:

[9] Fu C, Guo Q, Gast T, et al. A

Cell" Methods. Theory and

polynomial particle-in-cell method. ACM Trans. Graph. 2017; **36**: 222:1–222:12.

[10] Grigoryev YuN, Vshivkov VA, Fedoruk MP. Numerical "Particle-in-

Applications. Utrecht Boston. 2002.

Trans. Graph. 2015; 34:4.

dynamic discontinuities.

19:3: 74–86. (In Russian).

137–164.

**59**

[4] Tskhakaya D, Matyash K,
