CFD-DEM solution verification: Fixed-bed studies

Accepted Manuscript CFD-DEM solution verification: Fixed-bed studies

William D. Fullmer, Jordan Musser PII: DOI: Reference:

S0032-5910(18)30674-0 doi:10.1016/j.powtec.2018.08.044 PTEC 13629

To appear in:

Powder Technology

Received date: Revised date: Accepted date:

11 November 2017 28 June 2018 12 August 2018

Please cite this article as: William D. Fullmer, Jordan Musser , CFD-DEM solution verification: Fixed-bed studies. Ptec (2018), doi:10.1016/j.powtec.2018.08.044

This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

ACCEPTED MANUSCRIPT CFD-DEM solution verification: Fixed-bed studies

William D. Fullmera,b, Jordan Mussera National Energy Technology Laboratory, Morgantown, WV 26507, USA

b

AECOM, Morgantown, WV 26507, USA

IP

T

a

CR

Email address: [email protected] (William D. Fullmer),

US

[email protected] (Jordan Musser)

AC

CE

PT

ED

M

AN

August 16, 2018

ACCEPTED MANUSCRIPT Abstract Fixed (frozen) particle configurations are used to study CFD-DEM convergence behavior in the limit of infinite grid resolution. A

T

Gaussian filtering scheme based on solving a diffusion equation is used

IP

to interpolate particle information to the continuum grid, free of cell

CR

size limitations. A regression method is used for Richardson extrapolation of the discrete solutions to the grid-free solution. The

US

regression method is shown to predict accurate extrapolation curves in

AN

the presence of the observed oscillatory convergence behavior. The

M

extrapolation curve is used to complete the solution verification

ED

process by quantifying the numerical uncertainty present in the calculations. Several numerical and model parameters are studied

PT

including the spatial discretization method, the filter width and the

CE

mean drag law. A manageable level of discretization error is found

AC

with relative discretization errors generally less than 4%. Keywords: CFD-DEM, Solution verification, MFiX

ACCEPTED MANUSCRIPT 1 Motivation and Background Solution verification is an important aspect in the emerging field of verification, validation and uncertainty quantification (VVUQ). Solution verification deals with the

T

quantification or estimation of numerical errors that arise when a mathematical model is

IP

solved computationally. Often, discretization error is most important of the numerical

CR

errors and the present study has been devised to specifically isolate and study spatial

US

discretization error. The most general procedure is similar to informal convergence tests. First, a system response quantity of interest (SRQoI) is computed at several grid

AN

resolutions. Richard extrapolation is then used to approximate an exact or grid-free solution, i.e., what the SRQoI would be in the limit of infinite grid resolution. Finally, the

M

numerical error (signed) or numerical uncertainty (unsigned) can be quantified based

ED

primarily on the difference the computed solution on a given grid and the (approximate) exact solution. Other methods exist, e.g., numerical method refinement, error transport

PT

equation and adjoint methods, among others, [1], however, methods based on grid-

CE

refinement remain the most popular in scientific computing and are used in this work.

AC

Solution verification procedures have been primarily developed for and largely applied to single-phase flows and, more specifically, steady Reynolds-averaged Navier-Stokes (RANS)-type models. Even in single-phase flows, solution verification based on grid refinement is challenged in the case of Large-Eddy Simulations (LES). In LES, the subgrid scale (SGS) viscosity model typically depends on the grid size, so, changing the grid essentially changes the model [2].

ACCEPTED MANUSCRIPT A similar issue arises in the case of the multiphase numerical method known as computational fluid dynamics-discrete element method (CFD-DEM). The CFD-DEM method couples a traditional RANS- or LES-type CFD model for the gas-phase with a Lagrangian method for the solid particle phase where the position and velocity of each

T

particle are solved using Newton’s laws of motion and every collision is captured [3, 4].

IP

Filtering or interpolating schemes must be applied to transfer particle information, like

CR

drag and particle volume, to the continuum grid and vice versa. There are essentially two

US

“philosophical” approaches to the transfer function: grid-coupled or grid-independent.

In practice, most CFD-DEM models typically use grid-coupled methods, e.g., see Garg

AN

et al. [5] and references therein. However, the grid-coupled approach presents a unique

M

challenge to grid-refinement based solution verification methods. If the CFD grid is

ED

refined below the scale of the particle, the interpolated values are mapped to a physical space smaller than the representative particle itself. In the case of solids volume fraction,

