Granular and multiphase (gas–solids) kinetic theory-based models have emerged a leading modeling strategy for the simulation of particle flows. Similar to the Navier–Stokes equations of single-phase flow, although substantially more complex, kinetic theory-based continuum models are typically solved with computational fluid dynamic (CFD) codes. Under the assumptions of the so-called homogeneous cooling state (HCS), the governing equations simplify to an analytical solution describing the “cooling” of fluctuating particle velocity, or granular temperature. The HCS is used here to verify the implementation of a recent multiphase kinetic theory-based model in the open source mfix code. Results from the partial verification test show that the available implicit (backward) Euler time integration scheme converges to the analytical solution with the expected first-order rate. A second-order accurate backward differentiation formula (BDF) is also implemented and observed to converge at a rate consistent with its formal accuracy.

## Introduction

Understanding the behavior of solid particles suspended in a fluid is paramount to a variety of industrial and engineering fields. Unfortunately, gas–solid flows are incredibly complex, yielding a variety of structural patterns that can span several orders of magnitude in size, ranging from the system scale down to a few particles. Numerical strategies to track the motion of individual particles, such as computational fluid dynamics-discrete element model (CFD-DEM) or direct numerical simulation (DNS) (which also uses DEM for the particles), are able to faithfully reproduce the complicated behavior, but at a considerable computational cost. Continuum models that treat the particles as a continuous solids phase coexisting with and interpenetrating the gas phase (if considered) must be applied to simulate larger systems of engineering importance.

A common approach to derive continuum models relies on an analogy between particles and the molecules, which constitute traditional fluids. Starting from the Enskog kinetic equation (the Boltzmann equation equivalent for finite sized particles), a continuum theory for the solid phase complete with explicit closure relations can be derived employing a Chapman–Enskog expansion method. Such an approach can be used to derive both granular (no interstitial fluid) and multiphase models. As should be expected, kinetic theory-based models and their closures can be rather involved presenting a significant challenge to code verification.

One of the most well-known and arguably simplest cases of particulate flow is the homogeneous cooling state (HCSt). Originally proposed by Haff [1] for a granular system, the HCSt describes the cooling of a system of particles of infinite expanse in the absence of gravity or mean flow that is initially uniformly excited. In other words, the velocity of the particles is characterized by the initial granular temperature alone and the mean velocity in each direction is zero. The granular temperature is a measure of the fluctuating kinetic energy of the particles. An analytical solution to (continuum) kinetic theories can be obtained when the assumptions of the HCSt are applied, which will be used for the purposes of code verification in this work.

When solving the HCSt numerically, a finite domain size is used and periodic boundary conditions are applied in each direction. Although the HCSt and the homogeneous cooling system (HCSys) are occasionally used interchangeably (possibly because they have the same HCS acronym), there is a significant distinction, which is of importance here (hence the separate acronyms). Due to the dissipative nature of particle collisions [2] and/or the viscous dissipation of an interstitial fluid [3], the homogeneous cooling state is linearly unstable, which causes a system-size dependence [4,5]. Similar to the transition to turbulence in pipe flow, if the system is sufficiently small the instability is suppressed, but if the system is larger than a critical length scale, instabilities develop commonly referred to as clustering, see Ref. [6] for a review. As may be expected, once the instabilities become sufficiently large in amplitude, the inhomogeneous, clustered state no longer obeys the analytical solution to the HCSt [2,3,7–9]. Care is exercised here to avoid clustering so that direct comparison to an analytical solution is possible.

Both the HCSt and clustered (inhomogeneous) HCSys are commonly used in the granular and gas–solids flow communities. Some examples include:

- (1)
The HCSt is frequently used for the verification and validation of the transport coefficients derived from kinetic theory, e.g., see Refs. [10–12]. These studies typically utilize direct simulation Monte Carlo (DSMC) to generate the numerical data; the results can provide insight into the correctness of the derivation and the range of validity, but not on the validity of the statistical mechanics approach (kinetic equation) itself which DSMC also relies on.

- (2)
The critical length scale necessary for onset of instabilities in the HCSys has been used more recently to quantitatively assess the ability of kinetic theory-based models to predict instabilities, e.g., see Refs. [5], [13], and [14]. Comparisons are typically made to numerical data generated from DEM; unlike DSMC, DEM involves balances for each particle, instead of relying on a kinetic equation.

- (3)
Since the behavior of the granular HCSt and HCSys is relatively well understood for ideal particles (smooth, monodisperse, etc.), the HCSt and HCSys offer an attractive starting point to explore additional physics. Some examples include: variable coefficient of restitution [15], particle friction [16], relative impact of a viscous fluid [9], and irregularly shaped particles [17].

- (4)
The HCSt is commonly used as a benchmark problem for DEM simulations [9,18]. Although DEM provides more information than continuum models and does not rely on statistical mechanics, comparison to the kinetic theory-based analytical solution of the HCSt can provide some indication if numerical settings (e.g., spring constant, neighbor search frequency, etc.) have been chosen appropriately.

