Two-Phase Gas and Dust Free Expansion: Three-Dimensional Benchmark Problem for CFD Codes

Abstract

In the computational mechanics of multiphase dispersed flows, there is an issue of computing the interaction between phases in a mixture of a carrier fluid and dispersed inclusions. The problem is that an accurate dynamics simulation of a mixture of gas and finely dispersed solids with intense interphase interaction requires much more computational power compared to pure gas or a mixture with moderate interaction between phases. To tackle this problem, effective numerical methods are being searched for to ensure adequate computational cost, accuracy, and stability of the results at an arbitrary intensity of momentum and energy exchange between phases. Thus, to assess the approximation, dispersive, dissipative, and asymptotic properties of numerical methods, benchmark solutions of relevant test problems are required. Such solutions are known for one-dimensional problems with linear plane waves. We introduce a novel analytical solution for the nonlinear problem of spherically symmetric expansion of a gas and dust ball into a vacuum. Therein, the dynamics of carrier and dispersed phases are modeled using equations for a compressible inviscid gas. Solid particles do not have intrinsic pressure and are assumed to be monodisperse. The carrier and dispersed phases exchange momentum. In the derived solution, the velocities of gas and dust clouds depend linearly on the radii. The results were reproduced at high, moderate, and low momentum exchange between phases using the SPH-IDIC (Smoothed Particle Hydrodynamics with Implicit Drag in Cell) method implemented based on the open-source OpenFPM library. We reported an example of using the solution as a benchmark for CFD (computational fluid dynamics) models verification and for the evaluation of numerical methods. Our benchmark solution generator developed in the free Scilab environment is publicly available.

1. Introduction

The coupled dynamics of gas and solid particles underlies a large variety of natural phenomena and technology-driven processes. In this paper, we consider the dynamics of a mixture where the boundary between gas and solid inclusions is irrelevant; that is, solid particles are dispersed in a carrier gas. Such systems are called dispersed media. A number of mathematical models have been designed to investigate coupled gas and particles dynamics at different scales (e.g., [1,2,3,4]). In particular, macro and meso level models are used for dispersed media simulation. The description of a medium at a macro level involves the use of volume-averaged characteristics (density, velocity, internal energy) for both the carrier and dispersed phases. In this case, fluid dynamics (of the entire mixture or separately of the carrier and dispersed phases) are modeled using gas dynamics equations. A meso level description means applying kinetic equations to simulate gas and/or dispersed particles. In kinetic equations, as opposed to fluid dynamics equations, the velocity at any given point of space is not unique but has a spectrum of values.

Let us assume that both carrier and dispersed phases in the mixture are inviscid. If the mass or volume fraction of the dispersed phase is considerable (the mixture is dense), then the mixture dynamics are essentially impacted by the mutual exchange of energy and momentum between phases. Such exchange is modeled by introducing corresponding relaxation terms in the right-hand sides of the momentum and energy equations. Examples of such dense mixtures of solids and gas include fluidized bed chemical reactors (e.g., [1]), astrophysical protoplanetary disks, in which planets are formed (e.g., [5]), and other objects.

In the computational mechanics of dense inviscid dispersed media, computing the interaction between phases constitutes a problem. It is known that simulating the dynamics of a mixture of a carrier gas and finely dispersed solids with fast momentum and energy exchange between phases requires much more computational power compared to pure gas or mixture with moderate interaction between phases. To tackle this problem, effective numerical methods are searched for to ensure adequate computational cost, accuracy, and stability of the results at arbitrary interphase exchange rates. Thus, benchmark solutions of dynamics problems for mixtures with arbitrary interphase interaction are highly sought to evaluate the approximation, dispersive, dissipative, and asymptotic properties of numerical methods. A corresponding system of benchmark tests has been developing in recent decades. It includes problems of particle drift in either a stationary or rotating gaseous medium, classical inviscid pure gas dynamics generalized to a gas and dust medium, and solutions describing phenomena specific to a two-phase fluid. In what follows we review several solutions for gas and dust dynamics that are used to verify the computational fluid dynamics codes (primarily by the astrophysical community).

The authors of [6] introduced a solution for a single particle drifting under the gravitational force of a massive central body in an equilibrium gas cloud with Keplerian differential rotation. The particle does not affect the gas dynamics. This solution was further extended by [7] into the case of a gas cloud with non-zero radial velocity. The authors of [8,9] generalized this solution to a polydisperse dust cloud that is not a passive admixture but exchanges momentum with the carrier gas.

Analytical solutions of propagation of small amplitude acoustic waves in a gas and dust mixture are commonly used to assess the dissipative and dispersive properties of numerical methods. The authors of [10] derived the wave group velocity in a monodisperse dusty gas under the approximation of instant momentum exchange (relaxation equilibrium). The solution of a propagating monochromatic wave in a mixture of gas and monodisperse dust at an arbitrary relaxation time was introduced to the practice of astrophysical simulator verification in [11]. The authors of [8] extended this solution into the case of a mixture of gas and polydisperse dust. The authors of [12,13] generalized the monodisperse and polydisperse mixture solutions, respectively, to including particles buoyancy.

The Sod problem (Riemann problem) of the decay of a strong discontinuity in density, pressure, and internal energy represents a classical test for one-dimensional equations for inviscid gas dynamics. A numerical solution to this problem requires introducing an artificial viscosity or using Riemann solvers. This problem has a generalization to the case of a gas and dust mixture with negligible dust volume. Its non-stationary solution for a mixture under relaxation equilibrium was used as a code verification benchmark, for example, in [14,15] (for monodisperse mixtures) and [16] (for polydisperse mixtures). Additionally, the authors of [17] introduced a stationary solution of the Riemann problem for a monodisperse mixture as a benchmark for numerical code verification.

In the case of two dimensions, the numerical solution for inviscid dusty gas dynamics was tested by simulating the streaming instability that results from the differential rotation of a carrier gas and the momentum exchange between gas and dispersed particles [18]. The problem of streaming instability development in the shearing box approximation was used as a CFD verification test for monodisperse [19] and polydisperse [8] mixtures.

Spherically symmetric problems of inviscid gas dynamics, such as Noh, Guderley, and Sedov’s problems, are classical verification tests for three-dimensional CFD implementations. To our best knowledge, they have not been generalized to a gas and dust medium to date, even for the limiting case of relaxation equilibrium when the gas and dispersed phase velocities almost coincide. This paper handles the spherically symmetric problem of a mixture of inviscid gas and dispersed particles expanding into a vacuum, for which we managed to derive a new benchmark solution with arbitrary values of velocity relaxation time. This problem helps to verify a CFD code in three dimensions even at a coarse spatial resolution and to examine the asymptotic, dispersive, and dissipative properties of numerical methods for two-phase dispersed flows.

In Section 2, a brief overview of previous results of similar problems involving pure gases is provided (Section 2.1), the benchmark solution derivation for a mixture of gas and dust is described (Section 2.2), a code generating the benchmark solution is numerically introduced (Section 2.3), and the comparison of the newly derived solution to previous analytical results is presented (Section 2.4). Section 3 describes the SPH-IDIC method for gas and dust dynamics and its implementation based on the open-source OpenFPM library. In Section 4, we demonstrate that the derived solution can be reproduced at extreme and moderate relaxation time values. The conclusions are outlined in Section 5.

