Document [original]

Powder Technology 436 (2024) 119447

Available online 23 January 2024

nc-nd/4.0/).

DEM and DEM-CFD modeling of systems with geometric constrictions using

a new particle location based multi-level coarse graining approach

V. Brandt

, J. Grabowski

, N. Jurtz

, M. Kraume

, H. Kruggel-Emden

Technische Universit¨

at Berlin, Chair of Mechanical Process Engineering and Solids Processing, Ernst-Reuter-Platz 1, D-10587 Berlin, Germany

Technische Universit¨

at Berlin, Chair of Chemical and Process Engineering, Marchstraße 23, D-10587 Berlin, Germany

HIGHLIGHTS GRAPHICAL ABSTRACT

•Proposing a new particle location based

MCG approach with coarsening and

refinement.

•Validation of the new MCG approach in

two hopper systems analyzing discharge

rates.

•MCG is used in coupled DEM-CFD sim-

ulations of a Wurster coater for the first

time.

•Analyzing the impact of coarsening and

refinement zone placement in MCG

simulations.

•Comparing computation times of simu-

lations without CG, CG and MCG

simulations.

ARTICLE INFO

Keywords:

Discrete element method (DEM)

Computational fluid dynamics (CFD)

Coupled DEM-CFD

Multi-level coarse graining (MCG)

ABSTRACT

A popular Lagrangian approach for modeling granular systems is the discrete element method (DEM). Currently,

the simulation of industrial granular systems often exceeds modern computer capacities. Therefore, the coarse

graining (CG) approach bundles multiple particles into grains, which are less computationally demanding due to

the reduced number of entities considered. In systems, where the maximum grain size is limited by local ge-

ometry constrictions, uniform grain sizes lead to inefficient computations since far from the constrictions larger

grain sizes could be reasonably used. Multi-level coarse graining (MCG) addresses this by using multiple grain

sizes, which are adapted to the geometric constrictions. In this study, geometrical non-uniform hoppers and a

Wurster coater system were modeled using a new MCG approach. Results from literature regarding hopper

discharge were matched. For the Wurster coater system MCG outperforms conventional CG by providing a higher

resolution of the granular material close to the geometric constriction.

1. Introduction

Granular materials play a key role in a wide range of applications in

the chemical and pharmaceutical industry as well as in agriculture, food

and raw material processing [1]. During processing, granular materials

are either only mechanically or gravitationally agitated or they are often

also fluidized for transport as well as drying, coating or agglomerating

particles. When only being agitated mechanically or by gravity the

* Corresponding author.

E-mail address: [email protected] (V. Brandt).

Contents lists available at ScienceDirect

Powder Technology

journal homepage: www.journals.elsevier.com/powder-technology

https://doi.org/10.1016/j.powtec.2024.119447

Received 23 October 2023; Received in revised form 11 January 2024; Accepted 18 January 2024

Powder Technology 436 (2024) 119447

contact forces between particles or between particles and externals,

which result from particle-particle and particle-externals collisions are

relevant for the motion and behavior of the granular material. With a

fluid phase being additionally involved during fluidization, besides

contact forces the particle-fluid interaction forces can play a major role

in the motion and behavior of the granular material.

Today, there exist several approaches to model the complex physical

behavior of granular materials under mechanical or gravitational

agitation or when being fluidized. A popular Lagrangian approach for

modeling granular materials is the discrete element method (DEM), first

introduced by Cundall and Strack [2]. It is capable of tracking the tra-

jectory of each particle in a granular system and resolves the interaction

forces between particles or between particles and externals at a micro-

scopic level by solving contact force models. In contrast to other ap-

proaches that treat the granular material as a continuum [3], the DEM is

able to model non-spherical particles [4,5] and polydispersity. When a

fluid is additionally involved besides the granular material, direct nu-

merical simulations (DNS), coupled with the DEM to solve the particle

motion, can be performed to obtain highly resolved information about

the fluid flow and its interaction with the granular material at a

microscopic level [6–9]. In DNS, no additional closures are required to

model the interaction forces between the fluid and the granular phase.

However, a drawback is that this method is computationally highly

intensive and it is therefore reserved for the development of closures or

for performing reference simulations to benchmark other modeling ap-

proaches [10]. A less computationally intensive method is to couple the

DEM to the unresolved computational fluid dynamics (CFD) referred to

as DEM-CFD. In the DEM-CFD, interaction forces between the granular

and fluid phase are computed by drag force closures, while other

particle-fluid forces can also be described. A well-known drag closure is

the Gidaspow drag correlation [11], a combination of the Ergun drag

correlation [12] for dense regions and the Wen-Yu drag correlation [13]

for dilute regions of the modeled system.

Despite advances in CPU and GPU hardware and in parallelization

techniques, the application range of DEM and DEM-CFD simulations is

limited by its computational cost. Currently, the computation of gran-

ular systems with more than 106 particles quickly saturates the

computational capacity of modern computers, even with highly paral-

lelized and optimized DEM and DEM-CFD software on multi-core pro-

cessors [14,15]. Therefore, industrial-scale systems with up to

108−1012 particles easily exceed the limits of feasibility. To overcome

this problem, multiple particles can be bundled into large-sized parti-

cles, so-called grains. The resulting coarse grained DEM (CG-DEM) is less

computationally demanding due to the reduced number of entities

considered. However, scaling rules must be applied to the interaction

forces to obtain accurate results. There are many different CG ap-

proaches [16–21] to scale interaction forces, which were reviewed by Di

Renzo et al. [22]. Detailed benchmarking studies of CG models in

different systems have recently been published by other authors

[23–25]. Recently [23], CG scaling laws for contact forces and hydro-

dynamic forces have been benchmarked separately. For the contact

forces, the simplified dissipation scaling model of Benyahia and Galvin

[18] was identified as the most promising contact scaling approach,

which uses an additional scaling of the coefficient of restitution. For

further details on applicable approaches, see subsection 3.1. According

to Brandt et al. [23], the scaling of hydrodynamic forces is described

most accurately by the volumed-based scaling approach of Sakai and

Koshizuka [20] (see subsection 3.2 for details).

CG-DEM is typically limited to a maximum CG factor l (ratio of grain

to particle diameter) of l=3…5 [20,26–28] with some exceptions

[16,29]. As pointed out by Brandt et al. [23], the considered granular

system, the output parameters and the current system-to-grain size ratio

have a strong influence on an appropriate CG limit. Only uniform

cylindrical systems were analyzed by the authors, where the maximum

grain size is not restricted by local geometry constrictions. Hence, the

desired level of accuracy and the intended speed-up determines the CG

limit. In many industrial applications such as hoppers [30], coating

apparatus [31,32], circulating fluidized beds [33] and other systems

with geometric constrictions, the physical behavior is strongly depen-

dent on the size of the particles and in CG simulations on the size of the

considered grains. For a conical hopper, the maximum grain size is

limited by the size of the orifice to obtain reasonable results in a simu-

lation. At the same time, higher grain sizes are acceptable at a sufficient

distance from the orifice. Therefore, using a uniform grain size

throughout the system often results in an inefficient computation. Multi-

level coarse graining (MCG) addresses this issue by using more than one

grain size, defined by different CG factors in the simulation domain. In

MCG, the CG factor represents the local resolution of the granular ma-

terial, which is adapted close to geometric constrictions of a treated

system. The MCG approach was first considered by Queteschiner et al.

[34] and later by De et al. [15]. In the work of the former authors, the

refinement of the granular material was realized by monitoring cell

averaged properties such as the contact and kinetic stress and the

granular pressure on representative volume elements. Three

proportional-integral (PI) controllers per cell were used to monitor and

control the above mentioned properties before and after the refinement

of the granular material. A comparatively easier to implement MCG

approach was introduced by De et al. [15], where the refinement of the

granular material is determined only by its position. To date, all MCG

studies have dealt solely with hopper discharge.

The goal of this study therefore was to develop and validate a new

MCG approach, which differs from existing approaches in terms of the

coarsening and refinement and apply it to DEM-CFD simulations of

fluidized particle systems for the first time. The procedure was as fol-

lows: First, two hopper systems were analyzed to validate the new MCG

approach with two MCG studies [15,34] from literature. The discharge

rates were used as validation parameter. Thereafter, a Wurster coater

system was analyzed as a representative of a fluidized particle system

with geometry constrictions that make the granular flow dependent on

the size of particles or grains. The novelty of this work consists in the

application of MCG to a system described by DEM-CFD, where the res-

olution of the granular material is adapted close to geometric constric-

tions of the system. To achieve this, the granular entities were bundled

and refined multiple times in the circulating flow pattern of the Wurster

coater. For the Wurster coater, simulations without CG were compared

with CG and MCG simulations with respect to key evaluation parame-

ters. In addition, the influence of different CG contact force scaling

models on the results of MCG simulations was analyzed.

The paper is structured as follows: In section 2, the DEM-CFD method

is briefly outlined. In section 3, the CG approaches and the new MCG

approach used in the context of DEM and DEM-CFD are explained.

Subsequently, the simulation conditions and the numerical setup are

explained in section 4. Section 5 presents the results of the validation

simulations of our new MCG approach for hopper discharge. The results

of the Wurster coater using the new MCG approach under application of

a DEM-CFD are analyzed in section 6, followed by conclusions drawn in

section 7.

2. Method description

The governing equations of the fluid and solid phase were solved

separately in the coupled DEM-CFD method. To realize a regular

coupling, both sets of governing equations were solved with different

time steps one after the other. For the solid phase solution, at least drag

and pressure gradient forces must be considered (see subsection 2.1).

When solving the equations of the fluid phase, the local solids fraction

V. Brandt et al.

Powder Technology 436 (2024) 119447

and momentum exchange between fluid and solid phase must be taken

into account (see subsection 2.2).

2.1. Particle motion

The solid phase was solved using the DEM approach, where the

translational and rotational motion of each particle can be tracked. The

motion of each individual particle was determined by Newton's and

Euler's equations as follows:

d2x

→p

dt2=F

→c

p+F

→h

p+mpg

→,(1)

→=M

→c

p+M

→r

p+M

→h

p,(2)

with the particle mass mp, the particle acceleration d2x

→p/dt2, the con-

tact force F

→c

p consisting of normal and tangential components, the total

hydrodynamic force F

→h

p consisting of pressure gradient and drag forces,

and the gravitational force mpg

→, the angular acceleration

→, the

moment of inertia Ip, and the external moments M

→c

p, M

→r

p and M

→h

p which

result out of contact forces, rolling resistance or hydrodynamic forces.

Note that external moments Mh

p due to hydrodynamic forces were

neglected in the following, since only drag and pressure gradient forces

are considered to contribute to particle-fluid interactions, which do not

exert moments.

In this study, the linear-spring-dashpot model [2] was used as con-

tact force model. It determines the contact forces between particles or

particles and walls as a superposition of normal and tangential acting

forces as follows:

→n=knδ n

→+γnv

→n

rel,(3)

→t=min(ktξ

→t+γtv

→t

rel,

cF

→n)t