However, the HCSt is not commonly used to verify the codes that use kinetic theory-based continuum models. In this work, we show how the HCSt may be used as a partial verification test for CFD codes. It is our hope that this test will become commonplace among CFD codes that utilize such kinetic theory-based models.

## Governing Equations

The kinetic theory-based two-fluid model of Garzó et al. [19] is used in this work and hereafter referred to as the GTSH model. The theory is derived for monodisperse, spherical, frictionless particles immersed in a viscous fluid with a large density ratio, *ρ** = *ρ _{p}*/

*ρ*≫ 1. The Chapman–Enskog expansion of the kinetic equation is carried out to first-order resulting in a Navier–Stokes-order two-fluid model. Complete details of the GTSH model, its derivation, and complete closures can be found elsewhere [5,19]. In the HCSt, many terms of the model are null, namely all spatial derivatives, which will not be detailed here for the sake of brevity.

_{f}*ϕ*(

**x**,

*t*) =

*ϕ*, no mean flow of the solids (particles) or interstitial fluid in any direction, $Us(x,t=0)=Uf(x,t=0)=0$, and a uniform and nonzero initial granular temperature,

*T*(

**x**,

*t*= 0) =

*T*

_{0}. Under these assumptions, the GTSH equations reduce to simply a single ordinary differential equation (ODE) governing the decay of an initial granular temperature

*e*< 1) dissipation of particle collisions, and the thermal drag,

*e*, characterizes the degree of inelasticity of the particles and takes a value between zero (sticky limit) and unity (elastic limit). More specifically,

*e*gives the ratio of relative velocities between two particles post and precollision and is assumed to be a constant in the GTSH model. The particle diameter,

*d*, is also constant. The radial distribution function at contact,

_{p}*χ*, accounts for volume exclusion effects of the finite sized particles. Here, we use the model of Ma and Ahmadi [20],

where $ReT=\rho fdpT/\mu f$ is the thermal Reynolds number, which characterizes the ratio of fluctuating fluid inertial to viscous forces. Note that Eq. (7) is of a slightly different form than originally reported in Ref. [22]; the original DNS data were used to refit a curve with the proper limiting behavior *R*_{1} → 0 as *ϕ* → 0.

It should be noted that the HCSt greatly simplifies the GTSH model as the continuity and momentum equations of both phases drop out (identically satisfied). Therefore, comparisons to the HCSt cannot be considered a complete verification test in itself as the solution structure is not sufficiently complex to exercise the complete set of equations [23]. A partial verification test with the HCSt—particularly as a first step, however, is not without merit for several reasons: (i) the transport coefficients derived from the Enskog kinetic equation are significantly more complex than typically encountered in single-phase CFD, so it is worthwhile to perform reduced tests in order to be able to eliminate the portions of the model exercised by the HCSt as potential culprits if larger, model-wide verification tests fail; and (ii) some of the terms exercised here, namely *χ* and *a*_{2}, are also widely used in other transport coefficients not exercised by the HCSt; and (iii) if the HCSt verification test does fail, it is relatively straightforward to isolate the problem.

The exact solution of the HCSt, Eq. (11), will be used to verify code calculations.

## Numerical Method

^{2}For more general problems, mfix employs a finite volume discretization with staggered variables, upwind flux-limiter extrapolation schemes, and a SIMPLE-type algorithm for pressure correction. Because the HCSt reduces to a temporal ODE, only the time advancement schemes are of interest here. Two schemes are available in the 2015-1 version of mfix, the most widely used scheme is simply the (backward) Euler method,

A second method is also available in mfix 2015-1 referred to as two-step implicit Runge–Kutta in the user's manual and Crank–Nicolson in the actual source code. In actuality, the method is neither implicit Runge–Kutta nor Crank–Nicolson because the right hand side is only evaluated (implicitly) at the half time-step. A more accurate description of this method is a two-step leap frog scheme centered about the half time-step or an alternating backward/forward Euler half time stepping scheme. While the scheme does appear to be second-order accurate in time, there is no published documentation of the method and therefore it is not suitable for code verification. Furthermore, the nonconvergence of the scheme was recently observed in the mfix code in a more formal verification study using the method of manufactured solutions [24], which was attributed to the lack of an extrapolation step for the gas-phase pressure.

*n*− 1,

*n*− 2,…) [25]. In addition to higher order accuracy, BDF schemes are particularly desirable for CFD applications due to their relatively large region of absolute stability. The (formally) second-order accurate BDF,

*n*= 0 time level, no information of the

*n*− 1 state is known. The start-up issue is solved here by simply using a single backward Euler initial time-step.

