Gas-monodisperse dust mixtures in smoothed particle hydrodynamics: computing of stiff non-linear drag

Simulation the dynamics of gas-solid mixtures is crucial in many applications: chemical reactor design, evolution of circumstellar discs, etc. Such mixtures where gas is a carrier phase and solid grains are dispersed phase can be simulated as interpenetrating continuum media. The characteristic parameter of this problem is relaxation time between gas and dust velocities. In many applications this parameter varies significantly during the simulation (from small to unity). Moreover, the drag force can be non-linearly dependent on solids to gas relative velocity. These factors place high requirements on the numerical methods for such problems. We propose a novel non-iterative algorithm for Smoothed Particle Hydrodynamics for computing gas-solid mixtures with exchange momentum between phases. Forces (except drag force) are approximated explicitly, while drag force is linearized and velocity relaxation time is approximated explicitly while relative velocity is approximated implicitly. The algorithm was tested on dynamical problems for dusty gas mixtures. It was shown that in the developed algorithm for stiff nonlinear drag we can use temporal and spatial resolution independent of drag parameters.


Introduction
Mathematical modeling is widely used for studies of gas-solid mixtures. These studies solve scientific problems and technological tasks, e.g. reactor design for natural gas processing.
Extensive literature [1,2,3,4] is written about mathematical models of two-fluid gasdispersed phase (gas-dust and gas-droplet) laminar and turbulent flows. In particular, gasmonodisperse dust models were studied in [1,2,4], and gas-polydisperse dust models were proposed in [3] (in the local volume, the dispersed fluid is represented by a finite set of fractions, in each solid particles have one characteristic size). A review of current problems and progress of fluid dynamics and thermophysics of two-fluid flows is presented in [4]. In many developed models dynamics of dispersed phase are computed with equations of gasodynamic type separately from gas dynamics. In our paper we also consider this approach of describing gas-solid mixture dynamics where gas is carrier phase and dust is dispersed phase. The inter-phase interaction in gas-solid mixture is described by relaxation terms: velocity (dynamic) and thermal. The presence of these terms leads to the stiff equations if the relaxation time is much less than dynamic time of the problem. We will say that in this case relaxation time is a small parameter.
Finding numerical solution of problems with stiff relaxation terms is an issue. Several authors have shown that in order to reach acceptable accuracy of numerical solution it is necessary to adjust spatial and time steps according to small parameter value. For example simulations in [20] show that in case of smoothed particle hydrodynamics (SPH) where gas and dispersed phase are modelled with different particle sets there are following steps restrictions where t stop -time of particle to gas velocity relaxation (small parameter), τ -time step, h -smoothing length, defining spatial resolution in SPH method, c s -sound speed in the gas. Necessity to satisfy conditions (1)-(2) makes it impossible to run series of computational experiments for many applications.
In some cases, the search for a numerical solution of the original problem can be replaced with a solution of the simplified problem obtained from the original by expansion with respect to a small parameter (for example, [5]). Particularly it means that velocity of dispersed phase for gas-solid mixture dynamics is found not from initial differential equation but from algebraic relation where u -particles velocity, v -gas velocity, g rel -difference of accelerations that affect gas and dispersed phase. The simplified problem becomes non-stiff, but such a transition to the simplified problem is validated only for the case when the relaxation time is much less than the dynamic time of the carrier phase. Simplified problem does not allow simulating a transient regime with a relaxation time is close to the dynamic time.
At the same time practice of solving problems with stiff relaxation terms has shown that there are methods which do not require spatial and time steps dependant on a small parameter. Number of cases has succeeded in acquiring methods that reached appropriate solution accuracy with substantially bigger τ and h. Such methods for stiff or multi-scale problems are known as asymptotic-preserving schemes [9].
Main ideas for constructing finite-difference asymptotic-preserving schemes are provided in [9,10,11]. Finite-difference schemes for gas-monodisperse particles mixture dynamics with intense inter-phase interaction are reported and studied in [6,7], gas-polydisperse particles -in [8]In Euler finite-difference methods [6,7,8] the ability of scheme to keep asymptotic is achieved with the way of time approximation. In Lagrangian methods where gas and dust is simulated with different sets of particles the situation is more complex. To compute the relaxation term (drag force) we have to interpolate gas and dust quantities defined in different spatial points. So in two-fluid Lagrangian approaches keeping the asymptotic is defined not only by way of computing time differentials but also with space interpolation scheme [15]. For two-fluid SPH authors of [12,13,14] have suggested different ways of such interpolation. Their quantitative comparison when drag force linearly depends on the relative velocity between the gas and solids is given in [16].
It should be noted that for simulation of dusty gas dynamics with SPH authors of [17]- [20] have developed numerical schemes for an alternative mathematical model. Their model of the dynamics of a gas-dust medium is based on a one-fluid approach in contrast to the approach, where gas and disperse particles are modelled as interpenetrating continua and with their own system of equations for each continuum. One-fluid approach means equations are solved for the properties of the gas-dust medium as a whole, herewith the diffusion equation is solved for the concentration of the solid phase. In the one-fluid approach, each SPH particle is a carrier of the characteristics of both gas and dust. Such an approach allows computing intense inter-phase interaction using finite-difference scheme ideas [9]. Moreover it is free of numerical over-dissipation appearing in the implementation of several methods in the two-fluid approach. On the other side, in the one-fluid approach, the fulfillment of the mass conservation law for the dispersed phase is not guaranteed. In addition, there is a hard upper limit on the size of the simulated bodies as compared with the two-fluid approach.
This paper is devoted to the method for computing intense inter-phase interaction in a twofluid SPH. The approach we have presented in [14] is developed here for the case when drag force depends non-linearly on relative velocity between gas and solids. In the section 2 the statement of the test problem of the dynamics of the medium of gas-monodisperse dust with intense and moderate interphase interaction is given. Section 3 describes methods of computing inter-phase interaction: Monaghan -Kocharyan (MK) [13] and Implicit Drag in Cell (IDIC) [14]. Section 4 presents the results of computing the test problem using MK and IDIC methods. It was shown that for the intense inter-phase interaction MK method does not preserve assymptotics when IDIC does. It also shows that IDIC method gives appropriate solution accuracy of the non-linear problem at spatial and time steps determined by solving the equations of pure gas dynamics. Conclusions are given in the section 5.