→,(4)

with normal and tangential spring stiffness kn and kt, virtual overlap

δ, normal and tangential vectors n

→and t

→, normal and tangential

damping coefficients γn and γt, relative normal and tangential velocities

at the contact point v

→n

rel and v

→t

rel, tangential displacement ξ

→t and the

coefficient of Coulomb friction

c. The damping coefficients are

dependent on the coefficient of restitution according to [35]. When

rolling resistance was considered, it follows the model of [36]. The

moment M

→r

p, resulting from the rolling resistance can be calculated with

the rolling resistance coefficient

r, the normal force F

→n, the rolling

radius Rr and the relative angular velocity

→rel using

→r

p= −

r,pF

→nRr

→rel



→rel

.(5)

The combined drag and pressure gradient force acting on a single

particle can be calculated as

→h

p=Vpβ

f(1−

f)(u

→f−v

→p),(6)

where

f is the void fraction, u

→f and v

→p are the fluid and particle ve-

locities and β is the interphase momentum transfer coefficient [37]. The

Gidaspow correlation [11] is a widely used voidage dependent model to

calculate the interphase momentum transfer coefficient, which is a

combination of the Ergun [12] and the Wen-Yu correlation [13] to

distinguish dense and dilute regions in particle-fluid system:

β=⎧

⎪

⎨

⎪

⎩

150 (1−

f)2

