A Monotonic Method of Split Particles

Abstract

The problem of correct calculation of the motion of a multicomponent (multimaterial) medium is the most serious problem for Lagrangian–Eulerian and Eulerian techniques, especially in multicomponent cells in the vicinity of interfaces. There are two main approaches to solving the advection equation for a multicomponent medium. The first approach is based on the identification of interfaces and determining their position at each time step by the concentration field. In this case, the interface can be explicitly distinguished or reconstructed by the concentration field. The latter algorithm is the basis of widely used methods such as VOF. The second approach involves the use of the particle or marker method. In this case, the material fluxes of substances are determined by the particles with which certain masses of substances bind. Both approaches have their own advantages and drawbacks. The advantages of the particle method consist in the Lagrangian representation of particles and the possibility of” drawbacks. The main disadvantage of the particle method is the strong non-monotonicity of the solution caused by the discrete transfer of mass and mass-related quantities from cell to cell. This paper describes a particle method that is free of this drawback. Monotonization of the particle method is performed by spliting the particles so that the volume of matter flowing out of the cell corresponds to the volume calculated according to standard schemes of Lagrangian–Eulerian and Eulerian methods. In order not to generate an infinite chain of spliting, further split particles are re-united when certain conditions are met. The method is developed for modeling 2D and 3D gas-dynamic flows with accompanying processes, in which it is necessary to preserve the history of the process at Lagrangian points.

1. Introduction

Correct calculations of multi-material flows is the greatest challenge for ALE and Eulerian CFD codes, especially those using mixed cells at interfaces. There are two basic approaches to solving the advection equation for the multi-material case. In the first (grid-based) approach, interfaces are identified, and their position on the grid is tracked at each time step. The interface can be identified both explicitly, or it can be recovered based on the field of volume fractions. The latter algorithm serves as a basis for widely used methods, like the VOF method [1] (concentration method [2]). The second approach involves material particle methods first proposed by Harlow (the PIC method [3]). In this case, material fluxes from cells, including mixed ones, are controlled by particles, to which certain material masses are assigned.

Both approaches have their advantages and drawbacks. The advantages of the particle method consist in the Lagrangian representation of particles and the possibility of assigning material information to them. This minimizes the errors of solving the advection equation by the grid-based Eulerian methods. A number of modifications of the PIC method have been developed to improve its accuracy and extend the range of its applications [4, 5, 6, 7, 8, 9]. An overview of these methods is provided in [10].

The central drawback of the particle method is the highly non-monotonic character of the solution caused by the discrete transfer of mass and mass-related quantities from cell to cell. The corresponding error can be reduced in the most straightforward manner by increasing the number of particles in cells, but such an increase limits the method’s performance, especially in the 3D case. To minimize this drawback, a number of method modifications are employed. In [11, 12], for this purpose, the authors use particles having different masses. This approach, however, does not eliminate the need of involving a large number of particles. In a different approach, particles are used only in a limited part of the integration domain, for example, near interfaces [13]. As a result, only a small number of cells contain large quantities of particles. The rest part of the domain in this case is treated by the grid-based methods. Such a selective use of particles, however, does not eliminate the error of solving the advection equation by the grid-based methods and the need of remembering the history of a particular process in a large volume of the material.

This paper proposes a particle method that minimizes these drawbacks. Monotonization of the particle method is performed by particle splitting, so that the material volume flowing out of the cell corresponds to the volume calculated by schemes based on the grid approach. In order to prevent endless splitting, such split particles are further recombined under certain conditions. This approach allows us to do with a small number of particles in the cell, while delivering a monotonic solution.

2. Problem statement

The split-particle (SP) method has been implemented in a code called EGAK in the 2D approximation. In the source code, the major quantities for numerical solution of the multi-material gas dynamic equations include node-centered velocity vector components u_x and u_y and cell-centered thermodynamic quantities: density ρ_ξ, specific (per unit mass) internal energy e_ξ, and volume fractions β_ξ = V_ξ/V of the constituent materials.