Problem statement
We will study the properties of numerical schemes on one-dimensional shock tube problem for two-fluid gas-dust mixture: where ρ g and ρ d are volume densities of gas and dust, v and u are velocities of gas and dust, p is the gas pressure, f D is a drag force acting on a unit volume of the medium, e is an internal gas energy associated with pressure as follows: here γ is adiabatic exponent. Drag force is defined as where n d is a volume concentration of dust particles, F D is a drag force acting on a separate particle where s is a cross-section area of the particle (e.g. for a sphere of radius a s = πa 2 ), C D is a dimensionless drag coefficient, which is a function of two dimensionless quantities -Mach number here λ is mean free path of a gas molecule. Ma and Kn can be tied together in Reynolds number Following [21] we will use piece-wise approximation of drag coefficient We will make system (4)-(6) dimensionless by the values of length l * , mass m * and time t * , then the derived dimensional quantities have the form ρ The system will be solved in dimensionless variables. For the resulting system (4)-(5), (6)-(7) at the boundaries of the interval, we define the nonleakage conditions. At the initial moment of time, we set the initial zero dust and gas velocities, the gas pressure jump, the gas and dust density jump: γ l = γ r = 7/5.

Monaghan -Kocharyan method
The classical SPH approximation for the equations of motion (4)-(5) for the gas-dust medium is based on Monaghan-Kocharyan (MK) method [13]. We have implemented it in explicit way, where the term describing drag uses the values of velocities from the previous time step: dv n K aj = ρ n j,d ρ n a,g c n a,s where r ja = r j − r a , η is a clipping constant, η 2 = 0.001h 2 and σ is a factor, determined by the dimensionality of the problem; for 1D problems σ = 1. Here Π ab is the tensor of artificial viscosity [22], W n ab = W (h, r ab ) is the smoothing kernel. For each time moment let us split the whole computational domain into nonoverlapping cells, so that the union of these sets coincides with the whole domain. Let a certain cell contains N gas particles of the same mass m g and L dust particles with the same mass m d , with N > 0, L > 0.
For convenience, we introduce the relaxation time of the particle velocity to the gas velocity: where m p is a particle mass, ρ s is a dust particle material density. Then let us introduce volume-averaged values of t * stop and ρ * d , and set We assume, that in the computing of the drag force, affecting the dust from the gas, the gas velocity is the same within the whole cell and equal to v * , but the dust particles have different velocities (and vice versa). Besides that, we shall compute the drag coefficient and densities, using calculated values from the previous time step, while the relative velocity will be used from the next step. The scheme thus has the following form: If the time derivative in (18)- (19) is approximated by the first order finite difference, then a fast way for computing of u n+1 , v n+1 can be proposed.
We compute averaged ρ * d = L j=1 ρ d,j L and t * stop as follows where for C * D we use formula (13) with volume-averaged Re * , Kn * and Ma * . This quantities compute using u * , v * and c * s = N a=1 γ p a ρ g,a N .