p+1.75(1−

dpu

→f−v

→p,

f≤0.8

4Cd

f(1−

fu

→f−v

→p

f−2.65,

f>0.8

(7)

with the fluid viscosity

f, the fluid density

f and the particle diameter

dp, respectively. The drag coefficient Cd is given by

Cd=⎧

⎪

⎨

⎪

⎩

Rep(1+0.15⋅Re0.687

p),Rep≤1000

0.44,Rep>1000.

(8)

To calculate the drag coefficient, the particle Reynolds number Rep is

required, which can be determined as

Rep=u

→f−v

→p

fdp

.(9)

2.2. Fluid flow

The fluid flow was described by the CFD in an Eulerian formulation

by considering a volume averaging approach. The continuity and mo-

mentum equations were solved to describe the motion of the incom-

pressible fluid:

∂ε

∂

t+∇⋅(

→f)=0,(10)

∂

(

→f)

∂

t+∇⋅(

→fu

→f)= −

f∇p+∇⋅(

)+

→

β(¯

→p−u

→f).

(11)

Here p is the pressure,

 is the fluid viscous stress tensor, which is

related to the effective viscosity determined from the standard k−ϵ

turbulence model as in [32]. The solid medium has no effect on the

turbulence generation or dissipation rates of the fluid phase. The CFD

cell averaged interphase momentum transfer coefficient ¯

β and the CFD

cell averaged particle velocity ¯

→p are calculated as follows: ¯

β=∑βi/Np,

→p=∑vi

→/Np, where Np is the number of particles in a CFD cell. If DEM

and CFD cells are of different shape and size, it is necessary to use a two-

grid approach as in [38].

3. Coarse graining

Bundling a pre-set number of particles into a spherical solid grain is

the main principle of CG. From the solid phase perspective, the

Lagrangian DEM approach thus requires a reduced number of entities to

be tracked, resulting in a significant reduction in computational cost.

Ideally, CG simulations yield identical results to simulations without CG.

To maintain a consistent nomenclature, particles are referred to in the

following by the index p. The index g is used to refer to grains, which are

solids consisting of bundled particles. The ratio of the grain diameter to

particle diameter, called the CG factor l=dg/dp, is a fundamental

parameter in the context of CG. The CG number nCG represents the

number of particles per grain and is also related to the CG factor as nCG =

np/ng=l3. In general, CG is realized with respect to the conservation of

total mass Mtot =∑ngmg=∑npmp and the total volume Vtot =∑ngVg=

∑ngVp, which is linked to

p and mg=nCGmp. It is assumed that all

represented particles inside a grain have the same translational velocity

as a grain, such that vg

→=vp

→. The equivalence of the rotating kinetic

V. Brandt et al.

Powder Technology 436 (2024) 119447

energy requires that l⋅

→=

→, l⋅

→=

→and l⋅ φ

→g=φ

→p.

Applying CG in the context of a particle-fluid systems requires

scaling of the contact forces (see subsection 3.1), which can be sub-

divided into models ensuring constant absolute overlap (see sub-

subsection 3.1.1) and constant relative overlap (see subsubsection

3.1.2). Additionally, particle-fluid forces are scaled (see subsection 3.2).

Our approach to MCG is given in subsection 3.3 and the following

subsubsections 3.3.1 and 3.3.2.

3.1. Contact force scaling

3.1.1. Constant absolute overlap models

The constant absolute overlap model [20] (abbreviated as l3) as-

sumes that the total energy of the bundled particles and the represen-

tative grain is identical (see Di Renzo et al. [22] for details). The contact

between two or more grains is modeled as a simultaneous contact of all

particles within the colliding grains. Consequently, the contact force

must be l3- times larger than when two particles collide, resulting in the

normal and tangential components of the contact force as follows:

→n

g=l3F

→n

p=l3(kn

pδgn

→g+γn

→n

rel,g),(12)

→t

g=l3F

→t

p=min(l3(kt

pξ

→t

g+γt

pv

→t

rel,g),

c,pF

→n

g)t

→.(13)

Eqs. (12) and (13) can be interpreted as scaling laws for the normal

and tangential contact parameters of the linear-spring-dashpot model,

leading to kn

g=l3kn

p, kt

g=l3kt

p, γn

g=l3γn

p, γt

g=l3γt

p and

c,p=

c,g.

Grains maintain the same maximum absolute overlap during colli-

sions as particles due to the applied spring stiffness scaling. The ratio of

normal to tangential spring stiffness and damping coefficients remains

the same when the l3- contact force scaling is applied. Since the pro-

portion of mass, spring stiffness and damping also remains the same, the

collision time of grains tg is equal to the collision time of particles tp

according to

tg=



meff,g−(γn

2meff,g)2

√−1



l3kn

l3meff,p−(l3γn

2l3meff,p)2

√

−1

=tp,(14)

where meff,g and meff,p are the effective masses of the colliding grains

and particles, respectively. The restitution coefficient of the grains eg

remains unchanged compared to the restitution coefficient of the par-

ticles ep because the damping coefficient γn

p and the effective mass meff,p

are both scaled with l3, as

eg=e−γn

2meff,gtg=e

−l3γn

2l3meff,p

tp=ep.(15)

In the original CG model of Sakai and Koshizuka [20], no friction

scaling is applied, resulting in

c,p=

c,g. An additive friction scaling has

been proposed by Cai and Zhao [19] to preserve the dissipation during

sliding according to

l

√.(16)

Benyahia and Galvin [18] extended the constant absolute overlap

model in order to include the additive amount of dissipation caused by

intra-grain collisions. According to them, the coefficient of restitution is

scaled by the CG number, referred to as simplified dissipation scaling

(SDS), as

eg=e

nCG

√

pfor

nCG(γn

p)2

2kn

pmeff,p

≪1.(17)

Normal and tangential spring stiffness components are scaled by l3.

Tangential damping is scaled such that the ratio

=γn

g/γt

p of normal to

tangential damping coefficients remains constant when CG is applied.

No friction scaling is applied in the original model.

A modification of the model of Benyahia and Galvin combined with

the friction scaling of Cao and Zhao [19] was identified by Brandt et al.

[23] as a promising scaling rule in mechanical agitated and fluidized

particle systems as it preserve the dissipation during sliding and includes

the additive amount of dissipation caused by intra-grain collisions.

Therefore, it was also tested in this study, where it is called simplified

dissipation scaling with friction scaling (SDS-fric).

3.1.2. Constant relative overlap models

Motivated by the idea of scaling spring stiffness and damping pa-

rameters such that the maximum relative overlap during collisions be-

tween particles and grains is equal, Bierwisch et al. [39] and later Radl

et al. [21] proposed constant relative overlap models (abbreviated as l2)

for CG. The l2-dependence results in smaller spring stiffness and

damping parameters than the above discussed l3- scaling given by kn

lkn

p, kt

g=lkt

p, γn

g=l2γn

p and γt

g=l2γt

p while

c,p=

c,g. Due to the changed

relationship between spring stiffness, mass and damping, the collision

time scales linearly with the CG factor tg=ltp. This results in the same

coefficient of restitution for both particles and grains, eg=ep. A scaled

collision time for grains increases the integration time step size for CG

simulations, which leads to additional savings in computational time.

3.2. Hydrodynamic force scaling

Following section 2.1, the hydrodynamic force acting on a single

particle F

→h

p consists of the pressure gradient force F

→∇p

p and drag force

→d

p. For a grain, the pressure gradient force is scaled by volume. The

force acting on the grain is assumed to be the same as the force acting on

the bundled particles inside the grain [22]:

→∇p

g=l3Vp∇p=Vg∇p(18)

There is a large consensus in the literature (see e.g. [22,23,37]) to

scale the drag force in a similar way, such as

→d

g=l3Vpβp

(1−

f)(u

→f−v

→p)=Vgβp

(1−

f)(u

→f−v

→g).(19)

As a result, the combined hydrodynamic force acting on a grain is

then obtained as

→h

g=l3Vpβp

f(1−

f)(u

→f−v

→p)=Vgβp

f(1−

f)(u

→f−v

→g).(20)

3.3. Multi-level coarse graining

Two essential steps are performed as part of the multi-level coarse

graining (MCG). In the coarsening step, particles / grains are bundled

together. In the refinement step, the grains are refined to particles /

grains with a lower CG factor as already proposed by De et al. [15]. In

our MCG approach, the coarsening and refinement of particles / grains is

only dependent on the particle / grain positions similar to the MCG

approach of De et al. [15]. The coarsening and refinement zones are

predefined as 3-D domains in which the coarsening and refinement steps

are carried out. The process of locating the coarsening and refinement

zones requires detailed knowledge of the flow pattern of the granular

material and geometry constrictions. Careful placement of coarsening

and refinement zones is very important as it potentially affects the re-

sults and computation time of the MCG simulations (see subsection 5.2

and section 6). Both the coarsening and refinement steps are performed

periodically after a predefined number of DEM time steps, depending on

the particle dynamics of the treated system. When MCG is applied, the

contact and hydrodynamic forces are scaled according to subsection 3.1

V. Brandt et al.

Powder Technology 436 (2024) 119447

and subsection 3.2. In MCG simulations, grains with different CG factors

or grains and uncoarsened particles may collide. Commonly used con-

tact force scaling laws only provide the scaling of contact parameters for

a contact of same coarsened particles / grains. In the DEM, contact pa-

rameters must be however defined for all possible contact partners. In

the case of a contact between particles / grains of different CG factors,

the mean values of the respective contact parameters (spring stiffness,

damping and friction) for contacts with same-coarsened particles /

grains are used. Details on the coarsening and refinement step required

for MCG are explained as part of the following subsubsections 3.3.1 and

3.3.2.

3.3.1. Coarsening step

The coarsening step is performed in the coarsening zone and bundles

either particles or grains to obtain a higher CG factor. The bundling is

always done pairwise and can take place in multiple steps. The bundling

results in a reduced number of entities and a lower resolution of the

granular material. Each coarsening zone has its predefined target CG

factor and location. The coarsening zone is defined as a 3-D domain,

because it must contain multiple particles / grains to bundle them

together. If grains already have the target CG factor, they are excluded

from the coarsening step. When particles / grains to be bundled together

enter the coarsening zone, potential coarsening partners must be found.

Therefore, a nearest neighbor search is performed.

The nearest neighbor search considers only particles / grains with

the same CG factor. The search is done based on cells, which are larger

than the particles / grains and are provided as part of the efficient

realization of the contact detection among particles / grains and parti-

cles / grains and walls in the DEM. If a particle / grain is already

designated as a nearest neighbor, it will be neglected in the neighbor

search of other particles / grains. Therefore, it may occur that a particle

/ grain nearest neighbor is already assigned to another particle / grain.

In this case, the nearest available neighbor is determined. Nearest

neighbors are bundled together to create a new grain with a CG factor of

l=

2⋅nCG

√. Here, nCG is the CG number of the particles / grains to be

bundled. If a pair of nearest neighbors is bundled together while one or

both of them are in contact with surrounding particles, their contacts are

removed. When performing a nearest neighbor search, it may happen

that some particles / grains do not find a dedicated nearest neighbor. For

example, if there is an odd number of particles / grains with an identical

CG factor present in the coarsening zone at least one particle / grain will

not find an available neighbor to be bundled with. These particles /

grains will not be coarsened in the current coarsening step, but will be

checked again in following steps.

The coarsening procedure, which bundles particles with l=1 into a

grain with l=2 in three intermediate steps, each performed in separate

DEM time steps, using the nearest neighbor search, is discussed below as

an example. For a grain with a CG factor of l=2, the CG number is nCG =

8. Hence, this number of particles must be bundled together in the

coarsening step. Starting with 8 particles of l=1 as shown in Fig. 1, a

nearest neighbor search is performed. The nearest neighbor particles are

then bundled together, resulting in grains with an intermediate CG

factor of l=

√to preserve the particle mass. Thereafter, a nearest

neighbor search is performed on each of the created grains again. The

coarsening of 4 grains with l=

√results in two grains with an inter-

mediate CG factor of l=

√, which must be bundled again to complete

the coarsening step and reach the target CG factor of l=2.

When a grain is created, the initial translational velocity of the grain

is determined by averaging the translational velocities of the particles /

grains to be bundled. To obtain the initial rotational velocity of the

created grain

→it requires a scaling of the averaged rotational veloc-

ities of the particles / grains to be bundled according to

→=∑

i=1

p,i

→

2⋅ 

√,

→=∑

i=1

g,i

→

2⋅ 

√,(21)

with

p,i

→being the rotational velocity of the particles and

g,i

→the rota-

tional velocity of the grains to be bundled. Since the diameter of the

created grain increases by 

√, the average rotational velocity of the

particles / grains to be bundled is inversely scaled by 1/

√to maintain

the tangential velocity component resulting out of rotational motion at

the surface of the created grain.

To determine the new position of the created grain, it is intuitive to

average the positions of the particles / grains to be bundled as proposed

by De et al. [15]. However, this method can have inefficiencies, because

the averaged position of the particles / grains to be bundled can result in

the placement of grains with a large overlap with other neighboring

particles / grains or internal walls. Large overlaps are a problem in the

force based DEM. They can be addressed by either not placing the

created grain with its bundled particles / grains at all and looking for

other combinations to bundle particles / grains or by scaling down the

grain size and allowing it to grow to its regular size in consecutive DEM

time steps under a contact force triggered repositioning (see [15] or

[40]).

A simpler way to place the created grain is to replace the position of

one of the particles / grains to be bundled after removing both. This

ensures that the created grain is not strongly overlapping with neigh-

boring particles / grains or walls. However, some overlap with neigh-

boring particles / grains or walls can still occur, because the placed grain

has a larger diameter than the bundled particles / grains. To address this,

the created grain is placed with a scaled down diameter of ds. If ds is

chosen to be 1/

√of the regular diameter dg of the created grain, it is

ensured that no artificial overlap with surrounding particles / grains

prevails.

When created grains are scaled down as part of the coarsening step, it

is important that they grow afterwards to their regular diameter to

restore mass conservation. In the work of De et al. [15], the growth of

the diameter of scaled down grains is predefined by step, linear or

exponential time dependent functions. In our MCG approach, the grain

diameter can however grow instantaneously up to an predefined overlap

threshold with other particles, grains or walls as described by Bußmann

et al. [40]. When considering absolute values, larger overlaps are more

Fig. 1. Example of the coarsening step including its most delicate aspects,

starting with 8 particles with l=1 on the left, which are bundled in three in-

termediate steps to a grain with l=2 on the right.

V. Brandt et al.

Powder Technology 436 (2024) 119447

critical for small grains. Therefore, the overlap threshold is defined as a

percentage of the actual scaled down diameter ds. If the maximum

overlap threshold is reached, which predominantly occurs in denser

regions of a granular system, the grains do not grow further in that DEM

time step. Due to the resulting overlap, the surrounding particles /

grains are reorienting and creating new space for the created grain to

grow further in the next DEM time steps.

Using an overlap threshold instead of a predefined growth function

such as done by De et al. [15] has several advantages. If there is no

overlap with surrounding particles / grains or walls, the created grain

can grow instantaneously to its regular diameter and mass conservation

is at once ensured. In addition, predefined growth of the grain diameter

when not chosen carefully can lead to artificially high overlaps in denser

regions of a granular system, where created grains have no immediate

space to grow due to surrounding particles / grains or externals.

3.3.2. Refinement step

The refinement step is the opposite of the coarsening step. Thus, it

refines grains into particles / grains with a lower CG factor, when

entering the refinement zone. The refinement zone is placed in the

system to obtain a finer resolution of the granular material where

needed. Each refinement zone has its predefined target CG factor and

location. Whereas De et al. [15] refine grains that pass a predefined

height level, in this work, analogously to the coarsening step, a 3-D

domain is used to define the refinement zone. The reason for this is

that in more complex systems such as a Wurster coater, the refinement

zones are shaped according to the local geometry constrictions and

particles / grains can enter these zones from different directions. When

particles / grains enter the refinement zone and already exhibit the

target CG factor, they are excluded from the refinement step.

The refinement step is performed in a single step. In Fig. 2, the

refinement of a grain with l=2 into 8 particles with l=1 is shown. The

particles are placed in a cubic pattern around the center of the grain to

be replaced to ensure that the center of mass is unchanged. In this study,

the refinement step is used to refine grains from a CG factor of l=2 to

l=1, l=4 to l=2, l=6 to l=3, l=8 to l=4 and l=16 to l=8. This

means respectively, that grains are always refined into a cubic pattern of

8 particles / grains. Therefore, the pattern in which particles / grains are

placed will always look like shown in Fig. 2, even when for example a

grain with a coarse grain with l=16 is refined to 8 grains with l=8. The

procedure used here for placing particles / grains of lower CG factor can

be easily extended to arbitrary changes in CG factor if required. How-

ever, details on this are not given here.

All generated particles / grains are placed with the same trans-

lational velocity as the grain being replaced. Scaling rules inverse to the

coarsening step are applied to the rotational velocity of the created

particles / grains. The rotational velocity is set to l⋅

→with l being the

CG factor and

→being the rotational velocity of the grain to be

replaced. The placement of the particles / grains of lower CG factor in

total takes up more volume than the grain to be replaced needed.

Therefore, the placed particles / grains are positioned with a scaled

down diameter, which is adjusted ensuring a predefined overlap

threshold in the following DEM time steps. Thereby it is followed a

procedure as already described for the coarsening step in subsubsection

3.3.1. If particles / grains are placed as shown in Fig. 2 then a scaled

down diameter ds which is 2 −

√times smaller than the regular

diameter is chosen to ensure that particles / grain can initially be placed

without overlapping with surrounding particles / grains or walls.

4. Simulation conditions and numerical setups

The aim of this work was to apply the proposed MCG in different

systems and to evaluate its results in comparison to simulations without

CG and CG simulations with a fixed CG factor, which includes results

from literature.

For simulations without CG and CG as well as MCG simulations the

motion of the particles / grains was solved with a parallelized DEM code

used in previous studies (see e.g. [23,40]). In the DEM, the collision time

tp in case of simulations without CG is calculated using eq. (14). A

collision was then resolved in 25 DEM time steps. The DEM time step

applied can vary between simulations without CG and CG and MCG

simulations when, e.g., the l2- contact force scaling model is applied,

which scales the collision time of grains. Linear-spring-dashpot models

(see section 2) were utilized as contact force models for all considered

systems and simulations (no CG, CG and MCG). So called DEM cells were

used for the nearest neighbor search in the MCG and were also relevant

for the void fraction calculation for the DEM-CFD. For DEM cells, a fixed

cell to particle size ratio dx/dp of 2 was chosen for all simulations

without CG. When CG simulations were performed, the DEM cell size

was scaled to maintain the ratio dx/dg of 2. For MCG simulations dx/dg is

also 2, while dg is the largest grain size considered. When rolling resis-

tance was considered, the rolling resistance coefficient remained un-

changed when CG or MCG simulations were performed in accordance to

[34].

First, the new MCG was studied for dense particle flow in hoppers.

Here, the local maximum grain size is limited by the size of the orifice,

which makes hoppers predestined for the application of MCG. We

validated our MCG approach with the results of both previous MCG

studies performed for a rectangular [34] and a conical hopper [15],

where only the refinement step is required to simulate hopper discharge.

Second, a Wurster coater was analyzed. When applying MCG to the

simulation of a Wurster coater, both coarsening and refinement steps

were required to model the flow of granular material with changing

grain sizes. Thereby, a refinement was applied close to the geometric

constrictions of the system. More details on the system geometries and

simulation parameters for both hoppers and Wurster coater are given

below.

4.1. Rectangular hopper

In order to validate the proposed MCG approach with the results of

Queteschiner et al. [34] for discharge from a rectangular hopper, the

same geometry and simulation parameters were used as in the original

work. The geometry of the rectangular hopper (see [41] for more details

on the geometry) is shown in Fig. 3. The system consists of a box with a

height of 1.5 m and a cross section with an edge length of 0.04 m. The

particle outlet is a part of the bottom cross section with dimensions of

0.04 m ×0.01 m.

In our study, we compared the results of a simulation without CG as a

reference simulation with CG simulations using l=2 and l=4, with a

MCG simulation using l=4/2/1, and with results of Queteschiner et al.

[34]. The notation l=4/2/1 for the MCG means that three grain sizes

exist in parallel and l=4 grains are refined to grains with l=2, which

are then refined to particles with l=1. In the simulations of hopper

discharge, only the refinement step (see subsubsection 3.3.2) of the MCG

Fig. 2. Example of the refinement step including its most delicate aspects,

starting with a grain of l=2 on the left, which is refined into 8 particles with

l=1 on the right.

V. Brandt et al.

Powder Technology 436 (2024) 119447

was carried out. The refinement zone from l=4 to l=2 was placed at a

height of 0.02 m and the refinement zone from l=2 to l=1 was placed

at a height of 0.01 m (shown in Fig. 6 in accordance with [34]). As

contact force scaling model, the l2- scaling was used to achieve com-

parable results with the work of Queteschiner et al. [34]. Table 1 sum-

marizes the simulation conditions for the reference simulation without

CG, CG and MCG simulations.

Table 2 shows the contact parameters used in the DEM for the

simulation without CG. Except for the coefficient of friction, the contact

parameters for particle-particle collisions were identical to those for

particle-wall collisions. A DEM time step of 3.2⋅10−6 s was used in the

simulations without CG that is scaled by the CG factor l when using the

l2- contact force scaling according to subsubsection 3.1.2. In MCG sim-

ulations with l=4/2/1, the same DEM time step as in the simulation

without CG was used since particles with l=1 were present. Due to the

dense bed forming in the simulations, the overlap threshold for the

particle / grain growth was set to 1.5% of the diameter of the scaled

down particles / grains.

4.2. Conical hopper

The conical hopper considered in the work of De et al. [15] was

modeled with the same simulation parameters as in the original work to

validate the proposed MCG approach. The geometry of the conical

hopper is shown in Fig. 4. The hopper consists of a cylindrical part with a

length of 0.2 m and a diameter of 0.04 m. Below, the conical part of the

hopper begins, which has an opening angle of 75◦and a length of

0.0672 m. The orifice has a diameter of 0.004 m and a length of 0.025 m.

Since the simulation of 8.48⋅106 particles with a diameter of dp=

160⋅10−6m took more than six months to compute by De et al. [15], the

CG simulation with l=2 was considered as the reference simulation by

them. The same approach was used in this study. The reference simu-

lation with l=2 was used as comparison to CG simulations with l=4

and l=8 and to MCG simulations using l=8/4/21 and l=16/8/4 and

results of De et al. [15]. To predict the discharge rate in MCG

Fig. 3. Geometry of the rectangular hopper used in [34], including the geo-

metric properties of the system and its discharge zone in light gray.

Table 1

Simulation conditions for discharge from a rectangular hopper according to

[34].

Reference simulation (l=1):

Number of particles 960,000

Particle density in kg/m

2500

Particle diameter in m 0.001

Simulation time tsim in s 10

Contact force model Linear-spring-dashpot

CG and MCG simulations:

CG simulations l=2, l=4

MCG simulation l=4/2/1

Contact force scaling l2 [38]

Table 2

Contact parameters for particle-particle and particle-wall in-

teractions for modeling discharge from the rectangular hopper.

Normal spring stiffness kn

p in N/m 1000

Tangential spring stiffness kt

p in N/m 1000

Particle friction coefficient

c,p [−] 0.75

Wall friction coefficient

c,pw [−] 0.5

Coefficient of restitution ep [−] 0.9

Rolling resistance coefficient

r,p [−] 0.15

Fig. 4. Geometry of the conical hopper as used in [15], including its geo-

metric properties.

V. Brandt et al.

Powder Technology 436 (2024) 119447

simulations, again only the refinement step (see subsubsection 3.3.2)

was performed. The refinement zones were set according to the original

publication by placing the refinement zone from l=8 to l=4 at a height

of 0.055 m and the refinement zone from l=4 to l=2 at a height of

0.04 m (shown in Fig. 8). Another refinement zone placement was also

tested, which refined grains from l=8 to l=4 at a height of 0.05 m and

from l=4 to l=2 at a height of 0.035 m, referred as l=8/4/22. As

contact force scaling model, we used the l3- scaling to achieve compa-

rable results with De et al. [15]. Table 3 summarizes the simulation

conditions and the grain properties of the reference simulation for the

conical hopper.

The contact parameters for grain-grain and grain-wall collisions for

the simulation with CG using l=2 are given in Table 4. A DEM time step

of 2⋅10−6 s was used for all simulations (CG, MCG). According to De

et al. [15], no rolling resistance was applied. The overlap threshold for

the grain growth was set to 1.5% of the diameter of the scaled down

grains.

4.3. Wurster coater

The Wurster coater is a bottom spray granulator with an inserted

Wurster tube, which creates a circulating flow pattern of granular solid.

One of its applications is the coating of granular materials in the phar-

maceutical industry. The geometry of the Wurster coater considered

here was based on the work of Fries et al. [32] and is shown in Fig. 5.

For modeling the Wurster coater, the DEM was coupled to the

commercial CFD software Ansys Fluent 2022, where the fluid was

assumed to be incompressible. The fluid considered was air with a

density of

f=1.225 kg/m

and a viscosity of

f=1.7894⋅10−5 kg/m/s.

Two-way coupled DEM-CFD simulations were performed using the

particle-fluid interaction forces as the momentum source in a cell-based

way according to Kruggel-Emden and Oschmann [42]. The interaction

forces and the momentum source were calculated as described in sub-

sections 2.1 and 2.2. This requires knowledge about the particle and

fluid velocities, as well as the particle positions. User-defined-functions

provided by Ansys Fluent transfer the fluid velocities from the CFD to

the DEM. The mean weighted particle-fluid interaction forces and void

fractions are transferred to the CFD. For computational reasons, the void

fractions were calculated by assuming that the particles / grains are

volume equivalent porous cubes (see e.g. [43,44]). Obtained void frac-

tions were then mapped to the underlying CFD mesh. In the MCG sim-

ulations, scaled down particles / grains may be present due to the

coarsening and refinement of particles / grains as described in sub-

subsections 3.3.1 and 3.3.2. For the void fraction calculation, the actual

scaled down diameter ds was considered in this case.

A poly-hexcore CFD mesh of 67,500 cells was used for all simula-

tions. For the CG and MCG simulations of the Wurster coater, the CFD

cell size was not scaled because a lower mesh resolution led to conver-

gence issues due to the large inlet velocities. The cross section of the

used CFD mesh is shown in Fig. 5. The CFD time step was set to 2⋅10−4 s

using a fixed ratio of 1:100 with respect to the DEM time step.

Three velocity inlets are located at the bottom plate and one velocity

inlet is located at the nozzle tip (see Fig. 5). The nozzle is located inside

the Wurster tube. The inlet velocity values are shown in Fig. 5. To ensure

smooth fluidization of the granular solid and better convergence of the

CFD solver, the inlet velocity was ramped. A pressure outlet condition

was placed at the top of the system, while within the DEM a wall was

assumed to seal the outlet. All walls were represented as no slip

boundaries. Particles with a diameter of dp=0.002 m were considered.

As can be seen in Fig. 5, a gap of 0.01 m exists below the Wurster tube. As

granular solid must pass through this gap, the gap size strongly in-

fluences the recirculation of granular solid within the Wurster coater.

Larger gap sizes increase the amount of granular solid dragged through

the Wurster tube by the Venturi effect, because of the negative pressure

gradient towards the Wurster tube. Therefore, this system is predestined

for the application of MCG, because the gap size determines the local

maximum CG factor of the system.

For the Wurster coater, no comparison of the evaluation parameters

with the literature [32] was performed because the considered evalua-

tion parameters used in the original study are only determined for a

small scale system with 45,000 particles while our systems involves

150,000 particles. In our study, we therefore compared the results of the

simulation without CG as a reference with CG and MCG simulations for

several evaluation parameters. CG factors of l=2, l=4 and l=6 were

analyzed. MCG simulations with l=2/1, l=4/2/1, l=4/2 and l=6/3

were performed. The placement of the coarsening and refinement zones

is shown in Fig. 11e. For the Wurster coater, the CG contact force scaling

model SDS from Benyahia and Galvin [18] was chosen, because this

scaling model was identified by Brandt et al. [23] as the most accurate

one. Other contact force scaling models SDS-fric, l3, l2 (see section 3.1)

were also tested to evaluate their influence on the MCG results. The

hydrodynamic force scaling follows subsection 3.2 using the l3- scaling

of drag and pressure gradient force. The simulation conditions for the

Wurster coater and the particle properties of the reference simulation

are summarized in Table 5.

The contact parameters for particle-particle and particle-wall colli-

sions for simulations without CG are listed in Table 6, allowing a DEM

time step of 2⋅10−6 s as in the original work of Fries et al. [32]. All

simulations (without CG, CG and MCG) were performed with the afore

stated DEM time step. Due to the fluidization, the voidage is higher than

in the previous described hopper systems. Therefore, the overlap

threshold for the particle / grain growth was set to 0.1% of the diameter

of the scaled down particles / grains.

5. Validation of the new MCG approach for hopper discharge

In this section, we validate the derived MCG approach with the re-

sults of both previous MCG studies on discharge from a rectangular

hopper [34] (see subsection 5.1) and a conical hopper [15] (see sub-

section 5.2). For validation, we compared the discharge rates of simu-

lations without CG, CG and MCG hopper simulations.

5.1. Rectangular hopper

The particle / grain distribution during discharge from the rectan-

gular hopper at tsim =5 s is shown in Fig. 6. Looking at the simulations

without CG and the CG simulations with l=2 and l=4, the reduction of

the number of modeled entities with an increasing CG factor is visible.

For the MCG simulation with l=4/2/1, the refinement from l=4 to l=

2 can be observed in the specified refinement zones.

Table 3

Simulation conditions for discharge from the conical hopper according to [15].

Reference simulation (CG with l=2):

Number of grains 1,060,704

Grain density in kg/m

2500

Grain diameter in m 32 .10−6

Simulation time tsim in s 10

Contact force model Linear-spring-dashpot

CG and MCG simulations:

CG simulations l=4, l=8

MCG simulations l=8/4/21, l=8/4/22, l=16/8/4

Contact force scaling l3 [20]

Table 4

Contact parameters for grain-grain and grain-wall interactions for

modeling discharge from the conical hopper.

Normal spring stiffness kn

g in N/m 80

Tangential spring stiffness kt

g in N/m 80

Grain friction coefficient

c,g [−] 0.3

Coefficient of restitution eg [−] 0.9

V. Brandt et al.

Powder Technology 436 (2024) 119447

Obtained results of the discharge rates for the considered simulations

are shown in Fig. 7. In Fig. 7a, the discharge rate over a time of tsim =0 s

to 2 s is depicted to reveal details as in the original work [34]. To get

mean values of the discharge rate an averaging over 10 s is performed.

Doing so, the simulation without CG leads to a mean discharge rate of

m=0.17 kg/s with small oscillations (see Fig. 7a). When CG is applied,

the mean discharge rate decreases, due to the larger diameter of the

grains relative to the fixed size of the orifice. For a CG factor of l=2, a

mean discharge rate of ˙

m=0.126 kg/s is obtained. A further reduction

in the mean discharge rate occurs when a higher CG factor of l=4 is

used, resulting in a discharge rate of ˙

m=0.069 kg/s. For the CG sim-

ulations, higher fluctuations in discharge rates can be observed (see

Fig. 7a). For l=4 even temporary choking of the discharge is visible

around tsim =1.4 s, resulting in a mass flow of zero. The discharge rate of

the MCG simulation ( ˙

m=0.167 kg/s) however shows a good agreement

with the discharge rate and fluctuation of the simulation without CG

(see Fig. 7a).

In Fig. 7b the relative deviations δ of the CG and MCG simulations

and the results of Queteschiner et al. [34] to the simulation without CG

in form of δ=(x−xref )/xref are shown. Thereby xref =0.17 kg/s and x is

the respective mass flow rate for the CG or MCG simulation averaged

over 10 s. Because the discharge rate decreases for the CG simulations, a

negative relative deviation of δ= − 0.259 is obtained for l=2 and δ=

−0.594 for l=4 (see Fig. 7b). The discharge rate of the MCG simula-

tion results only in a small relative deviation of δ= − 0.018. In [34], a

MCG simulation with l=4/2/1 leads to a relative deviation of the

discharge rate in relation to their own simulation without CG (xref =

0.169 kg/s) of δ= − 0.006 (see Fig. 7b), which is in the same order of

magnitude as our MCG result.

It can be concluded that despite the differences in the MCG ap-

proaches, the results obtained were consistent with the work of Que-

teschiner et al. [34]. MCG was able to predict the discharge rate more

accurately than CG simulations, while still drastically reducing the

number of modeled entities.

5.2. Conical hopper

The grain distribution of the performed CG and MCG simulations of

the discharge from the conical hopper at tsim =5 s is shown in Fig. 8.

Note that the reference case considered for the conical hopper is a CG

simulation with l=2. With a rising CG factor from left to right, the

number of grains decreases. In the MCG simulation with l=8/4/21, the

refinement zones are set according to De et al. [15], as explained in

section 4. A further reduction in modeled grains can be achieved by

moving the refinement zones downward or by using higher CG factors.

This is tested in the MCG simulations with l=8/4/22 and l=16/8/4,

where the refinement zones are placed at heights of 0.05 m and 0.035 m.

At a time of tsim =5 s, the hopper fill level in the considered simulations

is unequal because the different discharge rates result in an inconsistent

cumulative discharge, as described in the following.

Our results of the obtained discharge rates over time of the different

CG and MCG simulations are shown in Fig. 9a for a time of tsim =0 s to 2

s as for the rectangular hopper. The mean values of the discharge rate

are then obtained by averaging over 10 s. The utilization of a CG factor

of l=2 leads to a stable discharge with a mean value of ˙

m=0.00366 kg/

s (see Fig. 9a). Similar to the results for the rectangular hopper, an in-

crease of the CG factor causes the discharge rate to decrease for CG

factors of l=4 and l=8, resulting in discharge rates of ˙

m=0.00288 kg/

s and ˙

m=0.00198 kg/s, respectively. In Fig. 9a, also higher fluctuations

Fig. 5. Geometry of the Wurster coater according to [32], including the values of the inlet velocities and the cross section of the CFD mesh.

Table 5

Simulation conditions for the Wurster coater.

Reference simulation (no CG):

Number of particles 150,000

Particle density in kg/m

1500

Particle diameter in m 0.002

Simulation time tsim in s 7

Contact force model Linear-spring-dashpot

CG and MCG simulations:

CG simulations l=2, l=4, l=6

MCG simulations l=2/1, l=4/2/1, l=4/2, l=6/3

Hydrodynamic force scaling l3 [37]

Contact force scaling SDS [18], SDS-fric [18,19], l3 [20], l2 [38]

Table 6

Contact parameters for particle-particle and particle-wall interactions

for simulations of the Wurster coater.

Normal spring stiffness kn

p in N/m 49,800

Tangential spring stiffness kt

p in N/m 15,920

Particle friction coefficient

c,p [−] 0.1

Coefficient of restitution ep [−] 0.8

Rolling resistance coefficient

r,p [−] 0.01

V. Brandt et al.

Powder Technology 436 (2024) 119447

in the discharge rates for increasing CG factors can be observed. For l=

8, a temporal choking of the discharge is observed at tsim =1.7 s (see

Fig. 9a), which is caused by the large grain size. When MCG is applied

with l=8/4/21, the refinement of the grains in two steps leads to a

similar discharge rate ( ˙

m=0.00364 kg/s) and fluctuations as for the CG

simulation of l=2 (see Fig. 9a). The downward displacement of the

refinement zones in the MCG with l=8/4/22 results in identical results,

confirming that this zone configuration is also reasonable. MCG with l=

16/8/4 leads to a similar discharge rate and resulting fluctuations as the

considered CG simulation with l=4, which can be seen in Fig. 9a.

Fig. 9b illustrates the relative deviations δ=(x−xref )/xref of the

mean discharge rate in relation to the reference simulation (xref =

0.00366 kg/s), which is the CG simulation with l=2 for the conical

hopper. For CG factors of l=4 and l=8, Fig. 9b shows rising relative

deviations of δ= − 0.213 and δ= − 0.459 for an increasing CG factor.

Both MCG simulations with l=8/4/21 and l=8/4/22 result in relative

deviations of only δ= − 0.005. A placement of the refinement zones too

close to the orifice can lead to an erroneous discharge rate and fluctu-

ations as described by De et al. [15]. In Fig. 9b, the CG simulation with

l=4 results in a similar relative deviation as the MCG with l=16/8/4

(δ= − 0.208). For comparison, Fig. 9b also illustrates the relative de-

viations of the MCG simulation of De et al. [15] in black related to their

results of the CG simulation with l=2 (xref =0.00322 kg/s). This leads

to a relative deviation of δ=0.003, which is in the same order of

magnitude as our MCG result.

In Fig. 10a, we exemplarily present the averaged grain stress and

averaged grain velocity magnitude over the height of the conical hopper

for the MCG simulation of l=8/4/21 and l=2. For the averaged stress

both simulations show a comparatively equal development, except from

Fig. 6. Distribution of particles / grains for discharge from a rectangular hopper after t

sim =5 s for simulations with l=1, l=2, l=4, l=4/2/1 also depicting the

locations of the refinement zones for the MCG.

Fig. 7. Results for the discharge from a rectangular hopper for simulations with l=1, l=2, l=4, l=4/2/1 a) over time b) relative deviations (x−x

ref )/xref in