2. Benchmark Solution for Gas and Dust Ball Expansion into Vacuum

2.1. Regular Expansion of a Gas Cloud into Vacuum—Brief Outline of Previous Results

The expansion of inviscid gas into a vacuum under its own pressure is a classical fluid mechanics problem. Solutions to this problem are examined with analytical and numerical methods. In particular, a family of solutions to this problem is well explored analytically that describes the velocity of expanding gas as a linear function of the distance from the symmetry center at rest. Such solutions are called regular, or solutions with a uniform velocity gradient. Solutions for the spherically symmetric gas expansion problem with a uniform velocity gradient have been known from the middle of the twentieth century [20,21,22,23]. Further on, the regular solution of the gas ball expansion problem was generalized to the cases of a three-axial ellipsoid [24] and a rotating two-axial ellipsoid [25]. The general approach to these solutions includes: (1) using spherical coordinates; (2) representing all gas macro parameters as products of functions exclusively depending either on time or on the position of a point in the expanding cloud; (3) integrating the initial partial differential equations by variable separation; (4) proceeding either to an ordinary second-order differential equation for the ball radius or to a system of ordinary differential equations for the radii of the ellipsoid principal axes. In the following years, solutions to the problem of a gas cloud expanding into a vacuum were systematically studied and applied in various subject areas. For example, in [26,27] solutions with a uniform velocity gradient were found for gas dynamics in a stationary self-consistent gravitational field for astrophysics (see reviews in [28,29]). A recent contribution [30] provided a general form of known regular solutions for equations of inviscid gas without additional external forces. A comprehensive overview of the literature on the subject, the consecutive derivation, and classification of the known families of regular solutions for a two-axial ellipsoid were provided.

In this paper, we introduce a spherically symmetric problem of expansion into a vacuum of a mixture of inviscid gas and dispersed particles. Using the variable separation method, we managed to construct a novel benchmark solution for arbitrary velocity relaxation time values. This setting generalizes the known solution of the dynamics of a pure gas ball expanding into a vacuum with a uniform velocity gradient and homogeneous density (see [23,30] p. III.B.1) to the case of a gas and dust mixture.

2.2. Problem Description and Derivation of the Analytical Solution

The diagram of gas and dust medium expansion is given in Figure 1.

A gas and dust medium with mass $M=M_{\mathrm{g}}+M_{\mathrm{d}}$, where $M_{\mathrm{g}}$ is the mass of gas and $M_{\mathrm{d}}$ is the mass of dust particles, is modeled as two concentric spheres with radii $R_{\mathrm{g}0}$ and $R_{\mathrm{d}0}$ respectively, centered at the origin. We assume that the gas ball has zero values of viscosity and heat conductance. For the cloud of dust particles, the averaged values of density and velocity are defined for any point inside the ball. Next, we assume that the gas and particles exchange momentum, however, there is no transfer of mass or energy between the phases. Let the distribution of gas pressure and velocities in the gas and particles be spherically symmetric at the initial time moment. Assume that the system expands freely into a vacuum. This means that the distribution of macro parameters remains spherically symmetric during the expansion. Hence, the dynamics of the system can be described by one-dimensional gas dynamic equations in the spherical coordinates:

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left({\rho }_{\mathrm{g}}{r}^{2}v\right)}\hfill & =& 0,\hfill \\ {\displaystyle {\rho }_{\mathrm{g}}\left(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}\right)}\hfill & =& {\displaystyle -\frac{\partial p}{\partial r}+\frac{{\rho }_{\mathrm{d}}(u-v)}{t_{\mathrm{stop}}},}\hfill \\ {\displaystyle {\rho }_{\mathrm{g}}\left(\frac{\partial e}{\partial t}+v\frac{\partial e}{\partial r}\right)}\hfill & =& {\displaystyle -\frac{p}{r^2}\frac{\partial }{\partial r}\left(r^{2}v\right),}\hfill \\ {\displaystyle \frac{\partial {\rho }_{\mathrm{d}}}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left({\rho }_{\mathrm{d}}{r}^{2}u\right)}\hfill & =& 0,\hfill \\ {\displaystyle {\rho }_{\mathrm{d}}(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r})}\hfill & =& {\displaystyle -\frac{{\rho }_{\mathrm{d}}(u-v)}{t_{\mathrm{stop}}},}\hfill \end{array}\right.$$

(1)

where t is the time, r is the distance to the center of the ball, ${\rho }_{\mathrm{g}}(r,t)$ is the gas density, ${\rho }_{\mathrm{d}}(r,t)$ is the volume density of dust particles (loose powder density), $e(r,t)$ is the internal energy of gas, $p(r,t)$ is the pressure, $v(r,t),u(r,t)$ are the gas and dust radial velocities, respectively, and $t_{\mathrm{stop}}$ is the time of relaxation of the particles velocity to the gas velocity. In the most general case, $t_{\mathrm{stop}}$ depends on the particle parameters (shape, size, and material), as well as on the parameters of the gas. To derive an analytic solution, we assume that the parameter $t_{\mathrm{stop}}$ is unchanged during the motion of the gas and dust cloud. It should be noted that the set of Equation (1) is consistent only at $R_{\mathrm{g}}\ge R_{\mathrm{d}}$, where $R_{\mathrm{g}}$, $R_{\mathrm{d}}$ are the radii of gas and dust balls, respectively. To close the system (1), we add the ideal gas equation of state

$${\displaystyle p={\rho }_{\mathrm{g}}e(\gamma -1),}$$

(2)

where $\gamma$ is the adiabatic exponent, which is a free parameter of the problem. We assume now that the pressure gradient and the drag force ensure the motion of the whole assembly, such that for any time moment the gas density in the gas ball and the dust density in the dust ball are independent on the space coordinates, while the ball radii are changing over time:

$${\rho }_{\mathrm{g}}(r,t)=\frac{3M_{\mathrm{g}}}{4\pi R_{\mathrm{g}}^3\left(t\right)}={\rho }_{\mathrm{g}}\left(t\right),$$

(3)

$${\rho }_{\mathrm{d}}(r,t)=\frac{3M_{\mathrm{d}}}{4\pi R_{\mathrm{d}}^3\left(t\right)}={\rho }_{\mathrm{d}}\left(t\right).$$

(4)

Then the continuity equation requires a linear dependence of gas and dust velocities on r. Indeed, from (3) the first equation in (1) can be cast as

$${\displaystyle \frac{1}{{\rho }_{\mathrm{g}}}\frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\frac{1}{r^2}\frac{\partial }{\partial r}\left(r^{2}v\right)=0,}$$

(5)

where the first term is a function of time and does not depend on the space point. From here it follows that the second term is also a function of time and independent of the coordinates, that is

$$\frac{1}{r^2}\frac{\partial r^{2}v}{\partial r}=P_{\mathrm{g}}\left(t\right).$$

(6)

After variables separation, we get the expression for the gas velocity

$${\displaystyle v(r,t)=P_{\mathrm{g}}\left(t\right)\frac{r}{3}+\frac{D_{\mathrm{g}}\left(t\right)}{r^2}.}$$

(7)

