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,79]. 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. (1)

    The HCSt is frequently used for the verification and validation of the transport coefficients derived from kinetic theory, e.g., see Refs. [1012]. 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. (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. (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. (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/ρf ≫ 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.

In the framework of the GTSH model, the HCSt specifies: a constant solids concentration, ϕ(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) = T0. Under these assumptions, the GTSH equations reduce to simply a single ordinary differential equation (ODE) governing the decay of an initial granular temperature 
dTdt=(ζ0+2γm)T
(1)
The two thermal sinks in Eq. (1) are the zeroth-order cooling rate, 
ζ0=8dpϕχ(1e2)(1+3a216)Tπ
(2)
which accounts for inelastic (e < 1) dissipation of particle collisions, and the thermal drag, 
γ=3πμfdp(R0*+ReTR1*)
(3)
which describes the viscous dissipation of particle granular temperature due to the interstitial fluid. The restitution coefficient, 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, dp, is also constant. The radial distribution function at contact, χ, accounts for volume exclusion effects of the finite sized particles. Here, we use the model of Ma and Ahmadi [20], 
χ=1+2.5ϕ+4.5094ϕ2+4.515439ϕ3(1(ϕ/ϕmax)3)0.678021
(4)
where ϕmax = 0.64356 is the solids concentration at the close packed limit. The kurtosis of the particle velocity distribution function is given by 
a2=16(1e)(12e2)8117e+30e2(1e)
(5)
and measures the deviation from an elastic, Maxwellian velocity distribution. The nondimensional concentration dependence of the thermal drag model is constituted via DNS. Here, the model of Koch and Sangani [21], 
R0*=1+3ϕ2+13564ϕlnϕ+4.6ϕχ+11.26ϕ(15.1ϕ+16.57ϕ221.77ϕ3)
(6)
with the first-order extension of Wylie et al. [22], 
R1*=0.3ϕ(1ϕ)3.6
(7)

where ReT=ρfdpT/μ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 R1 → 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 a2, 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.

Equation (1) can be cast into the form [9] 
dTdt=aT3/2bT
(8)
where 
a=8πdpϕχ(1e2)(1+3a216)+6πρfdp2R1*/m
(9)
and 
b=6πμfdpR0*/m
(10)
Equation (8) has an analytical solution of the form 
TT0=[ebt/2+abT0(ebt/21)]2
(11)

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

Numerical Method

The GTSH kinetic theory-based model is solved numerically using the open-source mfix code developed at the National Energy Technology Laboratory.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, 
Tn+1=Tn+h·f(tn+1,Tn+1)
(12)

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.

Due to the drawbacks of the available second-order time stepping method, we propose that the method be removed from the code and replaced with a more reliable, well-documented method. Although such a simple ODE problem as Eq. (11) could be solved with a very wide variety of methods, we limited our options to those which might be readily implemented in mfix in a more formal manner for all equations. The attractiveness of the current second-order scheme is that only one evaluation of the same right hand side function is required. Another method for obtaining higher order temporal accuracy with a single, implicit evaluation of the function at the new time is to use a backward differentiation formula (BDF). BDFs are a family of linear multistep methods that achieve higher order accuracy of the temporal derivatives by storing information from older time levels (i.e., 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, 
Tn+1=43Tn13Tn1+23h·f(tn+1,Tn+1)
(13)
was implemented into mfix for the granular temperature. In addition to the larger memory footprint, Eq. (13) shows another complication: BDFs also suffer from a start-up issue, i.e., at the initial 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/dp = 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=ρfdpT0/μf where T0 = T(t = 0), the density ratio, ρ* = ρp/ρf, the restitution coefficient, e, and the solids volume fraction, ϕ. Here, we use the conditions: ReT0 = 50, ρ* = 1 × 10+3, e = 0.85, and ϕ = 0.15. For completeness, we note that these conditions are obtained from dimensional properties: dp = 1 cm, ρp = 1 g/cm3, ρf = 1 × 10−3 g/cm3, μf = 2 × 10−3 g/cms−2, T0 = 1 × 10+4 cm2/s2. 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.

In addition to the analytical solution, Fig. 1 also shows several numerical solutions using the Euler time advancement scheme with decreasing time-step size. The error between the computed and exact solutions of the granular temperature will be calculated based on L1 and L2 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 
L1norm=1Nn=1N|TnT(tn)||T(tn)|
(14)
and 
L2norm=1Nn=1N(TnT(tn)T(tn))2
(15)

where Tn is the computed solution at the discrete time, tn, and T(tn) is the exact solution at the discrete time and N=tend*/h* is the number of discrete time steps.

The convergence rate of the L1 and L2 norms is shown in Fig. 2. The Euler method exhibits a first-order convergence rate as expected. The L1 and L2 norms are also reported in Table 1. The order of accuracy between two successive grid refinements, h1* and h2*, is approximated as [26] 
p=lnLi(h1*)lnLi(h2*)lnh1*lnh2*
(16)
is also reported for both norms, i = 1, 2.

The L1 and L2 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

     
  • a2 =

    velocity distribution function kurtosis

  •  
  • a, b =

    arbitrary coefficients

  •  
  • d =

    particle diameter

  •  
  • e =

    coefficient of restitution

  •  
  • h =

    time step size, h = tn+1 − tn

  •  
  • L =

    system size

  •  
  • m =

    particle mass, m=(π/6)ρsdp3

  •  
  • N =

    number of time steps, N = tend/h

  •  
  • p =

    order of accuracy

  •  
  • R0,1* =

    zeroth- and first-order dimensionless thermal drag functions

  •  
  • ReT =

    thermal Reynolds number, ReT=ρfdpT/μf

  •  
  • t =

    temporal dimension

  •  
  • T =

    granular temperature

  •  
  • U =

    mean velocity vector

  •  
  • x =

    spatial dimensions

Greek Symbols

    Greek Symbols
     
  • γ =

    thermal drag

  •  
  • ζ0 =

    zeroth-order cooling rate

  •  
  • μ =

    dynamic viscosity

  •  
  • ρ =

    density

  •  
  • ϕ =

    solids volume concentration

  •  
  • χ =

    radial distribution function at contact

Subscripts and Superscripts

    Subscripts and Superscripts
     
  • end =

    final time

  •  
  • f =

    fluid-phase

  •  
  • n =

    discrete time level

  •  
  • p =

    particle property

  •  
  • s =

    solids-phase

  •  
  • 0 =

    initial time, t = 0

  •  
  • * =

    nondimensional

References

References
1.
Haff
,
P. K.
,
1983
, “
Grain Flow as a Fluid-Mechanical Problem
,”
J. Fluid Mech.
,
134
, pp.
401
430
.
2.
Goldhirsch
,
I.
, and
Zanetti
,
G.
,
1993
, “
Clustering Instability in Dissipative Gases
,”
Phys. Rev. Lett.
,
70
(
11
), pp.
1619
1622
.
3.
Wylie
,
J. J.
, and
Koch
,
D. L.
,
2000
, “
Particle Clustering Due to Hydrodynamic Interactions
,”
Phys. Fluids
,
12
(
5
), pp.
964
970
.
4.
Garzó
,
V.
,
2005
, “
Instabilities in a Free Granular Fluid Described by the Enskog Equation
,”
Phys. Rev. E
,
72
(
2
), p. 021106.
5.
Garzó
,
V.
,
Fullmer
,
W. D.
,
Hrenya
,
C. M.
, and
Yin
,
X. L.
,
2016
, “
Transport Coefficients of Solid Particles Immersed in a Viscous Gas
,”
Phys. Rev. E
,
93
(
1
), p. 012905.
6.
Fullmer
,
W. D.
, and
Hrenya
,
C. M.
,
2017
, “
The Clustering Instability in Rapid Granular and Gas-Solid Flows
,”
Annu. Rev. Fluid Mech.
,
49
, pp. 485–510.
7.
Brito
,
R.
, and
Ernst
,
M. H.
,
1998
, “
Extension of Haff's Cooling Law in Granular Flows
,”
Europhys. Lett.
,
43
(5), pp.
497
502
.
8.
Pathak
,
S. N.
,
Jabeen
,
Z.
,
Das
,
D.
, and
Rajesh
,
R.
,
2014
, “
Energy Decay in Three-Dimensional Freely Cooling Granular Gas
,”
Phys. Rev. Lett.
,
112
(
3
), p. 038001.
9.
Yin
,
X.
,
Zenk
,
J. R.
,
Mitrano
,
P. P.
, and
Hrenya
,
C. M.
,
2013
, “
Impact of Collisional Versus Viscous Dissipation on Flow Instabilities in Gas-Solid Systems
,”
J. Fluid Mech.
,
727
, p.
R2
.
10.
Brey
,
J. J.
,
Ruiz-Montero
,
M. J.
, and
Cubero
,
D.
,
1999
, “
On the Validity of Linear Hydrodynamics for Low-Density Granular Flows Described by the Boltzmann Equation
,”
Europhys. Lett.
,
48
(
4
), pp.
359
364
.
11.
Montanero
,
J. M.
, and
Garzó
,
V.
,
2002
, “
Monte Carlo Simulation of the Homogeneous Cooling State for a Granular Mixture
,”
Granular Matter
,
4
(1), pp.
17
24
.
12.
Chamorro
,
M. G.
,
Reyes
,
F. V.
, and
Garzo
,
V.
,
2013
, “
Homogeneous Steady States in a Granular Fluid Driven by a Stochastic Bath With Friction
,”
J. Stat. Mech.
,
2013
, p. P07013.
13.
Mitrano
,
P. P.
,
Dahl
,
S. R.
,
Cromer
,
D. J.
,
Pacella
,
M. S.
, and
Hrenya
,
C. M.
,
2011
, “
Instabilities in the Homogeneous Cooling of a Granular Gas: A Quantitative Assessment of Kinetic-Theory Predictions
,”
Phys. Fluids
,
23
(
9
), p. 093303.
14.
Mitrano
,
P. P.
,
Zenk
,
J. R.
,
Benyahia
,
S.
,
Galvin
,
J. E.
,
Dahl
,
S. R.
, and
Hrenya
,
C. M.
,
2014
, “
Kinetic-Theory Predictions of Clustering Instabilities in Granular Flows: Beyond the Small-Knudsen-Number Regime
,”
J. Fluid Mech.
,
738
, pp.
R2
R12
.
15.
Brilliantov
,
N.
,
Saluena
,
C.
,
Schwager
,
T.
, and
Poschel
,
T.
,
2004
, “
Transient Structures in a Granular Gas
,”
Phys. Rev. Lett.
,
93
(
13
), p. 134301.
16.
Mitrano
,
P. P.
,
Dahl
,
S. R.
,
Hilger
,
A. M.
,
Ewasko
,
C. J.
, and
Hrenya
,
C. M.
,
2013
, “
Dual Role of Friction on Dissipation-Driven Instabilities in Granular Flows
,”
J. Fluid Mech.
,
729
, pp.
484
495
.
17.
Rubio-Largo
,
S. M.
,
Alonso-Marroquin
,
F.
,
Weinhart
,
T.
,
Luding
,
S.
, and
Hidalgo
,
R. C.
,
2016
, “
Homogeneous Cooling State of Frictionless Rod Particles
,”
Physica A
,
443
, pp.
477
485
.
18.
Liu
,
P.
,
Brown
,
T.
,
Fullmer
,
W. D.
,
Hauser
,
T.
,
Hrenya
,
C. M.
,
Grout
,
R.
, and
Sitaraman
,
H.
,
2016
, “
A Comprehensive Benchmark Suite for Simulation of Particle Laden Flows Using the Discrete Element Method With Performance Profiles From the Multiphase Flow With Interface Exchanges (MFiX) Code
,” National Renewable Energy Laboratory, Golden, CO, Technical Report No.
NREL/TP-2C00-65637
.https://www.nrel.gov/docs/fy16osti/65637.pdf
19.
Garzó
,
V.
,
Tenneti
,
S.
,
Subramaniam
,
S.
, and
Hrenya
,
C. M.
,
2012
, “
Enskog Kinetic Theory for Monodisperse Gas-Solid Flows
,”
J. Fluid Mech.
,
712
, pp.
129
168
.
20.
Ma
,
D.
, and
Ahmadi
,
G.
,
1988
, “
A Kinetic-Model for Rapid Granular Flows of Nearly Elastic Particles Including Interstitial Fluid Effects
,”
Powder Technol.
,
56
(
3
), pp.
191
207
.
21.
Koch
,
D. L.
, and
Sangani
,
A. S.
,
1999
, “
Particle Pressure and Marginal Stability Limits for a Homogeneous Monodiperse Gas-Fluidized Bed: Kinetic Theory and Numerical Simulations
,”
J. Fluid Mech.
,
400
(
1
), pp.
229
263
.
22.
Wylie
,
J. J.
,
Koch
,
D. L.
, and
Ladd
,
A. J.
,
2003
, “
Rheology of Suspensions With High Particle Inertia and Moderate Fluid Inertia
,”
J. Fluid Mech.
,
480
, pp.
95
118
.
23.
Roache
,
P. J.
,
2002
, “
Code Verification by the Method of Manufactured Solutions
,”
ASME J. Fluids Eng.
,
124
(
1
), pp.
4
10
.
24.
Choudhary
,
A.
,
Roy
,
C. J.
,
Dietiker
,
J. F.
,
Shahnam
,
M.
,
Garg
,
R.
, and
Musser
,
J.
,
2016
, “
Code Verification for Multiphase Flows Using the Method of Manufactured Solutions
,”
Int. J. Multiphase Flow
,
80
, pp.
150
163
.
25.
Drikakis
,
D.
, and
Rider
,
W.
,
2005
,
High Resolution Methods for Incompressible and Low-Speed Flows
,
1st ed.
,
Springer-Verlag
,
Berlin
.
26.
Burg
,
C. O. E.
, and
Murali
,
V. K.
,
2006
, “
The Residual Formulation of the Method of Manufactured Solutions for Computationally Efficient Solution Verification
,”
Int. J. Comput. Fluid Dyn.
,
20
(
7
), pp.
521
532
.