relation to the simulation without CG including the results of [34] in black for comparison.

V. Brandt et al.

Powder Technology 436 (2024) 119447

the refinement zones at a height of 0.055 m and 0.04 m, where small

peaks of the grains stress can be overserved. A reason for this might be

the defined overlap threshold during the grain growth in the refinement

step, which allows a reorientation and ensures a fast restoration of mass

conservation. For the averaged velocity magnitude smaller peaks close

to the refinement zones are visible. The elevated stresses close to the

refinement zones do not significantly lead to differences in the velocity

magnitudes since in the dense bed surrounding grains are damping the

resulted kinetic energy.

The distribution of the grain stress in Fig. 10b shows a similar

behavior as discussed above highlighting stress peaks at the refinement

zones in the MCG simulation with l=8/4/21. Overall, both simulations

show minor differences in the stress distribution.

Since De et al. [15] provide the computation times of the simulations

performed, the computation times for the simulations considered in this

study are listed in Table 7 for comparison. Normalized values tn=tl=x/

tl=2 of the computation time are given with simulations performed on a

single core.

For the CG simulation of l=4, the normalized computation time

obtained by De et al. [15] is about 0.1. In our study, we obtain a

normalized value of 0.125, which is slightly higher but still of the same

order of magnitude. Applying CG with a CG factor of l=8 further re-