Particles can also be specified for some materials (in a particular case, these can be all materials). Each particle (with index р) has its coordinates in space x_р(t), y_р(t) and velocity vector components u_xр(t), u_yр(t) (these are used in interim calculations in the Lagrangian step). In addition, all particles represent thermodynamic states of the corresponding material (density, specific internal energy, volume): ρ_ξр, e_ξр, V_ξр. Note that densities and volumes of particles can also give us their masses. Also note that in the method proposed particle velocities are obtained by interpolation between nodal velocities rather than “remembered” like in the classical PIC method.

Approximation of the corresponding equations is performed in two steps using a decomposition procedure. The first (Lagrangian) step involves calculations of the gas dynamic equations without convective members, i.e. gas dynamic equations in Lagrangian variables. In the second (Eulerian) step, a new grid is constructed (it generally coincides with the grid of the previous timestep), and the quantities are remapped onto it. As inputs, this step takes the outputs of the Lagrangian step. In turn, this step is divided into two sub-steps, i.e., approximation of the advection equation is done with decomposition in directions.

The difference equations below are presented in as much detail as needed to understand the algorithm of particle introduction; please refer to [14, 15] for a more detailed description of EGAK’s basic difference scheme. In what follows, if no confusion is possible, no subscripts or subscript n are used to denote the outputs of the previous timestep, and subscripts n+1/2 and n + 1 denote the outputs of the Lagrangian and the Eulerian step, respectively.

3. Lagrangian step

3.1 Approximation of the cell-centered quantities

The Lagrangian gas dynamic equations are approximated using EGAK’s standard scheme. As outputs, the Lagrange step delivers updated node-centered velocities, as well as densities, energies and volume fractions of each constituent material. This also applies to cells containing particles.

3.2 Definition of particle-specific quantities

In addition to the material-specific quantities, some particle-specific quantities are also defined for particles in the cells containing particles.

3.3 Updating of particle coordinates

Updated particle coordinates are found in two steps:

Step 1. At the Lagrangian step, particles are assumed to move together with the cell and inside the cell, without crossing its boundaries. The relative change in the particle position in the cell is associated with the difference in divergences (compression ratios) of different materials as a result of employing one closing model or another for the mixed-cell gas dynamic equations. In this study, we use only one assumption that the materials have equal divergences. This means that the sub-cell motion of particles does not change their position relative to the grid nodes.

Particle coordinates, ${\tilde{\mathrm{x}}}_{\mathrm{p}}^{\mathrm{n}+1/2}$, ${\tilde{\mathrm{y}}}_{\mathrm{p}}^{\mathrm{n}+1/2}$ are updated by bilinear interpolation between coordinates of cell nodes x^n+1/2, y^n+1/2, just as at time tⁿ.

Step 2. It is easy to show that the calculations of particle velocities by bilinear interpolation violate the law of conservation of momentum in the particle-containing cell. To ensure its conservation, the calculated particle velocities are corrected as follows:

1. Components of cell momentum are calculated:

   $${\mathrm{P}}_{\mathrm{cx}}=\frac{1}{4}\left({\mathrm{u}}_{\mathrm{x}0}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{x}1}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{x}2}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{x}3}^{\mathrm{n}+1/2}\right)\cdot \mathrm{M},$$E1
   $${\mathrm{P}}_{\mathrm{cy}}=\frac{1}{4}\left({\mathrm{u}}_{\mathrm{y}0}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{y}1}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{y}2}^{\mathrm{n}+1/2}+{\mathrm{u}}_{\mathrm{y}3}^{\mathrm{n}+1/2}\right)\cdot \mathrm{M}.$$

   Here, М is the cell mass and u_xi, u_yi (i = 0, 1, 2, 3) are the velocity vector components at four cell nodes.