Normally, mfix users allow the code to adjust the time-step automatically if convergence stalls or fails. Here, the stall detection algorithm is switched off and the time-step is set at constant values of *h** = 1.0 × 10^{+02}, 1.0 × 10^{+01}, 1.0 × 10^{+00}, 1.0 × 10^{−01}, and 1.0 × 10^{−02}, where time has been nondimensionalized by $h*=hT0/dp$. The criteria for convergence of the implicit evaluation are set to 1.0 × 10^{−10} and the tolerance of the linear solver is set to 1.0 × 10^{−11}. The iterative tolerance of 1.0 × 10^{−10} is several orders of magnitude lower than the default mfix value and also lower than tolerances used in typical applications. The low convergence criteria used here are necessary to observe the convergence of the discrete solution to the exact solution for the small time-step sizes. Although solution variables should be constant, a 3 × 3 × 3 cubic spatial grid is applied to a system size six particle diameters in each direction, *L** = *L*/*d _{p}* = 6, which is sufficiently small to avoid instabilities [5]. The simulation end time is set to $tend*=3800$, allowing the temperature to decay by eight orders of magnitude for the selected conditions, as described in Sec. 4.

## Results and Discussion

The conditions of the HCSt are defined by four nondimensional parameters: the initial thermal Reynolds number, $ReT0=\rho fdpT0/\mu f$ where *T*_{0} = *T*(*t* = 0), the density ratio, *ρ** = *ρ _{p}*/

*ρ*, the restitution coefficient,

_{f}*e*, and the solids volume fraction,

*ϕ*. Here, we use the conditions: Re

_{T}_{0}= 50,

*ρ** = 1 × 10

^{+3},

*e*= 0.85, and

*ϕ*= 0.15. For completeness, we note that these conditions are obtained from dimensional properties:

*d*= 1 cm,

_{p}*ρ*= 1 g/cm

_{p}^{3},

*ρ*= 1 × 10

_{f}^{−3}g/cm

^{3},

*μ*= 2 × 10

_{f}^{−3}g/cms

^{−2},

*T*

_{0}= 1 × 10

^{+4}cm

^{2}/s

^{2}. The exact solution of the HCSt, Eq. (11), with the preceding conditions, is plotted in Fig. 1 for the time period of 38 s nondimensionalized by $T0/dp$. There are two stages of decay, the first dominated by inelastic dissipation (from particle collisions) and the second dominated by viscous dissipation (from the fluid). The transition between the two stages occurs at approximately

*t** = 300 in Fig. 1.

*L*

_{1}and

*L*

_{2}norms. Since the granular temperature is decaying by several orders of magnitude, the error should be normalized based on the instantaneous value so that the norms are not dominated by the errors at the early times. Therefore, we take

where *T ^{n}* is the computed solution at the discrete time,

*t*, and

^{n}*T*(

*t*) is the exact solution at the discrete time and $N=tend*/h*$ is the number of discrete time steps.

^{n}*L*

_{1}and

*L*

_{2}norms is shown in Fig. 2. The Euler method exhibits a first-order convergence rate as expected. The

*L*

_{1}and

*L*

_{2}norms are also reported in Table 1. The order of accuracy between two successive grid refinements, $h1*$ and $h2*$, is approximated as [26]

*i*= 1, 2.

The *L*_{1} and *L*_{2} norms for the BDF scheme are also shown in Fig. 2 and provided in Table 2. The BDF scheme converges approximately at rate consistent with its formal second-order accuracy–albeit less uniformly than in the case of the Euler method. We note that one issue for the BDF scheme is that for the smallest time-step, the error is beginning to saturate to the minimum allowed by the set tolerance levels.

## Summary

This work demonstrated the use of the HCSt, for which an analytical solution is available, as a partial verification problem for kinetic theory-based (continuum) models. Although the HCSt greatly simplifies the complete governing equations, it is still a valuable tool for verification, particularly as a first step. Furthermore, it may be possible to combine the HCSt with other analytical solutions which isolate different terms (e.g., simple shear flow) in order to create a more systematic approach to verification for kinetic theory-based models. The GTSH model [19] considered here is a multiphase continuum, or two-fluid model, which was solved numerically with the open-source mfix code developed at the National Energy Technology Laboratory. The backward Euler method was verified to converge at the expected rate of first-order accuracy. A higher-order time stepping scheme is also available in the mfix 2015-1 code. However, the method was determined to be neither a Crank–Nicolson nor an implicit Runge–Kutta method as reported and no known published literature was found. Therefore, we recommend replacing the scheme with the second-order backward differentiation formula. The BDF2 scheme was implemented in mfix and observed to converge at a rate approximately equal to the formal order of accuracy. Although the new scheme was only implemented for the single equation (granular temperature) exercised by the HCSt, this study provides a blueprint for how the mfix code may achieve second-order temporal accuracy with minimal modifications to the main kernel of the code, i.e., the implicit evaluation of the right hand side function.

## Funding Data

Office of Fossil Energy (Grant No. DE-FE0026298).

## Nomenclature

*a*_{2}=velocity distribution function kurtosis

*a*,*b*=arbitrary coefficients

*d*=particle diameter

*e*=coefficient of restitution

*h*=time step size,

*h*=*t*^{n}^{+1}−*t*^{n}*L*=system size

*m*=particle mass, $m=(\pi /6)\rho sdp3$

*N*=number of time steps,

*N*=*t*_{end}/*h**p*=order of accuracy

- $R0,1*$ =
zeroth- and first-order dimensionless thermal drag functions

- Re
=_{T}thermal Reynolds number, $ReT=\rho fdpT/\mu f$

*t*=temporal dimension

*T*=granular temperature

**U**=mean velocity vector

**x**=spatial dimensions