duces the normalized computation time in our study to 0.014. In the

study by De et al. [15], no reference computation time value is given,

because the hopper was already chocked at this CG factor. The

normalized computation time of the MCG simulation with l=8/4/2

conducted by De et al. [15] is about 0.08. In the present study, an even

lower value of 0.04 for the normalized computation time could be

achieved. The downward displacement of the refinement zones in MCG

with l=8/4/22 further reduces the number of modeled grains, resulting

in a normalized computation time of only 0.022. The simulation with l=

16/8/4 reduces the normalized computation time to only 0.004 and still

obtains reasonable results. Again, no benchmark is given in De et al.

[15], because comparable MCG simulations were not performed. The

obtained results show that the CG approach used in this study scales the

computation time in the same order of magnitude as in the study of De

Fig. 8. Distribution of grains for discharge from a conical hopper after t

sim =5 s for simulations with l=2, l=4, l=8, l=8/4/21, l=8/4/22 and l=16/8/4 also

containing the locations of the refinement zones for the MCG.

Fig. 9. Results for the discharge from a conical hopper for simulations with l=2, l=4, l=8, l=8/4/2

1, l=8/4/22 and l=16/8/4 a) over time b) relative

deviations (x−xref )/xref in relation to the CG simulation with l=2 including results of [15] in black for comparison.

V. Brandt et al.

Powder Technology 436 (2024) 119447

et al. [15] and provides an even better scaling if MCG is applied.

In conclusion, the MCG simulations with l=8/4/21 and l=8/4/22

correctly predicted the discharge rate of the conical hopper and its

fluctuation. The relative deviation of the discharge rate of the MCG

simulation of l=8/4/21 showed a good agreement with results of De

et al. [15]. The resulting level of accuracy is determined by the lowest

CG factor used in the discharge zone. Therefore, an MCG simulation of

l=16/8/4 cannot result in more accurate results than a CG simulation

with a global CG factor of l=4. Optimizing the position of the refine-

ment zones can save additional computational cost by further reducing

the number of entities to be modeled.

6. Application of the new MCG approach to a Wurster coater

The results of the Wurster coater system are analyzed below. First,

results on the spatial distribution of particles / grains are presented (see

subsection 6.1), followed by results on further evaluation parameters

(see subsection 6.2) and the comparison of computation times (see

subsection 6.3).

6.1. Analysis of the spatial distribution of particle / grains in the Wurster

coater

Plots of the spatial distribution of particles / grains in the Wurster

coater are shown in Fig. 11 at tsim =5 s for simulations without CG, CG

and MCG simulations. For a better visualization, particles / grains with

|y| ≥0.02 m are blanked, to depict the inflow region of the Wurster tube.

Fig. 11a shows the particle distribution for a simulation without CG. The

particle flow inside the Wurster tube is directed upward. The particles

collide with the wall at the top of the Wurster coater (see section 4) and

then move downward into the denser particle bed, which is forming

around the Wurster tube. When CG is applied, the ratio between the

Wurster tube gap size and the particle / grain size changes, which can be

seen for CG factors of l=2 and l=4 in Fig. 11b and Fig. 11c. When the

grain size exceeds the gap size, the particle circulation in the Wurster

coater is fully suppressed which can be seen in Fig. 11d for a CG factor of

l=6. The distribution of particles / grains in the Wurster coater for MCG

with l=2/1 is shown in Fig. 11e. In the refinement zone with dz=0.12