2. Components of the total momentum of particles belonging to the cell are calculated:

   $${\mathrm{P}}_{\mathrm{px}}=\sum _{\mathrm{p}}{\tilde{\mathrm{u}}}_{\mathrm{xp}}^{\mathrm{n}+1/2}\cdot {\mathrm{m}}_{\mathrm{p},}{\mathrm{P}}_{\mathrm{py}}=\sum _{\mathrm{p}}{\tilde{\mathrm{u}}}_{\mathrm{yp}}^{\mathrm{n}+1/2}\cdot {\mathrm{m}}_{\mathrm{p}.}$$E2

   Here, the particle velocities are calculated using the particles’ coordinates determined by bilinear interpolation and previous particle coordinates:

   $${\tilde{\mathrm{u}}}_{\mathrm{xp}}^{\mathrm{n}+\frac{1}{2}}=\frac{{\tilde{\mathrm{x}}}_{\mathrm{p}}^{\mathrm{n}+1/2}-{\mathrm{x}}_{\mathrm{p}}^{\mathrm{n}}}{\mathrm{\tau }},{\tilde{\mathrm{u}}}_{\mathrm{yp}}^{\mathrm{n}+\frac{1}{2}}=\frac{{\tilde{\mathrm{y}}}_{\mathrm{p}}^{\mathrm{n}+1/2}-y_{\mathrm{p}}^{\mathrm{n}}}{\mathrm{\tau }},$$E3

   were τ = tⁿ⁺¹- tⁿ.

3. Coefficients ${\lambda }_{\mathrm{x}}={\mathrm{P}}_{\mathrm{cx}}/{\mathrm{P}}_{\mathrm{px}}$, ${\lambda }_{\mathrm{y}}={\mathrm{P}}_{\mathrm{cy}}/{\mathrm{P}}_{\mathrm{py}}$ are calculated.

The particles’ velocities and coordinates are updated using the resulting weights:

$${\mathrm{u}}_{\mathrm{xp}}^{\mathrm{n}+1/2}={\lambda }_{\mathrm{x}}\cdot {\tilde{\mathrm{u}}}_{\mathrm{xp}}^{\mathrm{n}+1/2},{\mathrm{u}}_{\mathrm{yp}}^{\mathrm{n}+1/2}={\lambda }_{\mathrm{y}}\cdot {\tilde{\mathrm{u}}}_{\mathrm{yp}}^{\mathrm{n}+1/2};$$E4

$${\mathrm{x}}_{\mathrm{p}}^{\mathrm{n}+1/2}={\mathrm{x}}_{\mathrm{p}}^{\mathrm{n}}+{\mathrm{u}}_{\mathrm{xp}}^{\mathrm{n}+1/2}\cdot \mathrm{\tau },{\mathrm{y}}_{\mathrm{p}}^{\mathrm{n}+1/2}={\mathrm{y}}_{\mathrm{p}}^{\mathrm{n}}+{\mathrm{u}}_{\mathrm{yp}}^{\mathrm{n}+1/2}\cdot \mathrm{\tau }.$$

3.4 Determination of particle velocity, density and energy

Changes in the relative density and energy of particles of a given constituent material are assumed to be equal to the corresponding relative changes in these quantities calculated for the respective material on average. This gives the following formulas:

$${\mathrm{\rho }}_{\mathrm{\xi p}}^{\mathrm{n}+1/2}={\mathrm{\rho }}_{\mathrm{\xi p}}^{\mathrm{n}}+\left({\mathrm{\rho }}_{\mathrm{\xi }}^{\mathrm{n}+1/2}-{\mathrm{\rho }}_{\mathrm{\xi }}^{\mathrm{n}}\right){\mathrm{\rho }}_{\mathrm{\xi p}}^{\mathrm{n}}/{\mathrm{\rho }}_{\mathrm{\xi }}^{\mathrm{n}},$$E5