PT

some cells will have unrealistically large values which may even exceed unity. This places a lower bound on how fine a grid-coupled interpolation-based CFD-DEM method

CE

can be refined which is approximately *  1 [6], where *  (Vtot / N xyz )1/3 / d p is the

AC

dimensionless grid size, Vtot is the volume of the computational domain, Nxyz is the total number of grid points (cells) in the domain and dp is the particle diameter. Superficially, setting limitations on the level of refinement that a solution verification study may explore is dissatisfying at best and, more critically, could be seriously problematic if the asymptotic region cannot be reached for a particular solution outside of this limitation. The concern is compounded by the fact that many studies have shown that the CFD-grid

ACCEPTED MANUSCRIPT must be refined to a level of *  2 in order to achieve grid-insensitive results, e.g., see [7, 8].

The troubling resolution limitations of grid-coupled filtering leads us to conclude that grid-independent filtering is the most appropriate method for depositing particle

IP

T

information on the Eulerian grid. A grid-independent filtering method uses a kernal

CR

specifying a weighted spatial distribution that is independent of the underlying CFD-grid. The graphical abstract qualitatively shows how the solids concentration changes as a

AN

filtering kernal is applied in the present work.

US

function of CFD grid resolution for the hypothetical bursting bubble case. A Gaussian

Solution verification of CFD-DEM is a largely unexplored topic for a variety of reasons.

M

For one, CFD-DEM is a computationally intensive method, limiting the amount of

ED

simulations practitioners are willing to undertake. Additionally, grid-coupled methods

PT

compound the challenge by i) limiting the range of viable grids and ii) effectively changing the solution method with each new grid. Finally, there is a general belief within

CE

the community, not without merit, that other uncertainties including steady drag force model, neglected interfacial force and energy transfer terms, collision model, particle

AC

shape and neglected micro-structure, etc., all contribute to larger uncertainties in CFDDEM results than the numerical grid–so long a reasonable grid is applied ( *  3 ).

In this work, we explore the topic of solution verification of CFD-DEM with gridindependent filtering. In this preliminary study, we avoid the complications associated with time-averaging by considering only fixed particle configurations (PCs), i.e., the particles are frozen in space and do not move with time. Therefore, a steady-state result is

ACCEPTED MANUSCRIPT computed and the change in the solution is attributable to the change in numerical (truncation) error associated with the change CFD-grid resolution. The numerical method briefly outlined in Sec. 2. The particle configurations are created in two ways: i) a prescribed, “fictitious” bursting bubble configuration by a particle generator, Sec. 3 and

T

ii) taken from instantaneous results of a previous, transient simulation, Sec. 4. Each case

IP

is studied on several different grids and a regression-based Richardson extrapolation

CR

method is used to calculate the grid-free solution and the numerical uncertainty [9].

US

Finally the conclusions and future outlook are briefly provided in Sec. 5.

AN

2 Methods

M

The results reported in this work were generated with the CFD-DEM model available in

ED

the open source MFiX code (https://mfix.netl.doe.gov/) version 2016.1. The governing equations and their numerical implementation are omitted for the sake of brevity and can

CE

aspects of the model.

PT

be found elsewhere [10, 11]. The remainder of this section outlines the most important

The MFiX CFD model is a finite volume implementation which uses a staggered grid and

AC

a SIMPLE-type iterative pressure-velocity coupling. Gradients are computed from face values extrapolated with a flux-limiter scheme [12], Superbee being considered the baseline method. Second-order center differencing is used for second-order derivatives. The time stepping is first order Euler, however, the time marching in this case is applied only to relax from an initial guess to a steady-state solution. In order to operate in a pseudo-steady mode, the call to the main DEM subroutine within the CFD iteration is

ACCEPTED MANUSCRIPT removed. This implementation results in three numerical tolerances: i) the linear equation solver, ii) the iterative time-step convergence, and, iii) the transient convergence. The linear equation solver tolerance is set to 10 – 8 for all equations except the Gaussian filter diffusion equation (see below) which is 10 – 10 . The value of the remaining tolerances are

T

set at 10– 5 following an extensive study which revealed that these settings produced

IP

negligible error relative to the discretization error for the range of CFD grids considered

CR

[13].

US

As mentioned in Sec. 1, we apply a grid-independent interpolation scheme, such as Gaussian filtering. While grid-independent filtering schemes do not have grid resolution