m and height hz=0.1 m, grains with l=2 are refined to particles of l=1

to resolve the critical inflow region of the Wurster tube with sufficient

high resolution. In the region outside of dz=0.12 m and of hz=0.1 m,

particles of l=1 are bundled together as part of the coarsening to form

grains with l=2, passing intermediate CG factors of l=1.26 and l=

1.58 thereby as shown in Fig. 1. Fig. 11f shows results of a MCG simu-

lation with three CG factors of l=4/2/1 applied, where the lowest

resolution of the granular solid with l=4 prevails at heights hz≥0.1 m.

Zones with l=1 and l=2 are identical to the MCG simulation with l=

2/1 (see Fig. 11e) for hz<0.1 m. In Fig. 11f upstream of hz=0.1 m,

some particles / grains with a lower CG factor than l=4 can be iden-

tified by their color and size, representing particles / grains that did not

find an association partner in the nearest neighbor search of their sur-

rounding DEM cells. This can also include grains with intermediate CG

factors in between l=1 and l=4. Particles / grains with l<4 are

further coarsened to grains with l=4, when other particle / grains with

the same CG factor enter their surrounding DEM cells. Particle / grain

distributions obtained for MCG with l=4/2 and l=6/3 are shown in

Fig. 11g and Fig. 11h, where the coarsening and refinement zones have

the same locations as in the MCG simulation with l=2/1. While the CG

simulation of l=6 out of obvious reasons fails to model the cyclic flow

pattern of granular solid through the Wurster coater, MCG simulations

with l=4/2 and even with l=6/3 are still able to lead to visually

reasonable results. This highlights the sense of using different grain sizes

Fig. 10. Comparison of the CG simulation with l=2 and the MCG simulation with l=8/4/2

1 of the conical hopper showing the a) averaged grain stress and

averaged grain velocity magnitude over height and b) distribution of the grain stress.

Table 7

Comparison of the normalized computation times tn=tl=x/tl=2 for the conical

hopper between the MCG approach of [15] and values obtained in the present

investigation using single core simulations.

De et al. [15] present investigation

l=4 0.095 0.125

l=8 – 0.014

l=8/4/2 0.083 0.040

l=8/4/22 – 0.022

l=16/8/4 – 0.004

V. Brandt et al.

Powder Technology 436 (2024) 119447

in a system with geometric constrictions. The local geometry to grain

size ratio has an influence on the obtained CG results [23]. At least when

the grain size approaches the size of the smallest cross-section of the

system, unphysical behavior compared to unscaled particles can occur

due to the artificially enlarged grains. This can be addressed by the local

refinement of grains in the MCG approach.

6.2. Analysis of the residence time in the spray zone and of the circulating

mass flow

For the comparison between the simulations without CG, CG and

MCG simulations, several evaluation parameters were considered. For

this reason, the residence time of particle / grains in a biconical spray

zone in between tsim =2 s to 7 s, as defined by Fries et al. [32] and the

circulating solids mass flow averaged between tsim =2 s to 7 s were

monitored.

The biconical spray zone is located above the nozzle inlet and rep-

resents the zone where particles would be coated. The residence time

distribution and the mean residence time of particles / grains in the

spray zone for simulations without CG, CG and MCG simulations are

compared. In coating processes, a narrow residence time distribution

resulting in a homogeneous coating layer is desired. The defined

biconical spray zone is shown in Fig. 12. The spray zone is located above

the nozzle inlet (see Fig. 5) with an opening angle of 40◦and a length of

0.045 m.

The results of the residence time distribution and the circulating

mass flow of the Wurster coater simulations are summarized in Fig. 13.

Since the CG simulation with l=6 failed to model the essential circu-

lating flow pattern of the Wurster coater, its results are not presented.

Fig. 13 shows the results of the residence time distribution of parti-

cles in the spray zone of the Wurster coater from tsim =2 s to 7 s for the

performed simulations. A typical residence time distribution curve of a

Wurster coater system can be observed for the simulation without CG

with l=1 (for comparison see [32]). Since only 5 s are evaluated, about

a third of the total number of particles do not pass the spray zone at all,

expressed by a residence time of zero. The first peak with a residence

time of 0.014 s represents the number of particles that pass the spray

zone once. Smaller peaks at multiples of this time period represent the

number of particles that pass the spray zone more than once. When

grains pass through the spray zone in CG and MCG simulations, the

residence time of all particles bundled within the grains is increased. For

the CG simulation of l=2, small deviations in the residence time dis-

tribution can be observed in comparison to the simulation without CG

(see Fig. 13). Applying CG with l=4 leads to more pronounced de-

viations of the residence time distribution as less grains pass the spray

zone once. The MCG approach with l=2/1 shows an almost identical

residence time distribution as the simulation without CG (see Fig. 13), as

both simulations use the same particle size in the critical inflow region of

the Wurster tube. In the MCG simulation with l=4/2/1, more

Fig. 11. Plots of the particle / grain distribution in the Wurster coater at t

sim =5 s for selected simulations with l=1, l=2, l=4, l=6, l=2/1, l=4/2/1, l=4/2

and l=6/3.

Fig. 12. Geometry of the biconical spray zone of the Wurster coater located

above the nozzle inlet according to [32].

V. Brandt et al.

Powder Technology 436 (2024) 119447

pronounced deviations are observed compared to the MCG simulation

with l=2/1, although only the grains moving upstream are further

coarsened to l=4 after passing the spray zone. This shows the vulner-

ability of the residence time distribution to small changes in predefined

CG factors in certain regions. For the MCG simulation with l=4/2, an

altered residence time distribution can be observed in Fig. 13. Overall,

fewer particles pass through the spray zone once, while the number of

particles with more than one cycle increases. This trend is further

developed when MCG is performed with l=6/3. While for grains with

l=6 alone no reasonable results are obtained, for l=6/3 a circulating

flow of granular solid is resulting. The simulation with l=6/3 provides

the lowest resolution of the granular solid in the MCG simulations and

therefore leads to results with the lowest level of accuracy.

In Fig. 14a, the relative deviations of the mean residence time related

to the simulation without CG in form of δ=(x−xref )/xref are shown.

Increasing relative deviations of the mean residence time of δ=0.072

and δ=0.267 with a rising CG factor are obtained for the CG simulations

with l=2 and l=4, respectively (see Fig. 14a). MCG with l=2/1 leads

to the overall most accurate mean residence time compared to the

simulation without CG with δ=0.054. When higher CG factors are

introduced in the MCG simulations, more pronounced relative de-

viations of the mean residence time can be observed for MCG with l=4/

2/1 and l=4/2 resulting in δ=0.09 and δ=0.126 in accordance with

Fig. 13. Since the finest resolution in MCG with l=4/2 is represented by

grains with l=2, this leads to higher relative deviations than l=4/2/1,

where particles with l=1 are present. MCG with l=6/3 leads to the

highest relative deviations of δ=0.342, confirming that for l=6/3 the

resolution is too low to obtain accurate results.

In Fig. 14b, the relative deviations of the mean circulating mass flow

related to the simulation without CG are presented. Since the residence

time in the spray zone and the circulating mass flow though the Wurster

tube are correlated, related results are expected. CG simulations with l=

2 and l=4 result in relative deviations of δ=0.08 and δ=0.212. MCG

with l=2/1 leads to the most accurate prediction of the mean circu-

lating mass flow with relative deviations of δ=0.03 compared to the

simulation without CG. Using higher CG factors in the MCG simulations,

the lower resolution of the granular solid leads to more pronounced

relative deviations for the MCG simulations with l=4/2/1 and l=4/2

of δ=0.072 and δ=0.106. Although for MCG with l=4/2 no particles

of l=1 are involved, the simulation still obtains reasonable results. For

the MCG simulation of l=6/3 the high relative deviations of δ=0.218

indicate that the resolution of the granular solid is too low to gain ac-

curate results.

In the following, the effect of the used CG contact force scaling model

on MCG simulations is analyzed. As a representative, the MCG simula-

tion with l=4/2 was chosen for this analysis as it provides reasonable

results and drastically reduces the number of modeled entities. As an

evaluation parameter, the relative deviations of the mean circulating

Fig. 13. Results of the residence time distribution in the spray zone obtained over 5 s for the simulations of the Wurster coater.

Fig. 14. Results of simulations for the Wurster coater for a) the relative deviation of the mean residence time obtained over 5 s related to the simulation without CG

and b) the relative deviation of the mean circulating mass flow averaged over 5 s with regard to the simulation without CG.

V. Brandt et al.

Powder Technology 436 (2024) 119447

mass flow for different CG contact force scaling models are considered,

which are shown in Fig. 15. The SDS-fric and l3- contact scaling models

tend to lead to more pronounced relative deviations of δ=0.132 and δ=

0.138, compared to the primarily used SDS- contact scaling model. The

friction scaling used in the SDS-fric scaling model [19] reduces the

friction coefficient with an increasing CG factor, resulting in a higher

circulating mass flow. The l2- contact scaling model shows slightly better

results (δ=0.099) than the primarily used SDS- contact scaling model.

Overall, a relatively small dependence on the results for the circulating

mass flow is observed for the different CG contact force scaling models

used.

In the following, the effect of the placement of coarsening and

refinement zones for MCG in the Wurster coater is analyzed. Therefore,

MCG simulations with l=4/2 are performed with different configura-

tions of coarsening and refinement zones analyzing the influence on the

relative deviations of the mean circulating mass flow related to the

simulation without CG (see Fig. 16 and Fig. 17).

In Fig. 16a, the spatial grain distribution after 5 s for the first

coarsening and refinement zone configuration is shown. Here, the

refinement zone is set at a height of hz=0.1 m, which is about the height

of the dense bed forming around the Wurster tube. In the refinement

zone, grains with l=4 are refined to grains with l=2. Outside the

refinement zone, grains with l=2 are coarsened to l=4. In Fig. 16b, the

refinement zone for grains with l=4 were set at a height of hz=0.05 m,

which is about the half of the dense bed height. Inside the Wurster tube

the grains are coarsened at a height of hz>0.1 m after the spray zone,

while around the Wurster tube grains are already coarsened at a height

of hz>0.05 m. The third configuration is shown in Fig. 16c. Here, the

initial configuration of coarsening and refinement zones as described in

section 6.1 were used. In Fig. 16d, the fourth configuration is depicted,

which is a combination of the second and third configuration, refining

grains at a height of hz=0.05 m and a diameter of dz=0.12 m. The fifth

configuration, where grains with l=4 are already refined at diameters

below dz=0.09 m and a height of hz=0.1 m is shown in Fig. 16e. The

sixth configuration is shown in Fig. 16f, where grains with l=4 are

refined directly before entering the Wurster tube at a height of hz=0.02

m and dz=0.09 m.

In Fig. 17, the relative deviations of the mean circulating mass flow

with regard to the simulation without CG are shown for the different

coarsening and refinement zone configurations depicted in Fig. 16. The

first configuration leads to the lowest relative deviation of δ=0.09 of all

tested configurations, because all grains in the dense bed are resolved