$${\mathrm{e}}_{\mathrm{\xi p}}^{\mathrm{n}+1/2}={\mathrm{e}}_{\mathrm{\xi p}}^{\mathrm{n}}+\left({\mathrm{e}}_{\mathrm{\xi }}^{\mathrm{n}+1/2}-{\mathrm{e}}_{\mathrm{\xi }}^{\mathrm{n}}\right),$$E6

$${\mathrm{V}}_{\mathrm{\xi p}}^{\mathrm{n}+1/2}={\mathrm{V}}_{\mathrm{\xi p}}^{\mathrm{n}}\left({\mathrm{V}}_{\mathrm{\xi }}^{\mathrm{n}+1/2}/{\mathrm{V}}_{\mathrm{\xi }}^{\mathrm{n}}\right).$$E7

It is easy to show that, when using (5)–(7), the particles’ total masses will remain unchanged, and the particles’ total internal energies will be equal to the energy calculated for the given material as a whole, i.e. the following relationships hold:

$${\mathrm{\rho }}_{\mathrm{\xi }}^{\mathrm{n}+1/2}{\mathrm{V}}_{\mathrm{\xi }}=\sum _{\mathrm{p}}{\mathrm{\rho }}_{\mathrm{\xi p}}^{\mathrm{n}+1/2}{\mathrm{V}}_{\mathrm{\xi p}},{\mathrm{e}}_{\mathrm{\xi }}^{\mathrm{n}+1/2}{\mathrm{m}}_{\mathrm{\xi }}=\sum _{\mathrm{p}}{\mathrm{e}}_{\mathrm{\xi p}}^{\mathrm{n}+1/2}{\mathrm{m}}_{\mathrm{\xi p}}.$$E8

4. Eulerian step

Major difficulties in implementing the particle method are associated with the Eulerian step, and namely, with calculations of mass and internal energy fluxes from cell to cell. In the PIC method, when a particle migrates to a neighbor cell, its mass and energy “migrate” with it. Because of the discrete (and, accordingly, non-monotonic) character of the transfer of mass and all the quantities defined per unit mass, this method delivers highly non-monotonic quantity profiles. Section 4 provides a detailed description of the monotonization algorithm for the PIC method.

In calculations, it is not always efficient to represent all constituent materials by particles because this requires extra calculations and computer memory. Therefore, it is reasonable to use particles for the constituent materials, for which the errors due to the solution of the advection equation are most essential, for example, for thin layers or materials, which require remembering the history of their Lagrangian particle.

As part of the proposed SP method, the following algorithms have been developed:

1. Interaction of materials described by particles and materials calculated by the code’s standard scheme;

2. Support of existence, creation and removal of particles only in the vicinity of the interface;

3. Particles merging;

4. Remapping of particles density and energy to the cell as a whole.

These algorithms are listed in the order of their execution in the Eulerian step after the monotonization algorithm.

5. Monotonization algorithm for the particle method

5.1 One-dimensional case

Let us discuss the concept of the algorithm as applied to a one-dimensional flow for a single particle migrating from cell to cell (Figures 1 and 2). The figures show two cells containing particles represented by dots and imaginary boundaries of volumes represented by dashed segments. Note that calculations by this technique require only numerical values of the volumes, not their layout.

Figure 1.

Illustration of the reason for the non-monotonic behavior: a) volume flux is smaller than the particle volume; b) volume flux is larger than the particle volume.

Figure 2.

Illustration of the monotonization algorithm in the 1D case: a) δV < 0; b) δV > 0.

The flow is directed from left to right as indicated by velocity vector (Figure 1). Volume ΔV flowing out of the left cell (in what follows we call it the volume flux, darker color) is then equal to the product of the cell’s lateral side length L and quantity S = u·τ:

$$\mathrm{\Delta V}=\mathrm{L}·\mathrm{u}·\mathrm{\tau }.$$E9

The non-monotonic behavior of the classical PIC method stems from the discrepancy between the real volume flux (and, accordingly, the mass flux) calculated by (9) and the volume of the particle crossing the cell side. In one case (Figure 1a), the volume moving from the left cell is smaller than the particle volume, and in the other (Figure 1b), it is larger than the particle volume. In the 2D case and in the case of several migrating particles, the situation remains fundamentally the same.