Results
Calculations were performed for the cases of small (a = 5 · 10 −6 ) and big (a = 5 · 10 −3 ) dust particles with a high content of dust in the gas (ε =   Figure 1 shows solution by MK scheme for small dust particles. The upper panels show the values for gas, the lower panels show dust. On the left panels are the densities, on the rightthe velocities. In this case, the time step was τ = 0.0001 (from stability condition (1)). It can be noticed that at h = 0.02 and h = 0.01 distinct diffuse spreading is present which disappears after reducing smoothing length to h = 0.001. Figure 2 shows solution by IDIC scheme for small dust particles. The arrangement of panels in this case is the same as for the previous figure. Time step is defined from Courant condition: with Courant parameter CFL = 0.1, cell size is h cell = 1 2 h. In the current calculation, there is no amplitude attenuation, but at h = 0.02 solution nonmonotonicity that disappears at h = 0.01 is visible. Figure 3 shows solution by IDIC scheme for the big dust particles case with h = 0.01 and h = 0.001. Time step is defined from Courant condition with CFL = 0.1, cell size is h cell = h. The upper panels show the values for gas (from left to right: density, velocity, pressure), the lower panels show density and velocity of dust, Reynolds numbers (also from left to right). The bottom right panel (Re numbers) shows that in this case nonlinear drag regimes are present. Indeed, f D is nonlinearly dependant on relative velocity as seen from formula (13) at Re > 1.
Thus, unlike MK, the IDIC method preserves the asymptotic behavior of the solution and is suitable for finding a solution for both small dust and large solids.

Conclusions
To simulate the dynamics of dusty gas mixtures with two-fluid SPH we propose fast and implicit method IDIC. We found that IDIC method for computing the inter-phase interaction is asymptotic preserving. This means that when computing the intense interaction, the spatial and temporal resolution can be chosen independently of the drag parameters. This property of IDIC method distinguishes it from previously proposed methods for computing inter-phase interaction in the two-fluid SPH. Moreover, we found that IDIC method allows to obtain appropriate accuracy for the case when the drag force non-linearly depends on the relative velocity between the gas and the bodies with spatial and time steps determined by the equations of pure gas dynamics. This means that the proposed method can be recommended for problems with intense and moderate inter-phase interaction and those where the size of the solid phase changes by orders of magnitude.