with l=2. Placing the coarsening and refinement zones as shown in

Fig. 15. a) Effect of CG contact force models (SDS, SDS-fric, l

3, l2) used in the

Wurster coater MCG simulations with l=4/2 on the relative deviations of the

mean circulating mass flow related to the simulation without CG averaged over

5 s.

Fig. 16. Distribution of grains in the MCG simulations with l=4/2 using different coarsening and refinement zone configurations (a-f) at t

sim =5 s.

V. Brandt et al.

Powder Technology 436 (2024) 119447

Fig. 16b and Fig. 16c, results in δ=0.108 and δ=0.106, respectively.

Both configurations reduce the number of modeled grains, while

maintaining the level of accuracy of the circulating mass flow. When the

resolution is further reduced and coarsening and refinement zones get

closer to the critical inflow region of the Wurster tube like in the fourth

configuration shown in Fig. 16d, a slightly elevated relative deviation of

δ=0.118 can be observed in Fig. 17. The coarsening and refinement

zone configuration shown in Fig. 16e leads to a more pronounced

relative deviation of δ=0.139. These results emphasize the sensitivity

of the Wurster coater MCG simulations to small changes in the place-

ment of the coarsening and refinement zones compared to the third

configuration. This trend is further developed for the sixth configuration

in Fig. 16e, where the number of modeled grains is drastically reduced

compared to the first configuration in Fig. 16a. The sixth configuration

results in the lowest level of accuracy with a relative deviation of δ=

0.166. When coarsening and refinement zones are placed close to critical

geometry constrictions the granular flow is disturbed as force chains

may be interrupted in the replacement process. Moreover, since parti-

cles / grains were placed initially shrunk they may prevail still shrunken

when facing the geometry constriction.

In can be concluded that the coarsening and refinement zones need

to be placed at a sufficient distance upstream from the critical geometry

constrictions. When coarsening and refinement zones are placed at a

distance <4…5 times the gap size of the Wurster tube (hg=0.01 m),

which is the case for the fifth and sixth configurations, the relative de-

viations are strongly rising. This agrees well with the findings of De et al.

[15], where the refinement zone was placed at a distance at least four

times the size of the orifice in the conical hopper.

To derivate general rules on the placement of coarsening and

refinement zones more detailed investigations on multiple granular

systems on different grain resolutions are required to elaborate on this.

This is something not aimed for in the scope of this study.

6.3. Comparison of computation times

In Fig. 18, the DEM-CFD and DEM computation times as well as the

summed up number of particles / grains normalized by the numbers of

particles of the simulations without CG are shown for the simulation

cases as depicted in Fig. 11, thereby relying on a coarsening and

refinement zone configuration for MCG as in configuration 3 (see Fig. 16

and Fig. 17). Normalized values of the computation time and of the

particles / grains are calculated in the form of tn=tl=x/tl=1 and nn=

nl=x/nl=1, respectively. Coupled single-core DEM and single-core CFD

simulations are performed to gain a representative comparison between

the computation times.

The normalized coupled DEM-CFD computation time of the l=2

simulation is about 0.55 in relation to the original particle simulation,

although the number of modeled particles / grains is reduced from

150,000 particles to 18,750 grains (nn=0.125, see Fig. 18). Since all

CFD simulations use the same mesh (see Fig. 5), the coupled DEM-CFD

computation time is limited to a normalized value about 0.55 even at

higher CG factors. When CG is applied, the computation time of the

DEM-CFD simulation is limited by the numerical overhead from the CFD

simulation. The DEM spends most of its time waiting for the CFD solu-

tion to proceed. The normalized DEM computation time is about 0.125,

which scales linearly with its normalized number of particles / grains.

For l=4, the normalized DEM computation time is further reduced

to 0.027 by reducing the number of particles / grains to 2344. For l=6,

only 694 grains need to be modeled, resulting in a normalized DEM

computation time of 0.008. This shows the high speedup potential of CG,

however no reasonable results were obtained for l=6.

In general, for MCG simulations, the number of modeled particles /

grains is variable depending on the location of the coarsening and

refinement zones. During the MCG simulations, the number of particles

/ grains is nearly constant as soon as the circulating flow pattern of the

granular solid is fully developed at tsim =2 s. When MCG is applied with

l=2/1, the mean number of particles / grains is reduced to 62,679,

resulting in a normalized DEM computation time of 0.57. Compared to

the l=2/1 simulation, the l=4/2/1 MCG simulation has a slightly

lower normalized DEM computation time, since the normalized number

of particles / grains is further reduced from 0.41 to 0.36. MCG with l=

4/2 immensely reduces the normalized DEM computation time to 0.09,

Fig. 17. Effect of different coarsening and refinement zone configurations

depicted in Fig. 16 used in the Wurster coater MCG simulations with l=4/2 on

the relative deviations of the mean circulating mass flow related to the simu-

lation without CG averaged over 5 s.

Fig. 18. Normalized values of the DEM-CFD and DEM computation time t

n=tl=x/tl=1 and the normalized number of particles / grains nn=nl=x/nl=1 for the Wurster

coater coupled single core DEM and single core CFD simulations.

V. Brandt et al.

Powder Technology 436 (2024) 119447

while the normalized number of particles / grains is reduced to 0.05 (see

Fig. 18). With a 20 times lower number of particles / grains, MCG with

l=4/2 still provides a relatively high level of accuracy with relative

deviations of the mean residence time and circulating mass flow about

δ≈0.1 (see Fig. 14). For MCG with l=6/3, the normalized number of

particles / grains is further reduced to 0.017, which corresponds to an

average number of 2498 particles / grains, resulting in a lower

normalized DEM computation time of 0.045.

When performing multi-core DEM-CFD simulations, the reduction of

computational cost of the DEM caused by the application of CG and MCG

can be used to redistribute cores from the DEM to the CFD. This pro-

cedure results in a faster CFD computation, leading to a rebalancing of

computational cost and an overall faster DEM-CFD computation.

7. Conclusions

In this study, MCG was applied to two hoppers and a Wurster coater

system and benchmarked against simulations without CG and CG sim-

ulations on key evaluation parameters. The newly developed MCG

approach was validated in two hopper systems, where CG led to an

erroneous prediction of the discharge rate, due to a shift in the ratio

between the orifice size and the grain size. The use of MCG allowed the

grain size to be adapted to the local geometry constrictions. In the

rectangular hopper grains with a CG factor of l=4 were used far from

the orifice to reduce the number of modeled entities, which were then

refined to grains with l=2 and particles with l=1 before entering the

discharge zone to produce reasonable results. The MCG simulation led to

a similar discharge rate as the simulation without CG. In the conical

hopper MCG was applied to refine grains with l=8 to grains with l=4,

which were further refined to grains with l=2 close to the discharge

zone. This resulted in an up to 45 times faster computation, while

maintaining the same level of accuracy in the discharge rate and its

fluctuations as the fine-scale simulation with l=2. The MCG approach

presented in this study led to a similar or even greater reduction in

computational cost as the approach of De et al. [15], while providing

results with the same accuracy as both previous MCG studies [15,34].

MCG was applied to a complex Wurster coater system [32] simulated

by DEM-CFD for the first time, where both coarsening and refinement

steps were required to be used. For CG simulations with l=4, the nar-

row gap size between the base plate and the Wurster tube led to a wrong

prediction of the residence time in the spray zone and the circulating

mass flow. For CG simulations with l=6, the circulating flow pattern

was suppressed, because the grain diameter exceeded the gap size of the

Wurster tube. In MCG simulations, the interplay between the coarsening

and refinement zones provided a higher resolved granular flow through

critical regions of the Wurster coater. MCG with l=2/1 was able to

predict the evaluation parameters with the highest level of accuracy

compared to all considered CG and MCG simulations. At the same time,

for l=2/1, the modeled number of particles / grains was the highest of

the CG and MCG simulations leading to an only 1.5 times faster DEM-

CFD computation. MCG with l=4/2 was still able to accurately pre-

dict the residence time and circulating mass flow, while drastically

reducing the number of modeled entities, resulting in a 10 times faster

DEM computation time with respect to the simulation without CG.

Although the CG simulation with l=6 failed to model the circulating

flow pattern, MCG with l=6/3 was still able to produce reasonable

results.

The accuracy of MCG simulations is among other things dependent

on the contact force scaling model and coarsening and refinement zones

used. The placement of the coarsening and refinement zones must be

done carefully as it will inherently affect the obtained results and

computation times. For the Wurster coater system, we find that the

coarsening and refinement zones should be placed at a distance >4…5

times the gap size of the Wurster tube (hg=0.01 m) to obtain accurate

results.

Our study demonstrated that the novel MCG approach is capable of

modeling complex systems with geometric restrictions with a high level

of accuracy, while drastically reducing the DEM computation time. The

wide range of applications and the ability to adapt the resolution of the

granular solid depending on the geometry constrictions makes this

approach very promising for future use in industrial-scale particle sys-

tems. Establishing general rules for the placement of the coarsening and

refinement zones to balance the achieved accuracy and computation

time will be a topic of future work on MCG. This would allow the MCG

algorithm to automatically detect the coarsening and refinement zones,

eliminating the need to manually define them.

CRediT authorship contribution statement

V. Brandt: Conceptualization, Formal analysis, Investigation,

Methodology, Software, Validation, Visualization, Writing – original

draft, Writing – review & editing. J. Grabowski: Writing – review &

editing. N. Jurtz: Conceptualization, Funding acquisition, Supervision,

Writing – review & editing. M. Kraume: Conceptualization, Funding

acquisition, Writing – review & editing. H. Kruggel-Emden: Concep-

tualization, Funding acquisition, Methodology, Resources, Supervision,

Writing – original draft, Writing – review & editing.

Declaration of competing interest

The authors declare that they have no known competing financial

interests or personal relationships that could have appeared to influence

the work reported in this paper.

Data availability

Data will be made available on request.

Acknowledgments

Financial support was provided by the Deutsche For-

schungsgemeinschaft (DFG, German Research Foundation) – Project-ID

456827728. Computing resources were funded by the DFG – Project-ID

463921749.

References

[1] J. Grace, X. Bi, N. Ellis, Essentials of Fluidization Technology, Wiley, 2020, https://

doi.org/10.1002/9783527699483.

[2] P.A. Cundall, O.D.L. Strack, A discrete numerical model for granular assemblies,

G´

eotechnique. 29 (1979) 47–65, https://doi.org/10.1680/geot.1979.29.1.47.

[3] A. Klimanek, W. Adamczyk, S. Kallio, P. Kozołub, G. Węcel, A. Szlęk, Experimental

and numerical study of pseudo-2D circulating fluidized beds, Particuology. 29

(2016) 48–59, https://doi.org/10.1016/j.partic.2015.09.009.

[4] T. Oschmann, J. Hold, H. Kruggel-Emden, Numerical investigation of mixing and

orientation of non-spherical particles in a model type fluidized bed, Powder

Technol. 258 (2014) 304–323, https://doi.org/10.1016/j.powtec.2014.03.046.