Let us introduce the following notation:

$$\mathrm{\delta V}=\mathrm{\Delta V}‐{\mathrm{V}}_{\mathrm{p}},$$E10

where ΔV is the volume flux calculated by (4), V_P is the volume of particle p migrating from one cell into another.

Below we explain the monotonization algorithm for both of these cases.

1. Volume flux ΔV is smaller than the migrating particle volume (δV < 0). This case is illustrated in Figure 2a. Left is the state at ${\mathrm{t}}^{\mathrm{n}+1/2}$; right, at ${\mathrm{t}}^{\mathrm{n}+1}$. In this case, the particle migrating from the donor to the acceptor cell splits into two parts, a mother and a daughter particle. The mother particle migrates into the acceptor cell and now has new coordinates corresponding to its velocity and a new volume equal to the volume flux leaving the donor cell ΔV. The daughter particle, whose volume is equal to the difference between the initial particle volume and volume flux ΔV, is placed in the donor cell and acquires coordinates on the cell side. The link between the mother and the daughter particle is indicated by a broken line.

2. Volume flux ΔV is larger than the migrating particle volume (δV > 0). This case is illustrated in Figure 2b. In this case, the missing volume of the migrating particle must be made up by forced transfer of some particles or particle fragments from the donor cell to the acceptor. To be split is the particle lying next to the side of these cells and not yet transferred to the acceptor cell. It produces a daughter particle of volume δV, which migrates into the acceptor cell with donor-acceptor side coordinates.

If more than one cell migrates from cell to cell, formula (10) will take the form of

$$\delta \mathrm{V}=\mathrm{\Delta V}-\sum _{\mathrm{p}}{\mathrm{V}}_{\mathrm{p}},$$E11

where ΔV is the volume flux calculated by (9); V_P is the volume of the particle with index p migrating from cell to cell; summing is performed for all transferred particles.

Here, let us describe the differences from the algorithm described above.

1. The volume flux ΔV is smaller than the migrating particles’ total volume (δV < 0). In this case, to be split are all the migrating particles. The mother particles migrating into the acceptor cell have volumes

   $${\mathrm{V}}_{\mathrm{p}}^{\mathrm{M}}={\mathrm{V}}_{\mathrm{p}}\left(\mathrm{\Delta V}/\sum _{\mathrm{p}}{\mathrm{V}}_{\mathrm{p}}\right).$$E12

   The daughter particles stay in the donor cell and have the following volume each:

   $${\mathrm{V}}_p^D={\mathrm{V}}_p-{\mathrm{V}}_p^M.$$E13

2. Volume flux ΔV is larger than the migrating particles’ total volume (δV > 0). The particle staying next to the interface on the donor side is split to fill the remaining volume δV. Note that if V_p < δV, the mother particle’s volume entirely goes to the daughter particles, and to be split is the next particle from the donor cell. If the donor cell is mixed, and the acceptor cell is pure and filled with the material present in the donor cell, the particles of this material will be split first.

5.2 Two-dimensional extension

Of particular interest in this case is the particle transition to the neighbor cell located diagonally from the donor cell (Figure 3). This case is special, because EGAK solves the advection equation using decomposition in directions, whereas no provision is made for diagonal cell-to-cell fluxes.

Figure 3.

Illustration of volume flux calculations in 2D case.

Consider the particular case depicted in Figure 3. Suppose only one node A moves to the new position B in the Eulerian step. The dashed lines in the figure show the locations of the cell sides, for which the grid node is a common vertex, at tn+1/2. C and D denote the points of intersection of straight lines AG and BE, and AF and BH, respectively.

In EGAK, the cell-to-cell volume fluxes are defined as follows. The volume flux corresponding to triangle ABG (for short, volume in triangle ABG) is transferred from cell (i, k) to cell (i + 1, k), that is:

$$\begin{array}{l}{V}_{i,k}^{n+1}=V_{i,k}^{n+1/2}-V_{\mathit{ABG}},\\ V_{i+1,k}^{n+1}=V_{i+1,k}^{n+1/2}+V_{\mathit{ABG}}.\end{array}$$E14

Accordingly, the following relationships apply to all the cells under consideration including all the fluxes:

$$\begin{array}{l}{V}_{i,k}^{n+1}=V_{i,k}^{n+1/2}-V_{\mathit{ABG}}-V_{\mathit{ABF}}=V_{\mathit{\text{BFKG}}},\\ V_{i+1,k}^{n+1}=V_{i+1,k}^{n+1/2}+V_{\mathit{ABG}}-V_{\mathit{ABE}}=V_{\mathit{\text{BGNE}}},\\ V_{i,k+1}^{n+1}=V_{i,k+1}^{n+1/2}-V_{\mathit{ABH}}+V_{\mathit{ABF}}=V_{\mathit{\text{BFLH}}},\\ V_{i+1,k+1}^{n+1}=V_{i+1,k+1}^{n+1/2}+V_{\mathit{ABE}}+V_{\mathit{ABH}}=V_{\mathit{\text{BHME}}}.\end{array}$$E15

Note that in accordance with (15) the volume of triangle ABC is included in the volume of cell (i + 1,k) twice – as part of triangles ABG and ABE – but in one case it is positive, and in the other, negative. Thus, in fact it is not included in the updated volume of this cell; but it will be included in the volume of cell (i + 1,k + 1). The same applies to the volume of triangle ABD, which will be included in the volume of cell (i + 1,k + 1) and not included in the volume of cell (i,k + 1).

In accordance with the above, when considering particle contributions, flux calculations between cell (i,k) and its non-diagonal neighbors assume that the particle lying in triangle ABC migrates into cell (i + 1,k), and the particle lying in triangle ABD, into cell (i,k + 1). Then, in calculations of the flux between cells (i + 1,k), (i + 1,k + 1) and (i,k + 1), (i + 1,k + 1), these particles migrate into cell (i + 1,k + 1). Therefore, when considering this process in terms of flux monotonicity, corresponding daughter particles are introduced as shown in Figure 4. Note that the mother particle’s position during the flux calculations is nevertheless defined in the true acceptor cell (i + 1,k + 1).

Figure 4.

Volume flux from a mixed cell.

The particle splitting is based on the following principles:

• Both particles produced from the particle being split share its thermodynamic state (to comply with the conservation laws);

• The index assigned to the daughter particle is the same as the index of its mother particle, which also indicates that the particle is a daughter;

• The mother particle “knows” nothing about its daughter particles;

• One mother particle may have several daughter particles;

• One daughter particle may have only one mother particle;

• The daughter particles may also split, but all the subsequent daughter particles will “remember” the index of their initial mother particle;

• The daughter particles are placed at the common donor/acceptor side to ensure their quickest possible combination with their respective mother particles.

6. Algorithm of interaction between particles and particle-free materials

The algorithm for volume flux calculations has a modification to deal with mixed cells containing materials with and without particles.

Consider the case of a cell filled with heterogeneous materials, one described only by grid quantities, and the other, by particles (1 and 2, respectively, in Figure 4). Suppose we need to divide the flux moving from left to right (shown with a darker color) between the materials. The volume flux leaving the cell is first filled with the volume of migrating particles. If the migrating particles’ total volume exceeds this volume flux, then the above particle splitting algorithm is engaged (see Section 4, case δV < 0).

Otherwise, the missing part of the outflowing volume flux is filled with the particle-free material. If there are several particle-free materials, the volume is distributed among them based on the VOF algorithm [1]. If the particle-free materials are still not enough, the remaining volume flux is filled with particles using the splitting algorithm for the case of δV > 0 from Section 4.

7. Near-interface algorithms