AN

limitations, there is a practical computational challenge associated with the inverse

M

relationship between the size of the grid and the number of cells contained in the stencils

ED

surrounding the particles. For distributed memory parallelism, the challenge is compounded by the relationship between the stencil and the domain decomposition

PT

strategy. Fortunately, grid-independent Gaussian filtering can be achieved by a diffusion method rather than an explicit interpolation stencil based on the distance between the

CE

CFD grid and the particle centroid, see [14, 15] for details. Briefly, data associated with

AC

the discrete particles, α, is deposited to the continuous fluid grid. In the present work, we simply use a particle-centroid deposition scheme. Then, the data is filtered by solving a diffusion equation,    f 2 , (1) 

ACCEPTED MANUSCRIPT τ is an artificial time and ν f is a viscosity associated with the diffusion length through the artificial time,





16ln 2

(2)

T

f f 

max  2f   2 , 0

IP

where δf is the Gaussian filter full width at half maximum (FWHM). In MFiX, Eq. 1 is

AN

3 Fictitious bursting bubble

US

considered as the baseline condition in this work.

CR

integrated in artificial time from τ = 0 to    f  1 . A filter width of  *f  6 is

M

The first case studied is a simple, hypothetical problem which resembles a dilute bubble

ED

about to burst at the surface of a dense bed. The problem is completely fictional and created using a particle generator script. The geometry of the domain is

PT

L*x  30, L*y  200 , and L*z  10 . The particles of diameter d p  0.01 cm only occupy the

CE

lower quarter of the domain and are separated into three regions: a dilute (1%) “freeboard” region above the surface, a “dense bed” (30%) region below the surface, and

AC

a single dilute (1%) “bubble” inside the bed near the surface. A visualization of the baseline randomization PC1 is provided in the Graphical Abstract. The x- and zdimensions are periodic BCs. In the y-direction, the inlet is a uniform mass inflow BC and the outlet is a pressure outflow BC. The inlet gas velocity is set to U 0 = 20 cm/s. The 3 4 gas-phase material properties are g  1  10 g/cm3 ,  g  2  10 g/cm-s. Gravity is set

to zero. A uniform initial condition is set with a solids concentration of 0.18 and gas

ACCEPTED MANUSCRIPT velocity of ug (t  0)  0, vg (t  0)  24.4 cm/s, and wg (t  0)  0 . Finally, we select as “baseline” conditions: i) the Superbee flux-limiter for its favorable properties with discontinuous data (Waterson and Deconinck 2007), ii) the DNS-based mean drag law of Beetstra, van der Hoef and Kuipers [16] (BVK), and, iii) a diffusion filter width (FWHM)

IP

T

of  *f  6 .

pg ( y*  y1* )  pg ( y*  y2* )

gU 02

(3)

AN

f  DP1* 2 

US

elevations averaged over the bed cross-sectional area,

CR

The SRQoI, f, is taken to be the nondimensional pressure drop between two vertical

M

We set y1*  4 and y2*  104 . The choice of y1* intentionally avoids the first row of grids

ED

cells adjacent to the inlet and does not lie along the cell centerline of any of the grids

PT

studied. Linear interpolation is used between the two rows of cells adjacent to y1* and y2* . DP1* 2 calculated on 18 different CFD-grids ranging in size from *  3.3 down to 0.16

CE

is displayed in Fig. 1. Perhaps the most important result of Fig. 1 is that the pressure drop

AC

does appear to converge to a unique solution in the limit of infinite resolution termed the grid-free or (approximate) exact solution





f 0  f *  0 .

(4)

From Fig. 1 it can be qualitatively observed that the exact solution, f 0 , occurs at DP1* 2  1159 . Another welcome result is that the maximum relative difference among all

solutions is less than 0.5% even though the grid has been refined from a modest CFD-

ACCEPTED MANUSCRIPT DEM grid size ( *  3 ) down to the boundary of DNS grid size ( *  1/ 6 ) [17]. Finally, it is observed that while the “global” convergence trend is well-behaved, the “local” or grid-to-grid behavior is oscillatory. An inherent grid dependency that arises from depositing particle data on the grid prior applying the diffusion operator may be

IP

T

responsible for the local oscillations.

CR

The oscillatory convergence behavior may likely be washed out by time-averaging in a transient simulation. However, the oscillatory pattern observed in the fixed bed studied

US