[5] H. Kruggel-Emden, K. Vollmari, Flow-regime transitions in fluidized beds of non-

spherical particles, Particuology. 29 (2016) 1–15, https://doi.org/10.1016/j.

partic.2016.01.004.

[6] H. Kruggel-Emden, B. Kravets, M.K. Suryanarayana, R. Jasevicius, Direct numerical

simulation of coupled fluid flow and heat transfer for single particles and particle

packings by a LBM-approach, Powder Technol. 294 (2016) 236–251, https://doi.

org/10.1016/j.powtec.2016.02.038.

[7] B. Kravets, H. Kruggel-Emden, Investigation of local heat transfer in random

particle packings by a fully resolved LBM-approach, Powder Technol. 318 (2017)

293–305, https://doi.org/10.1016/j.powtec.2017.05.039.

[8] J. Lu, M.D. Tan, E.A.J.F. Peters, J.A.M. Kuipers, Direct numerical simulation of

reactive fluid-particle systems using an immersed boundary method, Ind. Eng.

Chem. Res. 57 (2018) 15565–15578, https://doi.org/10.1021/acs.iecr.8b03158.

[9] N. Jurtz, M. Kraume, G.D. Wehinger, Advances in fixed-bed reactor modeling using

particle-resolved computational fluid dynamics (CFD), Rev. Chem. Eng. 35 (2019)

139–190, https://doi.org/10.1515/revce-2017-0059.

[10] A. Esteghamatian, M. Bernard, M. Lance, A. Hammouti, A. Wachs, Micro/meso

simulation of a fluidized bed in a homogeneous bubbling regime, Int. J. Multiphase

Flow 92 (2017) 93–111, https://doi.org/10.1016/j.ijmultiphaseflow.2017.03.002.

[11] D. Gidaspow, Multiphase Flow and Fluidization: Continuum and Kinetic Theory

Descriptions, Acad. Press, 1994.

V. Brandt et al.

Powder Technology 436 (2024) 119447

[12] S. Ergun, Fluid flow through packed columns, Chem. Eng. Prog. 48 (2) (1952)

89–94.

[13] Y.H. Wen, C.Y. and Yu, Mechanics of fluidization, Chem. Eng. Prog. Symp. Ser. 62

(1) (1966) 100–111.

[14] T. Tsuji, K. Yabumoto, T. Tanaka, Spontaneous structures in three-dimensional

bubbling gas-fluidized bed by parallel DEM-CFD coupling simulation, Powder

Technol. 184 (2008) 132–140, https://doi.org/10.1016/j.powtec.2007.11.042.

[15] T. De, J. Chakraborty, J. Kumar, A. Tripathi, M. Sen, W. Ketterhagen, A particle

location based multi-level coarse-graining technique for discrete element method

(DEM) simulation, Powder Technol. 398 (2022), https://doi.org/10.1016/j.

powtec.2021.117058.

[16] L. Lu, J. Xu, W. Ge, Y. Yue, X. Liu, J. Li, EMMS-based discrete particle method

(EMMS-DPM) for simulation of gas-solid flows, Chem. Eng. Sci. 120 (2014) 67–87,

https://doi.org/10.1016/j.ces.2014.08.004.

[17] C. Bierwisch, T. Kraft, H. Riedel, M. Moseler, Die filling optimization using three-

dimensional discrete element modeling, Powder Technol. 196 (2009) 169–179,

https://doi.org/10.1016/j.powtec.2009.07.018.

[18] S. Benyahia, J.E. Galvin, Estimation of numerical errors related to some basic

assumptions in discrete particle methods, Ind. Eng. Chem. Res. 49 (2010)

10588–10605, https://doi.org/10.1021/ie100662z.

[19] R. Cai, Y. Zhao, An experimentally validated coarse-grain DEM study of

monodisperse granular mixing, Powder Technol. 361 (2020) 99–111, https://doi.

org/10.1016/j.powtec.2019.10.023.

[20] M. Sakai, S. Koshizuka, Large-scale discrete element modeling in pneumatic

conveying, Chem. Eng. Sci. 64 (2009) 533–539, https://doi.org/10.1016/j.

ces.2008.10.003.

[21] S. Radl, C. Radeke, J.G. Khinast, S. Sundaresan, Parcel-based approach for the

simulation of gas-particle flows, in: 8th Interantional Conf. CFD Oil Gas, Metall.

Process Ind, 2011, pp. 1–10.

[22] A. Di Renzo, E.S. Napolitano, F.P. Di Maio, Coarse-grain dem modelling in fluidized

bed simulation: a review, Processes. 9 (2021) 1–30, https://doi.org/10.3390/

pr9020279.

[23] V. Brandt, J. Grabowski, N. Jurtz, M. Kraume, H. Kruggel-Emden, A benchmarking

study of different DEM coarse graining strategies, Powder Technol. (2023) 118629,

https://doi.org/10.1016/j.powtec.2023.118629.

[24] M.J.A. de Munck, J.B. van Gelder, E.A.J.F. Peters, J.A.M. Kuipers, A detailed gas-

solid fluidized bed comparison study on CFD-DEM coarse-graining techniques,

Chem. Eng. Sci. 269 (2023) 118441, https://doi.org/10.1016/j.ces.2022.118441.

[25] N. Jurtz, H. Kruggel-Emden, O. Baran, R. Aglave, R. Cocco, M. Kraume, Impact of

contact scaling and drag calculation on the accuracy of coarse-grained discrete

element method, Chem. Eng. Technol. 43 (2020) 1959–1970, https://doi.org/

10.1002/ceat.202000055.

[26] M. Sakai, Y. Yamada, Y. Shigeto, K. Shibata, V.M. Kawasaki, S. Koshizuka, Large-

scale discrete element modeling in a fluidized bed, Int. J. Numer. Methods Fluids

64 (2010) 1319–1335, https://doi.org/10.1002/fld.2364.

[27] J.E. Hilton, P.W. Cleary, Comparison of non-cohesive resolved and coarse grain

DEM models for gas flow through particle beds, Appl. Math. Model. 38 (2014)

4197–4214, https://doi.org/10.1016/j.apm.2014.02.013.

[28] E.S. Napolitano, A. Di Renzo, F.P. Di Maio, Coarse-grain DEM-CFD modelling of

dense particle flow in Gas–Solid cyclone, Sep. Purif. Technol. 287 (2022) 120591,

https://doi.org/10.1016/j.seppur.2022.120591.

[29] Y. Kosaku, Y. Tsunazawa, C. Tokoro, Investigating the upper limit for applying the

coarse grain model in a discrete element method examining mixing processes in a

rolling drum, Adv. Powder Technol. 32 (2021) 3980–3989, https://doi.org/

10.1016/j.apt.2021.08.039.

[30] H.P. Zhu, Z.Y. Zhou, R.Y. Yang, A.B. Yu, Discrete particle simulation of particulate

systems: a review of major applications and findings, Chem. Eng. Sci. 63 (2008)

5728–5770, https://doi.org/10.1016/j.ces.2008.08.006.

[31] L. Fries, S. Antonyuk, S. Heinrich, D. Dopfer, S. Palzer, Collision dynamics in

fluidised bed granulators: a DEM-CFD study, Chem. Eng. Sci. 86 (2013) 108–123,

https://doi.org/10.1016/j.ces.2012.06.026.

[32] L. Fries, S. Antonyuk, S. Heinrich, S. Palzer, DEM-CFD modeling of a fluidized bed

spray granulator, Chem. Eng. Sci. 66 (2011) 2340–2355, https://doi.org/10.1016/

j.ces.2011.02.038.

[33] F. Alobaid, J. Str¨

ohle, B. Epple, Extended CFD/DEM model for the simulation of

circulating fluidized bed, Adv. Powder Technol. 24 (2013) 403–415, https://doi.

org/10.1016/j.apt.2012.09.003.

[34] D. Queteschiner, T. Lichtenegger, S. Pirker, S. Schneiderbauer, Multi-Level Coarse-

Grain Model of the DEM, 2018, https://doi.org/10.1016/j.powtec.2018.07.033.

[35] J.M. Link, W. Godlieb, N.G. Deen, J.A.M. Kuipers, Discrete element study of

granulation in a spout-fluidized bed, Chem. Eng. Sci. 62 (2007) 195–207, https://

doi.org/10.1016/j.ces.2006.08.018.

[36] Y.C. Zhou, B.D. Wright, R.Y. Yang, B.H. Xu, A.B. Yu, Rolling friction in the dynamic

simulation of sandpile formation, Phys. A Stat. Mech. Appl. 269 (1999) 536–553,

https://doi.org/10.1016/S0378-4371(99)00183-1.

[37] M. Sakai, M. Abe, Y. Shigeto, S. Mizutani, H. Takahashi, A. Vir´

e, J.R. Percival,

J. Xiang, C.C. Pain, Verification and validation of a coarse grain model of the DEM

in a bubbling fluidized bed, Chem. Eng. J. 244 (2014) 33–43, https://doi.org/

10.1016/j.cej.2014.01.029.

[38] K. Vollmari, H. Kruggel-Emden, Numerical and experimental analysis of particle

residence times in a continuously operated dual-chamber fluidized bed, Powder

Technol. 338 (2018) 625–637, https://doi.org/10.1016/j.powtec.2018.07.061.

[39] C. Bierwisch, T. Kraft, H. Riedel, M. Moseler, Three-dimensional discrete element

models for the granular statics and dynamics of powders in cavity filling, J. Mech.

Phys. Solids 57 (2009) 10–31, https://doi.org/10.1016/j.jmps.2008.10.006.

[40] S. Bußmann, H. Kruggel-Emden, M. Reichert, Realizing fragment spawning and

fragment growth in the parallelized discrete element method (DEM) during

modelling of comminution, Adv. Powder Technol. 32 (2021) 2171–2191, https://

doi.org/10.1016/j.apt.2021.04.033.

[41] S. Schneiderbauer, S. Puttinger, S. Pirker, Comparative analysis of subgrid drag

modifications for dense gas-particle flows in bubbling fluidized beds, AICHE J. 59

(2013) 4077–4099, https://doi.org/10.1002/aic.14155.

[42] H. Kruggel-Emden, T. Oschmann, Numerical study of rope formation and

dispersion of non-spherical particles during pneumatic conveying in a pipe bend,

Powder Technol. 268 (2014) 219–236, https://doi.org/10.1016/j.

powtec.2014.08.033.

[43] J.M. Link, L.A. Cuypers, N.G. Deen, J.A.M. Kuipers, Flow regimes in a spout-fluid

bed: a combined experimental and simulation study, Chem. Eng. Sci. 60 (2005)

3425–3442, https://doi.org/10.1016/j.ces.2005.01.027.

[44] N.G. Deen, M. Van Sint Annaland, M.A. Van der Hoef, J.A.M. Kuipers, Review of

discrete particle modeling of fluidized beds, Chem. Eng. Sci. 62 (2007) 28–44,

https://doi.org/10.1016/j.ces.2006.08.014.

V. Brandt et al.