As part of the proposed method, we have developed an algorithm involving the particles located only near the interface. The region near the interface includes mixed cells and one layer of adjacent pure cells of each material on each side.

Figure 5 shows possible particle layouts relative to the interface. The dark and light cells (Figure 5a) are pure, and the intermediate-color cell (Figure 5b) is mixed. The particles in the cells are marked with a contrasting color. Figure 5, a illustrates the case when both materials are represented by particles only near the interface at t = 0, and Figure 5b shows the particle layout during the computation.

Figure 5.

Particle layout near the interface: a) t = 0, b) t > 0.

Thus, if the material is represented by particles only near the interface, then the same material can be represented both by particles and without particles depending on the interface location. To simplify the algorithm, the same material described by particles and without particles cannot occupy the same cell. This condition for each cell is formally given by

$${\mathrm{V}}_{\mathrm{\xi }}=\sum _{\mathrm{\xi p}}{\mathrm{V}}_{\mathrm{\xi p}},\text{or}\sum _{\mathrm{\xi p}}{\mathrm{V}}_{\mathrm{\xi p}}=0.$$E16

Let us consider two reasons leading to the case when two different representations of the same material are present in the same cell. Figure 6a,b shows cells filled with the same material. In each case, however, in one of the cells the material is represented by particles (black dots), and in the other, without particles. The darker color shows the volume flux relative to the cells’ total volume; the arrow indicates its direction. Below we describe the unwanted cases and the ways to avoid them.

Figure 6.

Conflicting cases of two representations of the same material: a) volume flux from a cell with particles to a particle-free cell; b) volume flux from a particle-free cell to a cell with particles.

Figure 6a shows a volume flux from a cell with particles to a particle-free cell. In this case, the particle splitting algorithm described in Section 4 stays unchanged in the donor cell, but the particles that were supposed to migrate into the acceptor cell are removed with an update of the thermodynamic state, and the particles staying in the donor cell become ordinary (they are no daughter cells any more if they were). Figure 6b shows a volume flux from a particle-free cell to a cell with particles. In this case, the acceptor cell receives a particle, the volume of which is equal to the volume flux and the state of which is the same as the donor-cell material parameters. The added particle then immediately combines with one of the particles in the acceptor cell. The combination rules are given below.

To preserve the layout, where particles are present only near the interface, the particles lying beyond this region are removed from the cells and new particles are created in the cells appearing in the region near the interface.

8. Particle combination algorithm

To balance the particle splitting algorithm (Section 4), we have developed a particle combination procedure. The latter serves to prevent uncontrolled multiplication of particles as a result of their splitting.

Two particles of the same material within the same cell must be combined if one of the following criteria is met:

• One of the particles is a daughter of the other one;

• Two daughters have the same mother;

• The particles have close (to within a constant) coordinates;

• The particle number exceeds the maximum number specified for the cell;

• One of the particles has a relatively small volume.

Particle combination rules:

• If a daughter is combined with its mother, the resulting particle inherits the mother’s coordinates;

• If two daughters of the same mother are combined (p₁ and p₂), the coordinates of the resulting particle (p) are chosen in accordance with their mass ratio:

$${\tilde{\mathrm{x}}}_{\mathrm{p}}={\mathrm{x}}_{\mathrm{p}1}+\left({\mathrm{x}}_{\mathrm{p}2}-{\mathrm{x}}_{\mathrm{p}1}\right)\left({\mathrm{m}}_{\mathrm{p}2}/\left({\mathrm{m}}_{\mathrm{p}1}+{\mathrm{m}}_{\mathrm{p}2}\right)\right),$$E17

$${\tilde{\mathrm{y}}}_{\mathrm{p}}={\mathrm{y}}_{\mathrm{p}1}+\left({\mathrm{y}}_{\mathrm{p}2}-{\mathrm{y}}_{\mathrm{p}1}\right)\left({\mathrm{m}}_{\mathrm{p}2}/\left({\mathrm{m}}_{\mathrm{p}1}+{\mathrm{m}}_{\mathrm{p}2}\right)\right);$$