here presents a challenge to traditional (determined) methods of Richardson extrapolation. If any two or three adjacent grid points are used to extrapolate the grid-free

AN

solution, f 0 , the result will vary wildly based on the data selected and the extrapolated

M

solution may be far from the computed results. Therefore, we employ the

ED

(overdetermined) regression based method of Eça and Hoekstra [9] in which the extrapolation variables are determined by an optimization problem minimizing the least-

PT

squares of the error between the regression curve and all of the data. Here, we restrict our

CE

focus to two regression functions: power-law,

(5)

AC

fˆ  f0  c(* ) p ,

and mixed-order,

fˆ  f0  c1*  c2 (* )2 ,

(6)

here integer first- and second-order. In both cases, we use an inverse grid size weighting function,

ACCEPTED MANUSCRIPT ng

and wi  Wi /  Wi ,

Wi  1/ *i

(7)

i 1

placing a higher value on the finer grid solutions. The solution procedure used to determine the regression constants is outlined in the Supplementary Material and more

IP

T

thoroughly covered in Appendix B of [9].

CR

Fig. 1 compares the results of the regression-based Richardson extrapolation curves to the discrete solution data. In the baseline case, the power-law and mixed-order methods yield

US

very similar results, giving approximate exact solutions of 1158.80 and 1158.79,

AN

respectively. Additionally, the convergence rate, p = 1.03, given by the power-law regression is within the expected range, taken to be 0.5 to 2.1 following Eça and Hoekstra

M

[9]. The regressions can then be used to quantify the numerical error, f 0  fi , present in

ED

each calculation. Here, we adopt a concise uncertainty estimator defined by,

i  Fs  f0  fi   RE  ,

PT

(8)

 f 3

ng

ng

AC

 RE 

CE

where Fs is a factor of safety, σRE is the standard deviation of the regression,

ng

i



2

 fˆi , (9)

i 1

ng is the number of grids and fˆi is the extrapolated solution at the i-level grid. The factor of safety is taken to be:





Fs  max 1.25,3 RE / R f ,

(10)

ACCEPTED MANUSCRIPT where

Rf 

max fi  min fi (11) ng  1

T

is the data range parameter. The ratio  RE / R f quantifies the goodness of fit, with values

IP

less than unity being considered good. In the following analysis, the regression method is

US

In addition to the baseline case, we have also studied:

CR

taken as the power-law when 0.5  p  2.1 and mixed-order otherwise.

(FOU) First order upwind spatial discretization,



(MUSCL) MUSCL flux limiter [18],



(  *f  4 ) Narrower diffusion filter width,



(  *f  8 ) Wider diffusion filter width,



(  *f  16 ) Much wider diffusion filter width,



(PC2) A second randomization of particle locations in the same general

CE

PT

ED

M

AN



configuration,

(PC3) A third randomization of particle locations,



(PC3t) A small “thermal” particle velocity is imposed on the previously static

AC



PC3, 

(GB) Applying the Gidaspow-blend drag law [19, 20],



(WY) Applying the Wen and Yu drag law [21],



(HKL) Applying the Hill, Koch and Ladd drag law [22–24].

ACCEPTED MANUSCRIPT The relative discretization error (RDE),

i

/ f 0 , results for all cases listed above have been

agglomerated into Fig. 2 with the baseline case. In general, most of the DP1* 2 results and their corresponding extrapolation curves are quite similar to the baseline case in Fig. 1 [13], making it difficult to distinguish the RDEs from one case to another. However, the

IP

T

key takeaway from Fig. 2 is the relatively small error observed in all cases with RDE below 4%. In most cases, grids with *  2 have less than 2% discretization error. The

CR

one major exception is the  *f  4 case; the FOU case also has slightly higher

US

discretization error than the other cases. Although for different reasons, these two cases

AN

present a similar regression challenge of a globally nonmonotonic convergence behavior while remaining locally oscillatory. Although the RDEs are quite similar in most cases,





M

the normalization by f 0 hides the differences between the cases. For instance, the relative

ED

difference, 2 f0a  f0b / f0a  f0b , is 1.5% between the cases PC1 and PC2, 5% between

PT

cases  *f  4 and  *f  8 , 14% between cases  *f  4 and  *f  16 , 30% between cases GB and HKL and nearly 50% between cases WY and HKL. It is important to note that

CE

solution verification (verification in general) only provides an estimate of a solution at

AC