To ensure that the gas velocity in the center is well defined, the constraint $D_{\mathrm{g}}\left(t\right)=0$ must hold, which implies a linear dependence of velocity on radius. The same is true for the cloud of dust particles. Thus we found that if we suggest uniform gas and dust density distribution throughout the expansion, we will derive uniform velocity gradient both for gas and dust. We note that the converse statement is not true, i.e., from uniform velocity gradient of gas and dust we can not derive a uniform density of gas and dust throughout the expansion.

The velocity of the ball boundary is the same as the time derivative of its radius:

$$v(R_{\mathrm{g}}\left(t\right),t)={\dot{R}}_{\mathrm{g}}\left(t\right),$$

(8)

$$u(R_{\mathrm{d}}\left(t\right),t)={\dot{R}}_{\mathrm{d}}\left(t\right).$$

(9)

Besides that, the spherical symmetry of the solution implies that the ball center is at rest. This means that, at any time moment, gas and dust velocities are the products of some functions of t and r:

$$v(r,t)=F_{\mathrm{g}}\left(t\right)r=\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}r={\dot{R}}_{\mathrm{g}}\frac{r}{R_{\mathrm{g}}},$$

(10)

$$u(r,t)=F_{\mathrm{d}}\left(t\right)r=\frac{{\dot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}r={\dot{R}}_{\mathrm{d}}\frac{r}{R_{\mathrm{d}}}.$$

(11)

Here and in what follows we write $R_{\mathrm{g}}$, $R_{\mathrm{d}}$ instead of $R_{\mathrm{g}}\left(t\right)$, $R_{\mathrm{d}}\left(t\right)$. Now we add the second and fifth equations from (1):

$${\displaystyle {\rho }_{\mathrm{g}}\left(\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}\right) + {\rho }_{\mathrm{d}}\left(\frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}\right) =-\frac{\partial p}{\partial r}.}$$

(12)

We express the derivatives of velocities along the flow through the products of t and r functions. Taking a derivative of (10) leads to:

$${\displaystyle \frac{\partial v}{\partial t}={\ddot{R}}_{\mathrm{g}}\frac{r}{R_{\mathrm{g}}}-{\dot{R}}_{\mathrm{g}}^2\frac{r}{R_{\mathrm{g}}^2},\frac{\partial v}{\partial r}=\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}},}$$

(13)

then

$$\frac{\partial v}{\partial t}+v\frac{\partial v}{\partial r}={\ddot{R}}_{\mathrm{g}}\frac{r}{R_{\mathrm{g}}}-{\dot{R}}_{\mathrm{g}}^2\frac{r}{R_{\mathrm{g}}^2}+\frac{r{\dot{R}}_{\mathrm{g}}^2}{R_{\mathrm{g}}^2}={\ddot{R}}_{\mathrm{g}}\frac{r}{R_{\mathrm{g}}}.$$

(14)

Similarly, from (11) it follows that

$${\displaystyle \frac{\partial u}{\partial t}+u\frac{\partial u}{\partial r}={\ddot{R}}_{\mathrm{d}}\frac{r}{R_{\mathrm{d}}}.}$$

(15)

Substituting (14) and (15) into (12) brings us to a differential equation for p that can be solved by separation of variables:

$${\displaystyle r\left({\rho }_{\mathrm{g}}\frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}+{\rho }_{\mathrm{d}}\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right) =-\frac{\partial p}{\partial r}.}$$

(16)

Taking into account the kinematic boundary condition at the gas–vacuum boundary

$$p(R_{\mathrm{g}},t)=0,$$

(17)

we find a unique solution of (16)

$${\displaystyle p(r,t)= \left(\frac{R_{\mathrm{g}}^2}{2}-\frac{r^2}{2}\right)\left({\rho }_{\mathrm{g}}\frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}+{\rho }_{\mathrm{d}}\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right) = (\gamma -1){\rho }_{\mathrm{g}}E(1-\frac{r^2}{R_{\mathrm{g}}^2}),}$$

(18)

from which it follows that $p(r,t)$ is proportional to $r^2$. Here

$${\displaystyle E=\frac{R_{\mathrm{g}}^2}{2(\gamma -1)}\left(\frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}+\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right) =E\left(t\right).}$$

(19)

From (18) and the equation of state (2), the internal energy can be cast as a product of two functions depending separately on t and r:

$$e=\frac{R_{\mathrm{g}}^2-r^2}{2{\rho }_{\mathrm{g}}(\gamma -1)}\left({\rho }_{\mathrm{g}}\frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}+{\rho }_{\mathrm{d}}\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right) =E(1-\frac{r^2}{R_{\mathrm{g}}^2}).$$

(20)

By substituting (20) into the third equation from (1), we obtain a differential equation for E:

$${\rho }_{\mathrm{g}}\left(\frac{\partial E}{\partial t}(1-\frac{r^2}{R_{\mathrm{g}}^2})+2E\frac{r^2{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}^3}-2\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}E\frac{r^2}{R_{\mathrm{g}}^2}\right) =-\frac{{\displaystyle {\rho }_{\mathrm{g}}(\gamma -1)(1-\frac{r^2}{R_{\mathrm{g}}^2})E}}{r^2}3\frac{r^2{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}},$$

(21)

which is equivalent to

$$\frac{\partial E}{\partial t}=-3\frac{(\gamma -1)E{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}.$$

(22)

After solving (22) by the variable separation method, we find

$$E=C^{*}{R}_{\mathrm{g}}^{-3(\gamma -1)},$$

(23)

where $C^{*}$ is a constant that is determined from the initial conditions (see (27)). Taking into account (23) and (19), we get a second-order differential equation for $R_{\mathrm{g}}$:

$${\displaystyle C^{*}{R}_{\mathrm{g}}^{-3(\gamma -1)}=\frac{R_{\mathrm{g}}^2}{2(\gamma -1)}(\frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}+\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}).}$$

(24)

Its right hand side contains terms with $R_{\mathrm{d}}$ and ${\ddot{R}}_{\mathrm{d}}$. To derive an equation for these variables, we substitute (10) and (11) into the fifth equation from (1):

$$\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}=\frac{1}{t_{\mathrm{stop}}}(\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}-\frac{{\dot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}).$$

(25)

Therefore, (24) and (25) constitute a closed system of ordinary second order differential equations for $R_{\mathrm{g}}$, $R_{\mathrm{d}}$. To find $R_{\mathrm{g}}$, $R_{\mathrm{d}}$ for a specific time moment, (24) and (25) must be completed with appropriate initial conditions and solved by numerical integration. Then from (10), (11), (18), and (20) with the determined radii the spatial distribution of macro parameters in a mixture can be found for any given time.

2.3. Benchmark Solution Generator

It is straightforward to find a numerical solution of the system (24) and (25) in order to derive $R_{\mathrm{g}}$, $R_{\mathrm{d}}$. However, to facilitate the use of this solution as a benchmark, we provide the implementation of the system solver for (24) and (25) as a program code $MonoDustyBall.sci$, that we developed in the $\mathit{Scilab}$ environment. The code (accessed on 16 December 2021) is available for download by https://github.com/MultiGrainSPH/OneDustBall_Scilab_OpenFPM.

The input parameters for $MonoDustyBall.sci$ routine are:

• physical parameters of the dusty gas medium $\gamma$, $M_{\mathrm{g}}$, $M_{\mathrm{d}}$, $R_{\mathrm{g}0}$, $R_{\mathrm{d}0}$, the initial radial velocity of the gas cloud boundary $v_0$, the initial radial velocity of the dust cloud boundary $u_0$, the internal energy of the gas cloud in its center at the initial time moment $E_0$,
• the time of integration $t^{*}$,
• the calculation parameter $\mathsf{\Delta }t$, which is the time step for integration of the system of ordinary differential equations.

The outputs of the program are:

• the spacial distribution of the dusty gas medium macro parameters ${\rho }_{\mathrm{g}}$, ${\rho }_{\mathrm{d}}$, v, u, p, and e at a given time $t^{*}$,
• the dependence of the ball radii on time $R_{\mathrm{d}}$, $R_{\mathrm{g}}$ at the time interval $[0,t^{*}]$.

Data are delivered as text and in the VTK format for visualization in the ParaView package.

To deliver the output data, the program numerically solves the following system of four ordinary first-order differential equations

$$\left\{\begin{array}{c}\frac{d{R}_{\mathrm{d}}}{dt}={\dot{R}}_{\mathrm{d}},\hfill \\ \frac{d{\dot{R}}_{\mathrm{d}}}{dt}=\frac{R_{\mathrm{d}}}{t_{\mathrm{stop}}}\left(\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}-\frac{{\dot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right),\hfill \\ \frac{d{R}_{\mathrm{g}}}{dt}={\dot{R}}_{\mathrm{g}},\hfill \\ \frac{d{\dot{R}}_{\mathrm{g}}}{dt}=2(\gamma -1)C^{*}{R}_{\mathrm{g}}^{2-3\gamma }-\frac{1}{t_{\mathrm{stop}}}\frac{M_{\mathrm{d}}{R}_{\mathrm{g}}^4}{M_{\mathrm{g}}{R}_{\mathrm{d}}^3}\left(\frac{{\dot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}-\frac{{\dot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}\right),\hfill \end{array}\right.$$

(26)

where $C^{*}$ is a constant that is determined from (23)

$${\displaystyle C^{*}=R_{\mathrm{g}0}^{3\gamma -3}{E}_0.}$$

(27)

This system is equivalent to the system of two second-order Equations (24) and (25). The initial conditions for (26) are

$$R_{\mathrm{d}}\left(0\right)=R_{\mathrm{d}0},{\dot{R}}_{\mathrm{d}}\left(0\right)=u_0,R_{\mathrm{g}}\left(0\right)=R_{\mathrm{g}0},{\dot{R}}_{\mathrm{g}}\left(0\right)=v_0.$$

(28)

The Cauchy problem (26)–(28) is solved by the explicit Euler’s method with $\mathsf{\Delta }t<t_{\mathrm{stop}}$.

2.4. Comparison of the New Solution for Dusty Gas with Previous Analytical Results for Pure Gas

In Section 2.1, a brief review of literature on the regular expansion of a medium into a vacuum is presented. As it follows from this review, an analytical solution was previously known for the expansion of a biaxial gas ellipsoid into a vacuum. In Section 2.2, we derive an analytical solution for the expansion of a gas and dust ball into a vacuum. It is clear that the expansion of a gas ball into a vacuum is a special case of both these solutions. In this subsection, we show that the equations describing these two special cases coincide.

When considering the solution for a gas ellipsoid, we rely on the work in [30]. In this paper, a nonlinear system of second-order ordinary differential equations is presented to find the two principal semiaxes of the ellipsoid $R_{\mathrm{r}}$, $R_{\mathrm{z}}$ at any moment of time. This system represents Equations (65), (66), (49), (72), and (73) from [30] and has the form

$$\left\{\begin{array}{ccc}{\displaystyle -\frac{R_{\mathrm{r}}^{1-2\gamma }{R}_{\mathrm{z}}^{1-\gamma }}{{\ddot{R}}_{\mathrm{r}}}}\hfill & =& {\kappa }_1,\hfill \\ {\displaystyle -\frac{R_{\mathrm{r}}{\ddot{R}}_{\mathrm{r}}}{R_{\mathrm{z}}{\ddot{R}}_{\mathrm{z}}}}\hfill & =& {\kappa }_3,\hfill \\ {\displaystyle {\kappa }_1=-\frac{R_{\mathrm{r}0}^{1-2\gamma }{R}_{\mathrm{z}0}^{1-\gamma }}{{\ddot{R}}_{\mathrm{r}0}},}\hfill \\ {\displaystyle {\kappa }_2=-\frac{R_{\mathrm{r}0}^{2-2\gamma }{R}_{\mathrm{z}0}^{-\gamma }}{{\ddot{R}}_{\mathrm{z}0}},}\hfill \\ {\displaystyle {\kappa }_3=\frac{{\kappa }_2}{{\kappa }_1},}\hfill \end{array}\right.$$

(29)

where $R_{\mathrm{r}0}$, ${\ddot{R}}_{\mathrm{r}0}$, $R_{\mathrm{z}0}$, ${\ddot{R}}_{\mathrm{z}0}$ are the initial radii of the semiaxes ellipsoid and initial accelerations along the major semiaxes. Assuming in (29) that $R_{\mathrm{r}}=R_{\mathrm{z}}=R_{\mathrm{g}}$, we obtain the equation

$${\displaystyle \frac{R_{\mathrm{g}}^{2-3\gamma }}{{\ddot{R}}_{\mathrm{g}}}=\frac{R_{\mathrm{g}0}^{2-3\gamma }}{{\ddot{R}}_{\mathrm{g}0}}.}$$

(30)

Now, let us turn to the system (24) and (25) describing the expansion of a gas and dust ball into vacuum. In (24) we assume that ${\rho }_{\mathrm{d}}=0$, thus obtaining the equation for the gas ball radius

$${\displaystyle \frac{R_{\mathrm{g}}{\ddot{R}}_{\mathrm{g}}}{2(\gamma -1)}=C^{*}{R}_{\mathrm{g}}^{-3(\gamma -1)}\stackrel{\left(27\right)}{=}{E}_0\frac{R_{\mathrm{g}0}^{3\gamma -3}}{R_{\mathrm{g}}^{3\gamma -3}}.}$$

(31)

From ${\rho }_{\mathrm{d}}=0$ and (19) it follows that

$${\displaystyle E_0=\frac{R_{\mathrm{g}0}{\ddot{R}}_{\mathrm{g}0}}{2(\gamma -1)},}$$

(32)

which means that (31) can be written in the form

$${\displaystyle \frac{R_{\mathrm{g}}{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}0}{\ddot{R}}_{\mathrm{g}0}}=\frac{R_{\mathrm{g}0}^{3\gamma -3}}{R_{\mathrm{g}}^{3\gamma -3}}.}$$

(33)

Now, the equivalence of (30) and (33) is clearly seen. Thus, the new solution for gas and dust expansion is consistent with previous results for pure gas.

3. Numerical Solution of the Gas and Dust Expansion Problem

3.1. Numerical SPH-IDIC Method

In this section, we show that the solution of the Couchy problem provided in Section 2 for inviscid gas and dust medium Equation (1) can be reproduced in numerical simulations. Let us note that this solution satisfies one-dimensional Equation (1) in spherical coordinates but is designed rather for three dimensional code verification. We numerically solve the three-dimensional equations of gas and dust dynamics in the Cartesian coordinates

$$\left\{\begin{array}{ccc}{\displaystyle \frac{\partial {\rho }_{\mathrm{g}}}{\partial t}+\nabla \left({\rho }_{\mathrm{g}}v\right)}\hfill & =& 0,\hfill \\ {\displaystyle \frac{\partial v}{\partial t}+(v·\nabla )v}\hfill & =& {\displaystyle -\frac{\nabla p}{{\rho }_{\mathrm{g}}}-\frac{{\rho }_{\mathrm{d}}(u-v)}{{\rho }_{\mathrm{g}}{t}_{\mathrm{stop}}},}\hfill \\ {\displaystyle \frac{\partial e}{\partial t}+(v·\nabla e)}\hfill & =& {\displaystyle -\frac{p}{{\rho }_{\mathrm{g}}}\nabla v,}\hfill \\ {\displaystyle \frac{\partial {\rho }_{\mathrm{d}}}{\partial t}+\nabla \left({\rho }_{\mathrm{d}}u\right)}\hfill & =& 0,\hfill \\ {\displaystyle \frac{\partial u}{\partial t}+(u·\nabla )u}\hfill & =& {\displaystyle \frac{u-v}{t_{\mathrm{stop}}},}\hfill \end{array}\right.$$

(34)

with initial conditions 

$$\left\{\begin{array}{ccc}{\rho }_{\mathrm{g}}(r,0)\hfill & =& {\displaystyle \frac{3M_{\mathrm{g}}}{4\pi R_{\mathrm{g}0}^3},r\le R_{\mathrm{g}0}}\hfill \\ {\rho }_{\mathrm{d}}(r,0)\hfill & =& {\displaystyle \frac{3M_{\mathrm{d}}}{4\pi R_{\mathrm{d}0}^3},r\le R_{\mathrm{d}0}}\hfill \\ v(r,0)\hfill & =& {\displaystyle v_0\frac{r}{R_{\mathrm{g}0}},r\le R_{\mathrm{g}0}}\hfill \\ u(r,0)\hfill & =& {\displaystyle u_0\frac{r}{R_{\mathrm{d}0}},r\le R_{\mathrm{d}0}}\hfill \\ e(r,0)\hfill & =& {\displaystyle E_0(1-\frac{r^2}{R_{\mathrm{g}0}^2}),r\le R_{\mathrm{g}0}}\hfill \end{array}\right.$$

(35)

and the equation of state (2). Here $r=\sqrt{x^2+y^2+z^2}$.

We approximate the system (34) using the SPH [31] (Smoothed Particle Hydrodynamics) method based on interpolating the values of functions and their derivatives on freely moving particles (nodes) using a smoothing kernel. This is a modern Lagrangian approach not relying on the Euler grid. It is widely used for the numerical simulation of a flow in highly deformed media, in particular, for material expansion or destruction. To compute the solution with a small $t_{\mathrm{stop}}$, we approximate the terms responsible for the momentum exchange between gas and dust particles, using the IDIC method (Implicit Drag In Cell) [16,32]. To calculate the drag force, in this method we introduce a spatial grid with cell values computed as the arithmetic means of SPH particle values within a cell. As opposed to approaches applying nonlinear smoothing kernels for drag calculation [33,34,35,36,37], the use of values obtained by linear averaging on a grid allows a fast implicit non-iterative algorithm to be built that ensures the conservation of the momentum in each computational cell and is stable at $\tau \gg t_{\mathrm{stop}}$ where $\tau$ is the time step.

The carrier and dispersed phase densities are calculated conventionally for the SPH method according to the physical definition of density

$${\rho }_{\mathrm{g}i}=\sum _j{m}_{\mathrm{g}}W({\overline{r}}_i-{\overline{r}}_j,h),{\rho }_{\mathrm{d}i}=\sum _j{m}_{\mathrm{d}}W({\overline{r}}_i-{\overline{r}}_j,h),$$

(36)

where ${\rho }_{\mathrm{g}}$ and ${\rho }_{\mathrm{d}}$ are the densities of gas and dust fractions, $m_{\mathrm{g}}$ and $m_{\mathrm{d}}$ are the masses of gas and dust particles respectively, h is the smoothing radius, $W(\overline{r},h)$ is the smoothing kernel, and indices i and j are used for numbering particles. We use here a standard kernel W as a cubic spline for a three-dimensional case [31]

$$W({\overline{r}}_i-{\overline{r}}_j,h)=\frac{1}{\pi h^3}\left\{\begin{array}{cc}1-\frac{3}{2}{q}^2+\frac{3}{4}{q}^3,\hfill & 0\le q\le 1,\hfill \\ \frac{1}{4}{(2-q)}^3,\hfill & 1<q\le 2\hfill \\ 0,\hfill & 2<q,\hfill \end{array}\right.,q=\frac{|{\overline{r}}_i-{\overline{r}}_j|}{h}.$$

(37)

The motion equations for a gas and dust medium have the following form within the SPH-IDIC approach:

$${\displaystyle \frac{v_i^{n+1}-v_i^n}{\tau }=-\sum _j{m}_{\mathrm{g}}\left(\frac{p_j^n}{{({\rho }_{\mathrm{g},\mathit{j}}^n)}^2}+\frac{p_i^n}{{({\rho }_{\mathrm{g},i}^n)}^2}+{\mathsf{\Pi }}_{ij}\right){▽}_i{W}_{ij}-\frac{{\epsilon }^{*n}}{t_{\mathrm{stop}}^{*n}}\left(v_i^{n+1}-u_{*}^{n+1}\right),}$$

(38)

$$\frac{u_j^{n+1}-u_j^n}{\tau }=\frac{1}{t_{\mathrm{stop}}^{*n}}\left(v_{*}^{n+1}-u_j^{n+1}\right),$$

(39)

where v, u are the gas and dust velocities, respectively. The artificial viscosity ${\mathsf{\Pi }}_{ij}$ in the SPH method prevents particle trajectory crossing [31]. During expansion the particles drift apart, and ${\mathsf{\Pi }}_{ij}=0$. Let us denote by $A_i$ the acceleration from all forces acting on a gas particle, except for drag (the first term in the right hand side of (38)):

$${\displaystyle A_i=-\sum _j{m}_{\mathrm{g}}\left(\frac{p_j^n}{{({\rho }_{\mathrm{g},\mathit{j}}^n)}^2}+\frac{p_i^n}{{({\rho }_{\mathrm{g},i}^n)}^2}+{\mathsf{\Pi }}_{ij}\right){▽}_i{W}_{ij}.}$$

(40)

We use here values that are averaged across cells and marked with an asterisk:

$${\epsilon }^{*}=\frac{m_{\mathrm{d}}L}{m_{\mathrm{g}}N},t_{\mathrm{stop}}^{*}=\frac{{\sum }_{j=1}^L{t}_{\mathrm{stop},\mathit{j}}}{L},v_{*}=\frac{{\sum }_{i=1}^N{v}_i}{N},u_{*}=\frac{{\sum }_{j=1}^L{u}_j}{L},A_{*}=\frac{{\sum }_{i=1}^N{A}_i}{N},$$

(41)

where N, L are the numbers of gas and dust particles in an IDIC cell. The average velocities in the layer $n+1$ in an IDIC cell are computed from the solution of the following system [32]:

$$\frac{v_{*}^{n+1}-v_{*}^n}{\tau }=A_{*}^n-{\epsilon }^{*n}\frac{v_{*}^{n+1}-u_{*}^{n+1}}{t_{\mathrm{stop}}^{*n}},\frac{u_{*}^{n+1}-u_{*}^n}{\tau }=\frac{v_{*}^{n+1}-u_{*}^{n+1}}{t_{\mathrm{stop}}^{*n}}.$$

(42)

Next, the velocities of all particles in the $(n+1)$-st layer are computed:

$$\left(\frac{1}{\tau }+\frac{{\epsilon }^{*n}}{t_{\mathrm{stop}}^{*n}}\right)v_i^{n+1}=\frac{v_i^n}{\tau }+\frac{{\epsilon }^{*n}}{t_{\mathrm{stop}}^{*n}}{u}_{*}^{n+1}+A_i^n,\left(\frac{1}{\tau }+\frac{1}{t_{\mathrm{stop}}^{*n}}\right)u_j^{n+1}=\frac{u_j^n}{\tau }+\frac{1}{t_{\mathrm{stop}}^{*n}}{v}_{*}^{n+1}.$$

(43)

The last step is to solve the internal energy equation

$${\displaystyle \frac{e_i^{n+1}-e_i^n}{\tau }=\frac{m_{\mathrm{g}}{p}_i^n}{{({\rho }_{\mathrm{g},i}^n)}^2}\sum _j\left(v_i^n-v_j^n\right){▽}_i{W}_{ij}+\frac{m_{\mathrm{g}}}{2}\sum _j{\mathsf{\Pi }}_{ij}\left(v_i^n-v_j^n\right){▽}_i{W}_{ij}.}$$

(44)

In the first step of the calculation, we split the entire computational domain into non-overlapping volumes, such that the concatenation of these volumes constitutes the whole domain. This step is performed once at the beginning of the main calculation cycle, as the splitting does not change with time. We use Cartesian cubic cells with linear dimensions equal to $\varkappa h$, where $\varkappa$ refers to a free numerical parameter of the problem.

Assuming that we know the coordinates, velocities, and densities of the gas and dust particles in the time layer n, we calculate these values for the layer $n+1$.

• Computing the acceleration $A_{*}^n$ (40) from all forces acting on gas particles in the layer n except for the drag force. This step is identical to that in the conventional SPH method and is performed in parallel. The acceleration computation procedures run concurrently to process the array of particles and are invoked in a cycle across neighboring cells of the auxiliary grid used for the search of neighbors.
• Computing average values within each cell in the layer n. This step requires a single run across all particles and is performed in parallel.
• Computing average velocities across cells in the layer $n+1$ is performed in parallel.
• Computing velocities in the layer $n+1$ is performed in parallel for each phase. The phases are handled consequently, that is, the calculation is first run in parallel for gas particles, and then for dust particles.
• Computing new coordinates, densities, pressures, and energies is the same as in the traditional SPH method and performed as parallel procedures.

3.2. Implementation of Boundary Conditions

One of the standard approaches to modeling the boundary conditions within the SPH method is to introduce boundary particles [38], for which the parameters are calculated in a specific way. In our case, gas and dust boundary particles were arranged around the gas and dust ball according to the same principle. They formed a layer of thickness above $2h$, and all the parameters were set according to the same expressions (3), (4), (10), (11), (18), and (20) for these particles. The boundary particles were treated as regular particles in the calculation but their motion equations were different.

As a result, the energy and pressure values for gas boundary particles turned to be negative, which allowed the pressure gradient at the vacuum boundary to be correctly calculated in a standard way within the SPH approach (38). Without boundary particles, the inner particles near the boundary would have an insufficient number of neighbors, and the density would numerically drop. Since boundary particles participated in the summation for density calculation in (36), this allowed such drop to be avoided, and constant density to be maintained across the space.

A substantial difference between the boundary and inner particles was in the velocity computation procedure. Boundary particle velocities were found by extrapolation from the velocities of inner particles. As a result, the outer layer of boundary particles expanded synchronously with the ball. The extrapolation consisted of finding the coefficients $F_{\mathrm{g}}$ (10) and $F_{\mathrm{d}}$ (11) in the expressions for velocities either as average values across inner particles, or using the least-squares method. Thus, we used the specific precondition of the solution, but the extrapolation could be successfully done from the nearest neighbors, without knowing a priori the expansion condition.

3.3. SPH-IDIC Implementation in OpenFPM

The SPH-IDIC algorithm was implemented using the open-source OpenFPM library [39]. The library represents a set of functions for parallel multiprocessor computing based on the MPI multiprocessor technology for the classes of hybrid particle-grid (SPH, Particle Mesh, etc.) and grid methods. The computational formulae of the method were implemented based on the OpenFPM data structures, and the library further handled parallel computations and provided the dynamic balancing according to its inner algorithms. In particular, the data structure aggregate implemented in OpenFPM used together with the structure aggregate, allows for a handy implementation of the SPH algorithm for gas dynamics with a dust component included. Data in this structure are automatically distributed in the computer memory for any number of processor cores involved in the simulation. Necessary macro parameters are referred to in a self-explanatory manner; for example, to handle the y component of the velocity of a gas particle with index a, gas.getProp<velocity>(a)[Y] should be written. Such transparency facilitates code support and modification.

To implement the SPH method, we need a tool for finding particles located within a given distance from a picked point. For this purpose, OpenFPM provides for a dedicated data structure that represents partition of the computation domain into cubic cells of a given size. In doing so, the list of particles located within any specific cell is known for any given time moment. In the course of simulation the particles can move between cells and the library adequately handles them performing necessary transfers between the MPI processes assigned for a specific subdomain. We use this data structure and its supporting methods to describe the IDIC cells.

The varying parameters in the numerical simulations [40] were: the particle number (from 1 million to 27 million for each phase), the number of MPI processes, and the smoothing radius h of the SPH method. In the computational complexity paradigm, the last parameter determines the cell size for the nearest neighbors search for a given particle: the larger is the cell size, the bigger is the number of neighboring particles to be handled for a given particle (this means that the run time scales as a square of the number of neighbor particles). For benchmark runs, we used a server computer with a 48 core AMD EPYC 7002 Series processor with a speed of 2.3–3.3 GHz (96 virtual cores with hyper-threading technology) and 192 Gb of RAM.

The benchmark tests have shown that the designed parallel implementation based on the OpenFPM library is a promising solution for numerical simulations with the number of model particles in the order of 100 million, run on a multiprocessor supercomputer with distributed memory and a large number (200 or more) of MPI processes.

The source code of this implementation (accessed on 16 December 2021) is available by https://github.com/MultiGrainSPH/OneDustBall_Scilab_OpenFPM.

4. Benchmark Problem for Arbitrary Relaxation Time for SPH-IDIC in Three Dimensions

4.1. Physical and Numerical Setup of Models

In this section, we show that the derived functions (3), (4), (10), (11), (18), and (20) with $R_{\mathrm{g}}$, $R_{\mathrm{d}}$ satisfying Equations (24) and (25) can be used as a simple and comprehensive tool to test various numerical methods of handling the system (34). For this purpose, we consider three models of expansion of a gas and dust ball into vacuum with different $t_{\mathrm{stop}}$. In all three models, we pinpoint $M_{\mathrm{g}}=1000$, $\gamma =4/3$, the initial radii of the gas and dust balls $R_{\mathrm{g}0}=R_{\mathrm{d}0}=1$, the initial radial expansion velocity ${\dot{R}}_{\mathrm{g}}={\dot{R}}_{\mathrm{d}}=v_0=u_0=1$, and the value of the internal energy in the center of the ball $E_0=18.61$.

The variable parameters in the models are the time of velocity relaxation and the mass of the dust ball. In the baseline run, at a moderate $t_{\mathrm{stop}}=0.1$ value (M model), we set $M_{\mathrm{d}}=0.04M_{\mathrm{g}}$, which corresponds to a mixture with low dust content. It is known (see [16,37,41]) that the computational cost of simulating gas and dust mixture dynamics especially increases in case of dense media with short relaxation times, compared to simulating pure gas dynamics. Therefore, for simulations with small $t_{\mathrm{stop}}={10}^{-3}$ (S model) and large $t_{\mathrm{stop}}={10}^3$ (L model), we set $M_{\mathrm{d}}=0.8M_{\mathrm{g}}$, which is realistic for a dense gas and dust medium.

The computational domain for all models is the cube ${[-1.8:1.8]}^3$. At the initial time moment, we locate the model particles in the nodes of a uniform Cartesian grid with a step $1/20$. The particles representing the gas ball (inner particles) are placed inside the domain $r\le 1$, and the particles modeling the gas–vacuum boundary occupy the area $1<r\le 1.21$. The mass of an inner gas particle is defined as $m_{\mathrm{g}}=M_{\mathrm{g}}/N_{\mathrm{g}}$, where $N_{\mathrm{g}}$ is the number of inner gas particles, and the mass of a boundary gas particle is set to be equal to $m_{\mathrm{g}}$. Dust particles are set up and their masses are assigned in a similar way.

In the baseline model M, let us assume that $\tau =0.002$ (this corresponds to the Courant–Friedrichs–Lewy number CFL = 0.2) and $\varkappa =1$. The numerical parameters for S and L models are selected from the following considerations. It is known that sufficient solubility conditions for the explicit classical method of computation of interphase interaction in the SPH method [15,35] are as follows:

$$\tau <t_{\mathrm{stop}},$$

(45)

$$h<c_{\mathrm{s}}{t}_{\mathrm{stop}},$$

(46)

where $c_{\mathrm{s}}$ is the sound speed in pure gas. If the condition (45) fails, the numerical scheme is unstable, and failure of the condition (46) results in an overdissipation of the solution. Satisfying (46) in three-dimensional simulations means that a tenfold decrease of a small value $t_{\mathrm{stop}}$ leads to an increase in the computational cost at a minimum by ${10}^3$, provided that other model parameters are kept constant. The IDIC method that we used in simulations ensures stability and helps prevent overdissipation of solutions, even if (45) and (46) fail. These properties were previously established [16,32] during numerical one-dimensional simulations of the propagation of linear and nonlinear waves in a gas and dust medium. To show that the IDIC approach holds the same stability properties for three-dimensional simulations, we use $h=0.1$ for all models, thus violating the condition (46) across the entire ball volume in model S and within the outer ball layers in models M and L. For models S and L, we use $\tau =0.005$, thus purposely violating (45) for model S. Furthermore, it was practically shown in [16,32] that the least noisy solutions can be obtained at $\varkappa =0.5$ for one-dimensional problems with extreme relaxation times. We pick this value of $\varkappa$ for three-dimensional simulations. All numerical and physical parameters of the M, S, and L models are summarized in

Table 1. Numerical and physical parameters of the models. N is the number of particles of one phase along the ball diameter. Parameters differentiating between the models are shown in bold, and identical parameters for all models are in normal font.

Model	N	h	$\mathit{\varkappa }$	$\mathit{\tau }$	${\mathit{t}}_{\mathbf{stop}}$	$\mathit{\gamma }$	${\mathit{M}}_{\mathbf{g}}$	${\mathit{M}}_{\mathbf{d}}$	${\mathit{R}}_{\mathbf{g}0}={\mathit{R}}_{\mathbf{d}0}$	${\mathit{v}}_0={\mathit{u}}_0$	${\mathit{E}}_0$
M	40	0.1	1.0	$\mathbf{0.002}$	${\mathbf{10}}^{-\mathbf{1}}$	4/3	1000	40	1.0	1.0	18.61
S	40	0.1	0.5	$\mathbf{0.005}$	${\mathbf{10}}^{-\mathbf{3}}$	4/3	1000	800	1.0	1.0	18.61
L	40	0.1	0.5	$\mathbf{0.005}$	${\mathbf{10}}^{\mathbf{3}}$	4/3	1000	800	1.0	1.0	18.61

.

Let us compute the initial density for the collocated particles from (36), and other macro parameters of the medium from (10), (11), (18), and (20).

4.2. Ball Expansion into Vacuum at Finite Velocity Relaxation Times

Foremost, let us ensure that the solution from Section 2 is reproduced in the regime with a moderate relaxation of time and small content of the dispersed phase. In computational multiphase fluid dynamics, this is the most simple case in terms of calculating the interphase interaction. We simulate the expansion of a gas and dust ball according to the M model until the time moment $t=0.25$. The results of this run are shown in Figure 2. The panels in the figure sequentially represent (from left to right and from top to bottom) gas density (a), dust density (b), gas pressure (c), gas velocity (d), dust velocity (e), and radii of the gas and dust balls, respectively (f). The open circles show the computed values for the time moments $t=0$, $t=0.05$, $t=0.1$, and $t=0.25$; the solid black line outlines the analytical values at the final time moment $t=0.25$. It is seen that for all time moments the computed density of gas and dust is uniform, and computed radial velocities of gas and dust are linear in radius. At the control time point $t=0.25$, the numerical values coincide with the analytical solutions for all macro parameters of the gas and dust medium. In this regime, since $t_{\mathrm{stop}}$ is comparable to the time of the ball expansion, the expansion velocity of gas turns to be higher than that of dust, and the radius of the gas ball increases faster compared to that of the dust ball. Specifically, this means in particular that the drag force undergoes a discontinuity at the boundary of the dust ball: indeed, the drag force inside the dust ball is non zero, while it vanishes outside the ball since ${\rho }_{\mathrm{d}}=0$. We note that the discontinuity in the drag force acting on a gas and dust mixture does not necessarily imply a weak or strong discontinuity of the macro parameters inside the domain filled with material. In future work, we shall investigate in detail the numerical challenges that may encounter during simulations of the dust ball boundary inside the gas ball.

It is seen in Figure 2b that, at the final time moment, the density of the dust ball is more contaminated with numerical noise compared to the dust density at earlier moments, and to the gas density at the final moment. On the other hand, the observable noise level is rather acceptable considering the coarse spatial resolution used. Besides that, the design specifics of the IDIC method prescribes the use of $\epsilon$ value from (41) in the calculation of dust and gas velocities, rather than the dust density computed from (36) and presented in Figure 2. Therefore the inaccuracy in the dust density calculation does not contribute to the noise contamination in the computation of gas and particles velocities, which can be seen in Figure 2d,e. To control the numerical noise of the SPH-IDIC method, the dispersion properties of this method will be analyzed in the following work.

4.3. Ball Expansion into Vacuum at Infinitely Small and Infinitely Large Velocity Relaxation Times

Let us show now that the devised solution can be reproduced in numerical simulations with extreme relaxation times. An infinitely large relaxation time corresponds to the frozen equilibrium mode with no drag between phases. An infinitely small relaxation time corresponds to the mode of relaxation equilibrium when the drag force is so large that the gas and dust propagate jointly as an integral medium.

First, we compute the expansion of a gas and dust ball according to model L with $t_{\mathrm{stop}}={10}^3$ until the time moment $t=0.25$. The simulation results are shown in Figure 3. The notations are the same as in the previous Figure 2. In Figure 3a,b it is seen that the gas and dust densities in the L model are uniform at all times, and the radial velocity of each phase depends on the radius linearly. At the reference time point $t=0.25$, the numerical values coincide with the analytical solutions for all macro parameters of the gas and dust cloud. It follows from Figure 3e,f that the dust ball is not affected by any forces, its center is at rest, and the dust particles of the cloud expand straight and with constant velocity. This means that frozen equilibrium conditions truly take place in model L at the times involved. The propagation of the gas ball material is affected by the pressure gradient, therefore it expands with acceleration (see the left and right panels in the second row). The comparison of the panels (f) in Figure 2 and Figure 3 shows that the gas balls expand with comparable velocities in models M and L. This is explained by the fact that the drag is almost absent in model L, and the mass content of the dust in model M is too small to exert a considerable impact on the momentum of the gas ball. However, dust balls expand with different velocities in models M and L. This is explained by the fact that, in the M model, the dust is dragged by the gas matter, so that the dust expansion velocity is higher in model M compared to model L.

Now we compute the expansion for the configuration of model S with $t_{\mathrm{stop}}=0.001$ until the same time moment $t=0.25$. The simulation results are shown in Figure 4 with the same legend as for Figure 2. From Figure 4e it is obvious that relaxation equilibrium conditions truly take place in model S, since the gas and dust ball radii coincide during the entire simulation. The panels (a), (b), (d), and (e) show that the gas and dust ball expands into vacuum in model S similarly as in M and L, with uniform densities and linear velocities maintained at all times. In the final time moment $t=0.25$, the computed values coincide with those from the analytical solution for all macro parameters of the gas and dust medium. This means that the IDIC method provides a high accuracy in 3D simulations of the interphase interaction in dense mixtures even if the conditions (45) and (46) fail. We note that previously, a three-dimensional model of gas and dust ball expansion with relaxation equilibrium was simulated using the SPH-IDIC method in [40,42]; however, no benchmark solution was provided for the simulated configuration to evaluate the quality of the results.

The spatial distribution of gas and dust SPH particles is shown in Figure 5 and Figure 6. The expansion of the gas–dust sphere is clearly visible, while the Cartesian uniformity of the distribution of particles is preserved. The color shows the pressure in gas particles and the velocity in dust particles.

Furthermore, the comparison of panels (f) in Figure 3 and Figure 4 shows that the ball in the relaxation equilibrium state expands slower compared to the gas ball under frozen equilibrium conditions. This dependence correlates with a known behavior of propagating sound waves in a gas and dust material. Under relaxation equilibrium conditions, the group velocities of plane waves in the pure gas $c_{\mathrm{s}}$ and in the gas and dust mixture $c_{\mathrm{s}}^{*}$ are related as

$${\displaystyle c_{\mathrm{s}}^{*}=\frac{c_{\mathrm{s}}}{\sqrt{{\displaystyle 1+\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}}}}}.}$$

(47)

This expression indicates that waves in a mixture propagate slower compared to a pure gas. Let us note that the corresponding factor ${\displaystyle (1+\frac{{\rho }_{\mathrm{d}}}{{\rho }_{\mathrm{g}}})}$ is present in Equation (19) under relaxation equilibrium, since it turns out in this case ${\displaystyle \frac{{\ddot{R}}_{\mathrm{g}}}{R_{\mathrm{g}}}=\frac{{\ddot{R}}_{\mathrm{d}}}{R_{\mathrm{d}}}}$.

4.4. Discussion, Results and Limits of the Study and Future Plans

In this paper, we have drawn a spherically symmetric solution for the problem of gas and dust ball expansion into vacuum. The problem was set for a system of two-fluid dynamics equations for a two-phase medium (1). For the gas–vacuum boundary, kinematic (8) and dynamic (17) conditions were defined, while for the dust ball boundary only the kinematic condition was set. The derived solution satisfies the specific initial conditions where the gas and dispersed phase densities are uniform across the space, the velocities of gas and the dispersed phase depend linearly on the radii and the dispersed phase pressure is a quadratic function of the radius. As far as it is known to the authors, no field measurements have been conducted to date that would validate the use of the two-phase fluid dynamics model for such a problem setting.

Nevertheless, considering the specific numerical simulation of ultradispersed medium dynamics, the derived solution is of interest for numerical code verification. A verification benchmark would help demonstrate that numerical methods of solving two-phase fluid dynamics equations (in the form of (1), (34), or another) may represent the analytical solution with an appropriate accuracy. Since the solution is presented for the first time, it needs to be carefully proved before it is recommended as a benchmark. In this paper, besides deriving it from the initial equations, as presented in Section 2.2, the following steps were performed to validate the solution. In Section 2.4 we show that a specific case of the derived solution correlates to a specific case of the known solution for the problem of gas ellipsoid expansion into vacuum suggested by other authors [30]. In Section 4.2 and Section 4.3 the values of the problem parameters are provided for which the derived solution is proved to be reproducible by numerical models. The numerical method used to reproduce the solution as well as its software implementation was developed by the authors of this paper. In the light of increased interest in the simulation of gas dispersion medium dynamics, we hope that this solution will be numerically reproduced in the nearest future by other authors using various numerical methods.

A promising development of this research would be a derivation of such benchmark solutions for the expansion of a two-fluid gas dispersion medium into a vacuum with no drag force discontinuity at the dust cloud boundary, as well as with no spherical symmetry requirement. It should be noted that a solution to the problem of the interpenetration of two separate fluids is known where the dispersed phase has its own pressure and buoyancy [43].

5. Conclusions

We introduced here the derivation of an analytical solution for the problem of spherically symmetric expansion of a gas and dust ball into a vacuum. A ball containing both gas and dust particles is modeled as a two-phase continuum with different velocities. The carrier phase in the system is a non-viscous and non-heat-conducting gas, and the dispersed phase represents a cloud of monodisperse particles with negligible volume. The carrier phase is described by the Euler equations and the dispersed phase by equations of continuity and motion of compressible gas without intrinsic pressure. The carrier and dispersed phases exchange momentum, and the exchange rate is described by the velocity relaxation time. The derived solution describes specific expansion conditions where radial velocities of the gas and dispersed phases depend linearly on the radii and the density of each phase is uniform across the space. The solution to the expansion problem with linear velocity and uniform density is derived for arbitrary relaxation time; hence, it extends the previously known solution into a more general case.

It is shown that the obtained solution is well reproduced in three-dimensional simulations at infinitely large, moderate, and infinitely small values of velocity relaxation time. The numerical simulations were performed with the SPH-IDIC method implemented in the open-source OpenFPM library. We suggest a practical use of this solution as a verification benchmark for computational flow dynamics models and investigating the properties of numerical methods. For convenience, we make the benchmark solution generator designed in the $\mathit{Scilab}$ environment available to the public.