• The parameters of the resulting particle are calculated subject to the laws of conservation of their mass, specific internal energy and volume:

$$\begin{array}{l}{\tilde{\mathrm{e}}}_{\mathrm{p}}=\frac{{\mathrm{e}}_{\mathrm{p}1}{\mathrm{\rho }}_{\mathrm{p}1}{\mathrm{V}}_{\mathrm{p}1}+{\mathrm{e}}_{\mathrm{p}2}{\mathrm{\rho }}_{\mathrm{p}2}{\mathrm{V}}_{\mathrm{p}2}}{{\mathrm{\rho }}_{\mathrm{p}1}{\mathrm{V}}_{\mathrm{p}1}+{\mathrm{\rho }}_{\mathrm{p}2}{\mathrm{V}}_{\mathrm{p}2}},\\ {\tilde{\mathrm{\rho }}}_{\mathrm{p}}=\frac{{\mathrm{\rho }}_{\mathrm{p}1}{\mathrm{V}}_{\mathrm{p}1}+{\mathrm{\rho }}_{\mathrm{p}2}{\mathrm{V}}_{\mathrm{p}2}}{{\mathrm{V}}_{\mathrm{p}1}+{\mathrm{V}}_{\mathrm{p}2}},\\ {\tilde{\mathrm{V}}}_{\mathrm{p}}={\mathrm{V}}_{\mathrm{p}1}+{\mathrm{V}}_{\mathrm{p}2}.\end{array}$$E18

9. Particle-to-cell density and energy remapping algorithm

For each cell containing particles, the quantities are remapped as follows:

$$\begin{array}{l}{\mathrm{\rho }}_{\xi }^{\mathrm{n}+1}=\sum _{\mathrm{p}}{\mathrm{\rho }}_{\mathrm{\xi p}}{\mathrm{V}}_{\mathrm{\xi p}}/\sum _{\mathrm{p}}{\mathrm{V}}_{\mathrm{\xi p}},\\ {\mathrm{e}}_{\xi }^{\mathrm{n}+1}=\sum _{\mathrm{p}}{\mathrm{e}}_{\mathrm{\xi p}}{\mathrm{\rho }}_{\mathrm{\xi p}}{\mathrm{V}}_{\mathrm{\xi p}}/\sum _{\mathrm{p}}{\mathrm{\rho }}_{\mathrm{\xi p}}{\mathrm{V}}_{\mathrm{\xi p}},\end{array}$$E19

where summing is performed for the particles of material ξ in the cell.

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: ${\mathrm{\rho }}_0=1$, ${\mathrm{e}}_0$ =0, ${\mathrm{u}}_{\mathrm{x}}=1$, ${\mathrm{u}}_{\mathrm{y}}=1$, no particles are specified; in subdomain 1: ${\mathrm{\rho }}_0=10$, ${\mathrm{e}}_0$ =0, ${\mathrm{u}}_{\mathrm{x}}=1$, ${\mathrm{u}}_{\mathrm{y}}=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 particle-free 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.

Figure 7.

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

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.

Figure 8.

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

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. 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: R₁ = 0.8, ρ₀ = 0.01, e₀ = 0, u₀ = 0, ideal gas with $\mathrm{\gamma }=5/3$. Domain 2: R₂ = 1, ρ₀ = 10, e₀ = 0, ${\mathrm{U}}_0^{\mathrm{R}}$ = − 1, equation of state of Mie-Grueneisen type with constants ρ₀ = 10, c₀ = 4, n = 5, γ = 2. Boundary condition at R2: pressure P = 0. Domain 3: vacuum with P = 0. Units of measurement: ρ - [g/cm³], 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.

Figure 12.

Geometry of test problem 4.

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).

Figure 13.

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.

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

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 shock-capturing 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.

Acknowledgments

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