infinite resolution. Verification does not give any insight into which solution (from the different modeling options) is the best or most appropriate. This is a validation question and beyond the scope of the present work.

4 SSCP-I In this section, we wish to determine if the conclusions reached in Sec. 3 hold for more realistic yet still frozen particle configurations. To achieve this, the particle locations are

ACCEPTED MANUSCRIPT taken here from a previous transient CFD-DEM simulation of the NETL Small Scale Challenge Problem (SSCP-I) [25]. The geometry is quite similar to the previous fictitious case except the vertical BCs are now no-slip walls and the particle count is almost 20 times larger (Np = 92,948). The simulated geometry, particle and fluid properties are

T

those measured experimentally [25]. The flow condition corresponds to the 2U mf case,

IP

U 0  219 cm/s. The original simulation used the Wen and Yu [21], the Superbee flux

CR

limiter, a grid-coupled interpolation method [11] and a CFD-grid of

US

N x  N y  N z  18  96  6 or *  3.89 .

AN

Only the particle positions and (frozen) particle velocities are used for the fixed bed simulations. The fluid domain is reinitialized to a uniform condition as in the fictitious

M

case. Six instances are considered corresponding t = 5, 6, 7, 8, 9 and 10 s. Several grid

ED

resolutions ranging from *  4.67 to 0.58 are studied using the same baseline condition

PT

as before in addition to the GB drag law [19, 20]. The nondimensional pressure drop, DP1* 2 , is again taken as the SRQoI, however, the locations y1* and y2* are set to the

CE

elevations of the pressure taps in the experiment [25]. Generally, the trends of the

AC

convergence behavior are quite similar to those observed in the fictitious bursting bubble case. However DP1* 2 does not flatten out quite as much, leading to slightly larger RDEs, though still less than 4% by and large. The discretization uncertainty for a typical CFDDEM resolution of *  3 is shown in Fig. 3 superimposed on the complete transient solution of the original simulation. In this scale, it can be observed that the discretization error is less impactful than the choice of a drag law (model form uncertainty) and the natural transient variations (related to time-averaging uncertainty [26]).

ACCEPTED MANUSCRIPT 5 Conclusions This preliminary study investigated grid refinement based solution verification techniques for fixed particle configurations likely encountered instantaneously in more

T

general (transient) CFD-DEM simulations. A diffusion filter was applied for particle

IP

interpolation to remove any grid refinement restrictions. The nondimensional, area-

CR

averaged pressure drop was shown to converge to a unique solution in the limit *  0 .

US

However, significant oscillatory behavior from one grid to another necessitated the use of regression-based Richardson extrapolation techniques [9] which proved to be very

AN

effective. Overall, the discretization was found to be rather small in these fixed bed cases with typical RDEs of only a few percent. Only a few instances studied had RDEs above

M

3%, namely the coarsest grids of the FOU scheme in Fig. 3 and at times t = 9 and 10 s in

ED

Fig. 3. Qualitatively similar convergence behavior can be observed with different

PT

modeling choices, quantitatively, however, the choice of drag model, filter width, etc.

*  4 .

CE

appear to have a bigger impact on the solution than the numerical grid, at least for

AC

In closing, we suggest exercising caution using this study as a more general path forward for CFD-DEM solution verification. That is, for a given problem: freeze the simulation at several random instances, perform grid refinement, calculate the uncertainty in the SRQoIs, average the uncertainties over an ensemble of frozen instances and attribute the resulting ensemble average uncertainty to the full, time-averaged solution. This may be an accurate means to capture spatial grid discretization error alone, however, such a method cannot quantify uncertainties related to temporal discretization nor time-

ACCEPTED MANUSCRIPT averaging [26]. The results of this could may be used in the ranking of numerical errors and modeling uncertainties in future CFD-DEM studies. In general, it is somewhat encouraging that modeling uncertainties, e.g., drag, appear to be much more significant than discretization, provided a modest ( *

4 ) resolution is applied. Future work will

T

aim to study solution verification of transient CFD-DEM and the challenge of calculating

IP

numerical error in the presence of time-averaging error. Finally, this work along with Sun

CR

and Xiao [15] and Patel et al. [27] have shown that CFD-DEM solutions depend on the

US

size of the diffusion filter. Future validation studies comparing CFD-DEM to DNS or high-fidelity experimental data to determine the most appropriate filtering method would

AN

be valuable result.

ED

M

Acknowledgment

This technical effort was performed in support of the National Energy Technology

PT

Laboratorys ongoing research under the RES contract DE-FE0004000.

CE

This project was funded by the Department of Energy, National Energy Technology

AC

Laboratory, an agency of the United States Government, through a support contract with AECOM. Neither the United States Government nor any agency thereof, nor any of their employees, nor AECOM, nor any of their employees, makes any warranty, expressed or implied, or assumes any legal liability or responsibility for the accuracy, completeness, or usefulness of any information, apparatus, product, or process disclosed, or represents that its use would not infringe privately owned rights. Reference herein to any specific commercial product, process, or service by trade name, trademark, manufacturer, or

ACCEPTED MANUSCRIPT otherwise, does not necessarily constitute or imply its endorsement, recommendation, or favoring by the United States Government or any agency thereof. The views and opinions of authors expressed herein do not necessarily state or reflect those of the United States

AC

CE

PT

ED

M

AN

US

CR

IP

T

Government or any agency thereof.

ACCEPTED MANUSCRIPT References [1] W. Oberkampf, C. Roy, Verification and Validation in Scientific Computing, Cambridge University Press, Cambridge, UK, 2010.

IP

T

[2] I. B. Celik, J. Li, Assessment of numerical uncertainty for the calculations of

CR

turbulent flow over a backward-facing step, International Journal for Numerical Methods

US

in Fluids 49 (2005) 1015–1031.

[3] M. van der Hoef, M. van Sint Annaland, N. Deen, J. Kuipers, Numerical simulation

AN

of dense gas-solid fluidized beds: A multiscale modeling strategy, Annual Review of

M

Fluid Mechanics 40 (2008) 47–70.

ED

[4] W. D. Fullmer, C. M. Hrenya, The clustering instability in rapid granular and gassolid flows, Annual Review of Fluid Mechanics 49 (2017) 485–510.

PT

[5] R. Garg, C. Narayanan, D. Lakehal, S. Subramaniam, Accurate numerical estimation

CE

of interphase momentum transfer in Lagrangian-Eulerian simulations of dispersed two-

AC

phase flows, International Journal of Multiphase Flow 33 (2007) 1337–1364.

[6] Z. Peng, E. Doroodchi, C. Luo, B. Moghtaderi, Influence of void fraction calculation on fidelity of CFD-DEM simulation of gas-solid bubbling fluidized beds, AIChE Journal 60 (2014) 2000–2018.

[7] J. Capecelatro, O. Desjardins, R. O. Fox, On fluid-particle dynamics in fully developed cluster-induced turbulence, Journal of Fluid Mechanics 780 (2015) 578–635.

ACCEPTED MANUSCRIPT [8] P. Liu, C. Q. LaMarche, K. M. Kellogg, C. M. Hrenya, Fine-particle defluidization: Interaction between cohesion, Young’s modulus and static bed height, Chemical Engineering Science 145 (2016) 266–278.

[9] L. Eça, M. Hoekstra, A procedure for the estimation of the numerical uncertainty of

IP

T

cfd calculations based on grid refinement studies, Journal of Computational Physics 262

CR

(2014) 104–130.

US

[10] R. Garg, J. Galvin, T. Li, S. Pannala, Open-source MFIX-DEM software for gassolids flows: Part I-Verification studies, Powder Technology 220 (2012) 122–137.

AN

[11] R. Garg, J. Galvin, T. Li, S. Pannala, Documentation of open-source MFIX-DEM

M

software for gas-solids flows, Technical Report

ED

https://mfix.netl.doe.gov/documentation/dem_doc_2012-1.pdf, National Energy Technology Laboratory, Morgantown, WV USA, 2012.

PT

[12] N. P. Waterson, H. Deconinck, Design principles for bounded higher-order

AC

182–207.

CE

convection schemes - a unified approach, Journal of Computational Physics 224 (2007)

[13] W. D. Fullmer, J. M. Musser, An Investigation into Solution Verification for CFDDEM, Technical Report NETL-PUB-21504, DOI: 10.2172/1427020, National Energy Technology Laboratory, Morgantown, WV USA, 2017.

[14] J. Capecelatro, O. Desjardins, An euler-lagrange strategy for simulating particleladen flows, Journal of Computational Physics 238 (2013) 1–31.

ACCEPTED MANUSCRIPT [15] R. Sun, H. Xiao, Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD-DEM, International Journal of Multiphase Flow 72 (2015) 233–247.

[16] R. Beetstra, M. A. van der Hoef, J. A. M. Kuipers, Drag force of intermediate

IP

T

Reynolds number flow past mono- and bidisperse arrays of spheres, AIChE Journal 53

CR

(2007) 489–501.

US

[17] W. D. Fullmer, G. Liu, X. Yin, C. M. Hrenya, Clustering instabilities in sedimenting fluid–solid systems: critical assessment of kinetic-theory-based predictions using direct

AN

numerical simulation data, Journal of Fluid Mechanics 823 (2017) 433–469.

M

[18] B. van Leer, Towards the ultimate conservative difference scheme. V. A second-

ED

order sequel to Godunovs method, Journal of Computational Physics 32 (1979) 101–136.

[19] J. Ding, D. Gidaspow, A bubbling fluidization model using kinetic theory of

CE

PT

granular flow, AIChE Journal 36 (1990) 523–538.

[20] D. Lathouwers, J. Bellan, Modeling of dense gas–solid reactive mixtures applied to

AC

biomass pyrolysis in a fluidized bed, International Journal of Multiphase Flow 27 (2001) 2155–2187.

[21] C. Y. Wen, Y. H. Yu, Mechanics of fluidization, Chemical Engineering Progress Symposium 62 (1966) 100–111.

[22] R. J. Hill, D. L. Koch, A. J. Ladd, The first effects of fluid inertia on flows in ordered and random arrays of spheres, Journal of Fluid Mechanics 448 (2001) 213–241.

ACCEPTED MANUSCRIPT [23] R. J. Hill, D. L. Koch, A. J. Ladd, Moderate-Reynolds- number flows in ordered and random arrays of spheres, Journal of Fluid Mechanics 448 (2001) 243–278. [24] S. Benyahia, M. Syamlal, T. J. O’Brien, Extension of Hill-Koch-Ladd drag correlation over all ranges of Reynolds number and solids volume fraction, Powder

IP

T

Technology 162 (2006) 166–174.

CR

[25] B. Gopalan, M. Shahnam, R. Panday, J. Tucker, F. Shaffer, L. Shadle, J. Mei,

US

W. Rogers, C. Guenther, M. Syamlal, Measurements of pressure drop and particle velocity in a pseudo 2-D rectangular bed with Geldart Group D particles, Powder

AN

Technology 291 (2016) 299–310.

M

[26] M. Syamlal, I. Celik, S. Benyahia, Quantifying the uncertainty introduced by

ED

discretization and time-averaging in two-fluid model predictions, AIChE Journal 63 (2017) 5343–5360.

PT

[27] R. G. Patel, O. Desjardins, B. Kong, J. Capecelatro, R. O. Fox, Verification of

CE

Eulerian-Eulerian and Eulerian-Lagrangian simulations for turbulent fluid-particle flows,

AC

AIChE Journal 63 (2017) 5396–5412.

Fig. 1 (Color online.)

Convergence behavior of the bed nondimensional pressure drop

in the fictional bursting bubble case for the baseline model parameters compared to the predicted Richardson extrapolation curves for power-law (red) and mixed-order (blue).

ACCEPTED MANUSCRIPT Fig. 2 (Color online.)

Relative discretization error in the fictional bursting bubble case

for: baseline (black  [filled]), FOU (black ×), MUSCL (black +),  *f  4 (red  ),

 *f  8 (red ×),  *f  16 (red +), PC2 (green  ), PC3 (green ×), PC3t (green +), GB (blue

Fixed bed nondimensional pressure drop with discretization

CR

Fig. 3 (Color online.)

IP

T

 ), WY (green ×), and HKL (green +).

uncertainty for the six frozen particle configurations on a grid of *  3 for the BVK

US

(blue) and GB (red) drag laws compared to the complete transient of the original

AC

CE

PT

ED

M

AN

simulation (black line).

ACCEPTED MANUSCRIPT

Highlights Solution verification based on grid-refinement is considered for CFD-DEM



Fixed-bed particle arrays and grid-independent Gaussian filter interpolation is used



The pressure drop converges to a unique solution in the limit of infinite resolution



Regression-based Richardson extrapolation works well with oscillatory convergence



Discretization error of the fixed-beds studies is quite small, generally <4%

AC

CE

PT

ED

M

AN

US

CR

IP

T



Figure 1

Figure 2

Figure 3