- Email: [email protected]
A multiscale DEM-VOF method for the simulation of three-phase flows
Accepted Manuscript
A Multiscale DEM-VOF method for the simulation of three-phase flows Gabriele Pozzetti, Bernhard Peters PII: DOI: Reference:
S0301-9322(16)30732-7 10.1016/j.ijmultiphaseflow.2017.10.008 IJMF 2664
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
11 December 2016 9 October 2017 9 October 2017
Please cite this article as: Gabriele Pozzetti, Bernhard Peters, A Multiscale DEM-VOF method for the simulation of three-phase flows, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.10.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
A Multiscale DEM-VOF method for the simulation of three-phase flows Gabriele Pozzettia,∗, Bernhard Petersa a
Campus Kirchberg, Universit du Luxembourg 6, rue Richard Coudenhove-Kalergi L-1359 Luxembourg
CR IP T
Abstract
A novel multiscale approach for three-phase flows is presented. The goal of the proposed method is to solve arbitrary three-phase flow configurations in a computationally efficient way, and in particular taking into account the effects of different length scales while sharply reducing the computational burden. This is particularly important in chemical, environmental, and process engineering, where large-scale quantities are
AN US
normally of interest, but small-scale dynamics cannot be neglected. The method is based on the definition of two different length scales: the bulk scale, and the fluid fine scale. A dual-grid approach is adopted in order to resolve the bulk scale with information from the fluid fine scale. It is shown that the described method succeeds in delivering more accuracy than a standard approach based on the volume averaging technique, still, it remains suitable for the solution of real interest problems. The method is shown to successfully satisfy
M
experimental results presented in the literature. Examples of three-phase flows simulations are provided to show how the proposed numerical approach can describe the physics of particle-laden, free surface flows with
ED
competitive computational cost. It is shown how the proposed method can naturally extend the DEM-VOF method to the presence of complex interface dynamics.
PT
Keywords: Multiphase flows, CFD-DEM coupling, Grid convergence, Multiscale simulation, DEM-VOF
CE
method
Highlights
AC
• A multiscale DEM-VOF coupling is proposed to simulate three-phase flows • It produces grid-convergent solutions for the fluid phase • It can reproduce droplet formation and interface perturbations from a single particle • It is validated against three-phase flows presented in the literature
∗ Principal
corresponding author. Email: [email protected] Email address: [email protected] (Gabriele Pozzetti)
Preprint submitted to Elsevier
October 12, 2017
ACCEPTED MANUSCRIPT
1
1. Introduction Gas-liquid-particles flows are extremely common in industrial and environmental applications, and repre-
3
sent a standard configuration for chemical and process engineering. The engineering community is currently
4
demanding numerical tools able to predict the behavior of such complex flows at reasonable computational
5
cost. In particular, a direct numerical approach to all the phenomena that arise in these flows is unbearable
6
for nowadays computational capacities. Therefore, a considerable modeling effort is required to provide
7
reliable solutions with competitive cost.
CR IP T
2
The dispersed solid phase often plays a key role in the physics of the problem, being responsible for
9
the main chemical and dynamical interactions of the system. A method capable of describing the discrete
10
dynamics of the granular phase could, therefore, be extremely interesting [28]. For this reason, CFD-DEM
11
approaches have become very popular as they are able to combine the accuracy of the Discrete Element
12
Method with the highly developed CFD techniques for the study of complex flows.
AN US
8
Recently, significant efforts have been made for coupling DEM solvers with multiphase fluid solvers using
14
the Volume of Fluid (VOF) technique. This approach has been shown promising to provide reliable numerical
15
solutions of the three-phase flow dynamic and its influence on relevant chemical processes [20, 26, 11, 10, 17].
16
The VOF method is a relatively old approach for solving free surface flows, that after several modifications,
17
is still used today [22]. This gave rise to the development of multiple VOF schemes that try to match
18
a compromise between computational efficiency, and accuracy [22]. Generally, a VOF scheme aims to
19
reconstruct the position of the interface location. The interface reconstruction is particularly critical as it
20
is used to calculate the local curvature which is proportional to the surface tension force that will appear
21
as a contribution in the momentum equation. Furthermore, within a VOF scheme, a volume fraction field
22
is transported by the velocity field in order to define the new interface location.
CE
PT
ED
M
13
Two approaches, commonly referred to as algebraic and geometrical schemes, are widely used today.
24
Geometrical VOF techniques reconstruct the interface using volume fraction at each cell and its neighbors,
25
so to approximate interface information such as area and intersections with cell faces. Many variations of
26
27
AC
23
this strategy have been presented in the literature, and an extended description of various VOF schemes can be found in [22]. Algebraic advection methods perform a discretization of the advection equation by
28
reconstructing volume fluxes from the solution of the velocity field. This is used to compute the mixture
29
volume fraction and its flux over cells. A commonly used algebraic advection scheme is implemented in
30
R [25]. Various modification to the method have been proposed, the algebraic scheme adopted OpenFOAM
31
in this contribution is more extensively described in section 4. 2
ACCEPTED MANUSCRIPT
Depending on the technique chosen to couple the Eulerian and the Lagrangian domains, CFD-DEM
33
couplings can lead to extremely high computational cost [3, 23, 15]. In order to address industrially significant
34
problems, a widely accepted solution is the volume averaging approach[3]. It consists of the projection
35
of Lagrangian fields into Eulerian ones through the usage of volume-averaged variables[3]. CFD-DEM
36
approaches based on the volume averaging technique aim to perform a zero-dimensional coupling between
37
CFD and DEM solution, using cell values or locally interpolated values for the evaluation of the fluid-particle
38
interactions. This significantly reduces the computational burden and represents one of the main advantages
39
of the method. The VOF method, and in general interface-capturing methods, are particularly interesting
40
to be used within these couplings, as they allow solving the well-known single-phase equations within each
41
single fluid phase, smearing the interface region between two phases over a few cells. This can represent a key
42
advantage for a volume-averaged CFD-DEM coupling, as it allows adopting the highly developed empirical
43
models for the coupling between particles and single-phase fluid within the major part of the computational
44
domain.
AN US
CR IP T
32
Due to their attempt to perform inexpensive couplings, volume averaging approaches naturally operate
46
on a length scale at which the exchange of momentum between particles and fluid can be considered as
47
zero-dimensional. We will refer to this scale as the bulk or coupling scale of the simulation. We will instead
48
refer to fluid fine scale as the smallest length scale that needs to be resolved to produce grid independent
49
results for the fluid.
ED
M
45
In standard CFD-DEM couplings based on the volume averaging technique, the bulk scale is automat-
51
ically chosen as the one at which the fluid equations will be written and solved. Therefore, this technique
52
naturally applies to cases in which the fluid fine scale is significantly larger than the particles characteristic
53
dimensions. For this reason, the volume averaging approach was originally proposed for couplings in which
54
the CFD cells can be kept significantly larger than the largest particle diameter. However, it nowadays
55
became a standard to apply this coupling technique with cell-particle diameter ratios close to the unity, by
56
using complex averaging procedures [18].
58
59
CE
AC
57
PT
50
In a VOF based simulation, the fluid fine scale is associated with the interface dynamic, that must be
correctly captured to achieve grid convergence. Generally, the ability of a VOF scheme to reconstruct the interface is bounded to the grid spacing. The particular requirement will depend on the specific scheme
60
adopted, but in all the cases a correct reconstruction of the interface is mandatory as it will affect the cal-
61
culation of the surface tension terms, and therefore the determination of the velocity fields, and finally the
62
advection of the volume fraction. In addition to that, in algebraic VOF schemes, the usage of a sufficiently 3
ACCEPTED MANUSCRIPT
fine discretization is mandatory in order to deal with the problem of spurious currents. This can lead to an
64
erroneous calculation of the fluxes, and therefore a poor conservation of the overall scheme [2, 4]. Further-
65
more, correct grid resolution is particularly important in case of high density ratios. In these configurations,
66
small errors in the volume fraction can easily cause major problems for momentum conservation, leading in
67
the end to a poor accuracy of the solution [4]. For this reason, to extend DEM-VOF couplings to complex
68
interface scenarios, it is important to develop algorithms that can ensure the correct grid resolution in order
69
for the VOF scheme to properly converge.
CR IP T
63
In purely CFD cases, a higher grid resolution is directly linked to an increase in computational cost, that
71
represents a major challenge for the numerical solution of these configurations. In a DEM-VOF coupling
72
based on the volume averaging technique, the scale separation between the simulation domain (or integral
73
scale), the bulk scale for CFD-DEM coupling, and the fluid fine scale, can represent an even bigger problem.
74
This is due to the fact that the bulk scale is the smallest one that a single-scale volume averaging coupling
75
can resolve, but the solution of the fluid fine scale is needed to provide unique solutions.
76
77
AN US
70
The main objective of this contribution is to propose a multiscale method that overcomes this problem. The method we propose uses a dual-grid approach where:
• A coarse grid is used to perform the CFD-DEM coupling at the bulk scale
79
• A fine grid is used to resolve the equations governing the fluid phases at the fluid fine scale.
ED
M
78
Those two scales can have significantly different characteristic lengths, and a correct modeling for both of
81
them is mandatory in order to reach a compromise between accuracy and computational feasibility. This
82
dual-grid multiscale approach aims to capture the most significant phenomena of the two different scales,
83
while limiting as much as possible the computational burden.
PT
80
This method is meant to bring the most advantages when the separation between those two scales is more
85
marked. That happens when the smallest structures of the fluid flow that must be resolved to obtain grid
86
independent results, are significantly smaller than the ones determining the bulk fluid-particle momentum
88
AC
87
CE
84
exchange. In section 2, we refer to the Stokes number as a parameter to identify this scale separability. Compared to the existing methods, our novel multiscale approach is thought to have several advantages.
89
First, it allows achieving grid-convergence of macroscopic flow quantities in cases in which the single-scale
90
DEM-VOF method cannot. Second, it allows more accurate predictions of the interface dynamics that are
91
fundamental to capture the most significant fluid dynamic instabilities. This leads to a highly versatile
92
numerical method that can be applied to several flow configurations. The method can also be naturally 4
ACCEPTED MANUSCRIPT
extended to other CFD-DEM couplings once the length scales are properly individuated and defined. This
94
paper is organized as follows: First, the novel dual-grid coupling procedure is proposed in its more general
95
configuration. Then, the governing equations for both the fluid and the granular phases are presented.
96
The model for the coupling terms between continuum and discrete phases is subsequently proposed with
97
particular reference to the phase interchange problem, and to the interpolation strategies adopted. Some
98
reference cases are finally discussed. To the best of our knowledge the multiscale DEM-VOF method is the
99
first that allows capturing the effect of a single particle on a liquid interface, still keeping the computational
100
cost of a volume-averaged simulation.
101
2. Method description
CR IP T
93
A fundamental nondimensional number in the study of particle-laden flow is the Stokes number. It
103
represents the ratio between the particles and the fluids characteristic time scales and, when the particle
104
Reynolds number is less than unity, can be expressed as: St =
AN US
102
ρp d2p |uf | , dI 18µf
(1)
with ρp the density of the particles, dp their diameter, uf the fluid velocity, dI the hydraulic diameter,
106
µf the viscosity of the fluid. In this paper, we refer to the case of high Stokes number flows, for which a
107
single particle has relatively large mass and diameter and therefore its single action plays an important role
108
in the process dynamic. In this kind of flow, a particle’s motion is mainly influenced by large structures
109
of the fluid due to its inertia. Therefore, the coupling between the CFD and DEM can be operated on a
110
relatively coarse scale, where only the bulk structures of the flow are resolved. Nevertheless, when the fluid
111
phases are characterized by high-density ratios and complex interface dynamics, their bulk behavior cannot
112
be univocally determined without resolving small-scale dynamics. In particular, to correctly predict the
113
phase distribution, the interface must be correctly advected, and this requires detailed information on the
114
velocity field. Furthermore, the very solution of the velocity field is influenced by the local behavior of the
116
117
ED
PT
CE
AC
115
M
105
interface position and curvature. For these reasons, in order to obtain grid independent solutions, one often has to resolve a fluid scale smaller than the one used for the CFD-DEM coupling. In this contribution, we propose a dual-grid algorithm for the CFD-DEM coupling that overcomes this issue, and allows obtaining
118
grid convergent solutions for the fluid while keeping the computational advantages of a volume-averaged
119
simulation.
120
The dual-grid algorithm proposed to achieve this is described in figure 2, and structured as follows. A 5
ACCEPTED MANUSCRIPT
Table 1: List of variables
Symbol Variable
Symbol Variable
uf
fluid velocity
α
µf ρf dI Ffpi
uc Γ ψ C
TΓ
fluid viscosity fluid density hydraulic diameter volumetric source of momentum due to interaction with particles surface tension force
Fb p
body force pressure
σf τf
npc
number of particles in a cell
g
CR IP T
h
Dirac’s delta function describing the interface fluid stress tensor deviatoric part of the stress tensor gravitational acceleration modulus hight
kf
free surface compression factor
Φm Φjm
mass flux partial mass flux (relative to a sub-cycling step) uf F · AF
δTtot
face Area face Area normal vector (pointing outward the owner cell) compression velocity face value = uc · AF simulation time step
δTi
sub-cycle time step
aN
H
non diagonal contribution of the velocity matrix for momentum predictor
ρp dp
particle density particle diameter
AN US
ucF
η
M
AF AF
aP
ED
χ
P porosity field: i Vpi /Vcell with i = 1....npc phase indicator
PT
δΓ
volume fraction of the liquid phase compression velocity Surface Tension particle volume concentration Curvature
m
particle mass
Mext
I Fc Fij
particle moment of inertia contact force acting on a particle contact force acting on the particle i due to the collision with particle j force rising from the particlefluid interaction particle-fluid interface normal vector
up Ap φ
gravitational force torque acting on a particle due to collision events external torque acting on a particle particle velocity particle surface area particle orientation
ω
particle angular velocity
particle based Reynolds number interface local velocity interface normal vector
ρh β CD
density of the heavy phase drag factor drag coefficient
AC
CE
Fg Mc
the matrix coefficient related to the centroid (velocity equation) the coefficients of neighbor centroids (velocity equation)
Ff n
Re uint nint
6
ACCEPTED MANUSCRIPT
Bulk Scale
Fluid 2
CR IP T
Fluid 1
Fluid Fine Scale
AN US
Fluid 2
Fluid 1
M
Figure 1: Different length scales in high-Stokes three-phase flows: bulk (coarse) scale and fluid fine scale.
first coarse mesh is used to construct averaged fields from the granular-phase distribution, according to the
122
standard volume averaging technique. This mesh is used to solve the fluid-particle coupling consisting in
123
the estimation of the drag force, the calculation of the porosity field, and the definition of the particle-
124
related source of momentum in the fluid. Different strategies can be chosen to perform the averaging of
125
particle-related information to the coarse mesh. This operation is also referred to as backward interpolation.
126
Similarly, multiple choices are possible to interpolate the flow fields at the location of the particles. In
127
this contribution, the procedure proposed in [27] is adopted to perform the backward interpolation, while
128
R functions are adopted for the forward interpolation. Also, the semi-implicit the standard OpenFOAM
129
treatment algorithm for the calculation of the fluid-particle interaction terms of [27] is used.
131
PT
CE
AC
130
ED
121
After the fluid-particle interaction problem is solved, we project the updated variables in a second mesh
where the fluid fine scale is resolved. We will refer to this second mesh as the supporting domain of our
132
approach. This procedure is described in section 6. The fluid fine-scale problem is in this way initialized
133
with particle-related information provided by the coarse-scale solution. In the obtained set of equations,
134
the variables derived from the Lagrange fields projection are defined on an averaged scale. The usage of a
135
supporting domain for the fluid equations solution forces us to adopt those averaged fields, but enables us 7
ACCEPTED MANUSCRIPT
-Mapping Lagrangian fields to Eulerian. -Solving fluidparticles interactions
CR IP T
Fluid Solution
Particle related fields
Coarse Mesh
- Solving Fluid equations
AN US
Fine Mesh
Figure 2: Diagram of the solution procedure of the two different length scales for the simulation. The two boxes represent the different models adopted, while the arrows show schematically the communication between the scales. A coarse mesh (top) is used to map the Lagrangian field of the DEM into an Eulerian reference and to solve the fluid-particle interaction. Particle-related fields are mapped to the supporting domain (bottom) where a finer mesh is used to solve the fluid equations.
to sharply reduce the computational burden. The discretization of the supporting domain can be done in
137
accordance with the fluid system requirements, and in general can be totally independent of the particles’
138
characteristics. The usage of sub-grid modeling at the fluid fine scale is also in theory possible, but goes
139
beyond the purpose of this paper, in which we will assume the effect of scales smaller than the fluid fine one
140
as negligible.
PT
ED
M
136
After the set of equations of the supporting domain has been solved, the fluid variables are mapped to the
142
coarse domain where averaged variables are reconstructed from the up to date data. It must be noted that
143
using a finer discretization for the solution of the fluid equations does obviously increase the computational
144
burden with respect to a classical DEM-VOF simulation. However, CFD-DEM couplings based on the
145
volume averaging technique are normally applied in simulation featuring O(103 ) − O(106 ) particles. In this
147
148
AC
146
CE
141
range, the DEM often represents the more expensive part of the simulation. Therefore, as better shown in section 7.3, the additional cost related to the usage of a finer grid can, in many cases, be negligible. 3. Granular phase governing equations
149
The Discrete Element Method (DEM) is a well-established approach for granular materials [6]. It has
150
been proven to be capable of successfully addressing the physics of granular flows in an efficient way. Its 8
ACCEPTED MANUSCRIPT
151
advantages gave rise to multiple attempts to extend DEM solvers by coupling them with flow fields resolved
152
by CFD [14, 10, 13, 16]. In this work, the DEM is used to evolve a set of particles describing physical entities
153
in the presence of a multiphase fluid. Defining with xi the positions, mi the masses, φi the orientation, and
154
Ii the tensorial moment of inertia of the particles one can write d2 xi dt2 d2 Ii 2 φi dt
155
= Mc + Mext ,
where Fc is the collision force Fc =
X
Fij (xj , up j , φj , ωj ),
i6=j
(2)
CR IP T
= Fc + Ff + Fg ,
mi
P
(3)
(4)
with up the particle velocity, and ω the angular velocity. In particular,
157
particles different from particle i, while Fij is the interaction force between particle i and particle j. At the
158
same time, the term Mc indicates the torque acting on the particle due to collisions Mc =
X
AN US
156
i6=j
represents the sum over all
Mij (xj , up j , φj , ωj ),
i6=j
(5)
with Mij the torque acted from particle j to particle i. A detailed description of all the models present in
160
the literature for the equations 4 and 5 is presented in [28], and citations within. In equation 2 the term Fg
161
represents the gravitational force and Ff takes into account the force rising from the interaction with the
162
fluid. An analytic expression for Ff can be obtained as
PT
ED
M
159
Ff =
Z
S
σf · n dγ,
(6)
with σf the fluid stress tensor, n the normal vector to the particle-fluid interface, S the particle surface.
164
The specific treatment of this term will be addressed in section 5.1.
165
4. Governing equations of the fluid scale
167
168
AC
166
CE
163
In this contribution, we refer to an unsteady incompressible, immiscible, two-fluid system, with forcing
terms arising from the third phase in accordance with section 5.1. The equations governing this system are the Navier-Stokes equations in the form ∂ρf uf + ∇ · (ρf uf uf ) = −∇p + ∇ · µf ∇uf + ∇T uf + TΓ + Fb + Ffpi , ∂t 9
(7)
ACCEPTED MANUSCRIPT
169
with TΓ = ΓCnδΓ the surface tension force, Fb a generic body force, Ffpi the fluid-particle interaction force
170
as shown in section 5.1. Within each phase, the incompressibility constraint must be fulfilled in the form ∇ · uf = 0.
(8)
In diluted flow configurations, the presence of the discrete phase is commonly taken into account by the
172
term Ffpi [1], while in configurations characterized by dense granular phase (like for instance, packed beds)
173
porous equations are normally preferred. If one defines the porosity as = Vsolid /Vcell , the equations take
174
the form
CR IP T
171
∂ρf uf + ∇ · (ρf uf uf ) = −∇p + ∇ · µf ∇uf + ∇T uf + TΓ + Fb + Ffpi , ∂t
(9)
where the terms Fb , and Ffpi , must be changed accordingly [5]. Similarly, the incompressibility constraint
176
imposes that, within a control volume, changes in the void fraction are taken into account leading to
AN US
175
∇ · uf =
ρf (x) = ρ1 α(x) + ρ2 (1 − α(x)),
(11)
µf (x) = µ1 α(x) + µ2 (1 − α(x)),
(12)
ED
178
with α the volume fraction defined by
=
χ =
Z 1 χ(x)dx, V V 1 if first fluid, 0
(13)
(14)
if second fluid.
The volume fraction is considered a scalar transported by the fluid flow for which
AC
180
CE
PT
α
179
(10)
Local density and viscosity are dependent on the fluid phase and can be written in the form
M
177
∂ . ∂τ
∂α + ∇ · (αuf ) + ∇ · (α(1 − α)uc ) = 0, ∂t
(15)
must hold, with uc the relative velocity between the two-phases referred to as compression velocity. The
181
third term is introduced in order to avoid an excessive numerical dissipation.
182
4.1. Finite-volume discretization of the fluid equations
183
R [25], an extensively-used openThe computations of the flow fields are performed using OpenFOAM
184
R structure, the control source library for field operation and manipulation. According to the OpenFOAM
10
ACCEPTED MANUSCRIPT
volumes can be of general polyhedral shape, all depended variables share the same control volumes, and
186
equations are treated in a segregated way. We adopt, and here briefly recall, the VOF scheme described in
187
[4]. In particular, the changes proposed by [10] are used for the porous cases. The finite volume method
188
(FVM) is used to discretize the conservation equations. For the source terms as well as the transient term
189
the midpoint rule is adopted and an integration over the cell volume is performed. Convective and diffusive
190
terms are treated converting the volume integrals into surface integrals (Gauss theorem). The integration
191
is then performed over the surface values. Gradients are finally estimated with a linear face interpolation,
192
with a correction term that takes into account the non-orthogonal contributions. For the convective term,
193
a normalized variable diagram (NVD) based face interpolation is chosen, with a high-resolution differencing
194
scheme proposed in [9], in order to ensure boundedness in the discretization. An Euler implicit scheme is
195
used for the temporal discretization [8]. As proposed in [8], the method expresses the face values required
196
by the discretization of diffusive and convective terms with the new time-level cell values. This results in
197
a first-order accurate method that guarantees boundedness of the solution. In particular, the fluid-particle
198
interaction term is treated in a semi-implicit way according to [27]. As proposed in [4] the term uc of equation
199
15 is evaluated as a conservative volume flux derived from the pressure-velocity coupling. In particular, its
200
expression is
M
AN US
CR IP T
185
204
205
ED
PT
203
normal face flux calculated from the gradient of the volume fraction as ξF =
∇α · AF , ||∇α| + ∆|
(17)
with ∆ a stabilization factor introduced for non uniform grids, evaluated on the mesh characteristic by
CE
202
(16)
with η the face volume flux, kf the free surface compression factor here set to 1. Finally, ξF is the face
AC
201
η , max η , uc,F = ξF min kf |AF | |AF |
∆ = P
B
Nc
Vi
Nc
13 ,
(18)
with B a small parameter (here set to 10−7 ) and N the number of computational cells. A temporal sub-
cycling strategy for the advection of the phase is adopted. This is motivated by the fact that for VOF-based
206
methods the stability of the solution procedure is very sensitive to the phase advection problem. In order to
207
address this problem, in addition to the usage of bounded discretization schemes for the divergence terms,
208
the solution of the volume fraction equation is divided into multiple sub-cycles. The time step of the single 11
ACCEPTED MANUSCRIPT
209
sub-cycle is obtained as δTi =
δTtot , Nsc
(19)
with Nsc the number of sub-cycles, δTtot the simulation time step. Within a sub-cycle, the volume fraction
211
α is calculated, and from this the total mass flux Φm is obtained as ΦG m = ρf uf · AF =
Nsc X δTi j Φ , δ tot m j=1 T
212
with Φjm the partial mass flux, Nsc the number of sub-cycles.
213
4.2. PIMPLE loop for the porous equations
CR IP T
210
(20)
R , applied on the We here briefly describe the standard PIMPLE loop as implemented in OpenFOAM
215
porous equation 9. First, the porosity field is mapped from the coarse mesh into the centroids points of the
216
fine mesh, and the porosity values on the cell faces are reconstructed. Then, the matrix for the velocity
217
equation is assembled, and the momentum predictor is solved in the form
AN US
214
au P ufP = H − ∇p − gh∇ρf + ΓC∇α + Ffpi , where au P is the matrix coefficient related to the centroid point, and
M
218
(21)
ED
H = ∇ · µf ∇uf + ∇T uf
−
X
au N ufN ,
(22)
N
with aN the coefficients for the neighbor centroids. The velocity equation is relaxed to ensure stability and
220
solved accordingly. The pressure-correction is then performed as follows. From equation 21, a predicted
221
velocity is obtained as
PT
219
CE
−1 ufP = [au [H − ∇p − gh∇ρf + ΓC∇α + Ffpi ] , P]
(23)
and the face values of the velocity ufF is estimated through face interpolation. The pressure equation is
223
then solved as
AC
222
∂ −1 ∇ · 2 [au ∇p = ∇ · uf + , P] ∂t
(24)
224
that is the discretized version of equation 10, where the time derivative of the porosity fields takes into
225
account the changes of the cell volume due to particle presence. The velocity field is then corrected with
226
the resolved pressure field, so to obtain conservative fluxes, and therefore a divergence-free velocity field.
12
ACCEPTED MANUSCRIPT
227
5. Fluid particle system
228
5.1. Momentum coupling
230
As a standard procedure in fluid dynamics, the fluid stress tensor in equation 6 can be divided into pressure and residual strain contribution, leading to Ff = −
231
Z
S
p I · n dγ +
Z
S
τf · n dγ = Flif t + Fdrag ,
CR IP T
229
with τfij = σfij + pδij ,
(25)
(26)
and p = −1/3 tr(σf ) the fluid pressure, δij the Kronecker delta, I the identity matrix. If the fluid-particle
233
interaction is aimed to be solved on an averaged scale, it is possible to consider the pressure gradient of a
234
region of fluid surrounding the particle as constant. In this scenario, it is convenient to express the first
235
integral as Flif t = −
Z
S
AN US
232
p I · n dγ =
Z
Vp
∇p dV,
(27)
where the Green theorem has been used, Vp is the volume of the particle. Considering the pressure gradient
237
as constant within the integration volume, one obtains
M
236
ED
Flif t = −∇p Vp ,
(28)
with p pressure of the liquid. In a gas-liquid-particles system, the pressure gradient between the liquid and
239
the gas region is normally not negligible due to the difference in density typical of a gas-liquid configuration.
240
This effect must, therefore, be taken into account.
243
244
CE
242
A well-stated method to study the transport of particulate matter into a fluid flow is via the usage of a drag law that takes into account the average effect of the fluid stress on the solid
AC
241
PT
238
Fdrag =
Z
S
τf · n dγ .
(29)
Many authors have focused on the description of the drag force [28, 20] that is normally provided in the form of a semi-empirical law Fdrag = β(uf luid − up ),
(30)
β = β(uf luid − up , ρf , ρp , dp , Ap , µf , ),
(31)
13
ACCEPTED MANUSCRIPT
245
with uf luid , uparticle the fluid and particle velocity respectively, ρf , ρp the respective densities, dp , Ap the
246
particle characteristic length and surface, µf the fluid viscosity. For the sake of generality, we took β as in
247
[19]
Cd =
0.5
d2p 4 π Cd ρf
24(1+0.15Re Re
0.687
)
|up − uf | ,
(32)
if Re < 1000
CR IP T
β=
0.44
(33)
if Re > 1000.
It has to be pointed out that, in the empirical correlation for the drag force of equation 32, the pressure
249
contribution cannot be separated as proposed in equation 25. A contribution caused by the pressure dis-
250
tribution around the particle due to creeping flow is taken into account in those models. However, in a
251
volume-averaged coupling this interaction is not resolved due to the choice of the spatial coarse scale. At
252
this scale, the main assumption is that the particles can be considered as zero-dimensional entities, and
253
therefore the effects of the creeping flow are modeled and not resolved. Therefore, equation 25 holds on the
254
coarse scale since here the fluid variables have to be seen as bulk pressure and bulk velocity. A further as-
255
sumption consists of considering the angular momentum exchange between particles and fluid as negligible.
256
The results of this work focus on spherical particles. However, we want to clarify that in case of particles
257
with a pronounced nonsphericity, the angular momentum exchange between fluid and particle should not
258
be neglected.
259
5.2. Phase interchange strategies
AC
CE
PT
ED
M
AN US
248
Fluid 2
Fluid 2
α=0
α=0
α=1
α=1
Fluid 1
Fluid 1
Figure 3: Particle crossing the interface between two fluids.
260
As previously mentioned, the VOF (and in general the interface capturing methods) can be extremely
261
interesting for a CFD-DEM coupling, as it uses a single velocity and single pressure formulation. This allows 14
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 4: Phase replacement technique, particle selected for conversion to equivalent droplet. Coarse grid in red, fine grid in black.
262
simplifying the model adopted for the fluid-particle coupling since for the most part of the computational
263
domain (the one not interested by the presence of the interface) it will be reduced to a single-phase coupling.
264
Nevertheless, the interaction between particle and interface remains a complex phenomenon that needs to
265
be addressed.
First, as shown in figure 3, when a particle crosses the interface between a Fluid 1 and a Fluid 2 it will
267
not occupy volume in the Fluid 2 but instead a volume in Fluid 1. This phase interchange phenomenon
268
must be taken into account and is particularly important in dense granular flows. Then, the particle motion
269
will induce a perturbation of the interface. This phenomenon is particularly interesting to be studied
270
since a local perturbation of the interface can have a major influence on the global two-phase configuration.
271
Finally, capillary forces arise at a gas-liquid interface that can influence the particle dynamic. At a gas-liquid
272
interface, a difference of pressure will be established in the form
274
ED
PT
CE
∆p = ΓC,
(34)
with Γ the surface tension, C the interface local curvature. The capillary force will act in the normal
AC
273
M
266
direction to the surface, namely Fcap = ΓC n.
(35)
275
As mentioned in section 2, this work focuses on high Stokes number flows for which the effects of capillary
276
forces are assumed being negligible with respect to the inertia of the particles. This is reasonable when
277
operating with common values of surface tension (indicatively below 100 (10−3 N )/m). 15
ACCEPTED MANUSCRIPT
278
We here propose three different strategies for handling the phase interchange with reference to the
279
multiscale DEM-VOF method. In [10] the authors propose a porous approach for the phase interchange
280
problem in a single-scale DEM-VOF coupling. This approach consists of modifying the volume fraction
281
transport equation in order to take into account the presence of the granular phase as a porosity field,
282
leading to
CR IP T
∂α + ∇ · (αuf ) + ∇ · (α(1 − α)uc ) = 0. ∂t
(36)
The authors showed that the method can provide good results in case of packed beds of solid particles
284
interacting with a steady flow under the effect of gravity. This approach can be seen as a natural extension
285
to multiphase environments of the classic volume averaging technique for dense granular phases. As shown
286
in [10] for a sedimentation problem, this approach leads to errors in the volume conservation in the order of
287
the classical algebraic two-phase VOF on a coarse grid. This method has an obvious limitation in case of
288
high porosity (dilute) configurations since the interface motion will tend to be not affected by the particle
289
presence. In a multiscale DEM-VOF coupling, this phase interchange strategy can be used as well. The
290
porosity field defined on the coarse scale will be used as well in the modified advection equation 36. In
291
addition, one can expect the VOF scheme to be less dissipative when operated on a finer grid [7].
M
AN US
283
Another possible phase interchange strategy purely based on the definition of the color function consists
293
of re-defining the standard two-phase color function while keep using equation 15 for the phase advection.
294
In a two-phase environment, the color function represents the ratio
ED
292
Vliquid , Vtot
(37)
Vliquid . Vliquid + Vgas
(38)
297
CE
296
or
αtwo phase =
When particles enter in a two-phase environment, the total volume can be expressed as
AC
295
PT
αtwo phase =
Vtot = Vliquid + Vgas + Vparticles .
(39)
If one keeps operating on the coarse scale one can define a particle volumetric concentration as ψ=
Vparticles . Vtot
16
(40)
ACCEPTED MANUSCRIPT
298
For gas-liquid-particles mixtures, it can be useful to define a heavy-phase volume concentration as αthree phase =
299
Vliquid + Vparticles . Vliquid + Vgas + Vparticles
(41)
From the previous definitions, one can easily obtain the relation αthree phase = αtwo phase (1 − ψ) + ψ.
CR IP T
(42)
300
In this way, the volume occupied by the particle is substituted by a ghost fluid that enters directly in the
301
phase advection problem. The phase advection equation 15 will remain unchanged, while the variable itself
302
is substituted. As shown in section 7.2, this approach is extremely sensitive to the grid resolution.
Finally, using a dual-grid multiscale algorithm for the DEM-VOF coupling allows adopting a modified
304
version of the previous strategy that will be referred to as phase replacement technique. This method is
305
based on the replacement between DEM particle and an equivalent droplet. In [12] the authors, study
306
atomization processes with a VOF technique. In order to save computational cost, the authors propose to
307
substitute small droplets with Lagrangian point-particle entities, when they are far from the interface. The
308
idea we propose here is to substitute the DEM particle with a liquid droplet of the same volume when it
309
impacts the interface. In particular, as shown in figure 4, particles that are in a part of the domain interested
310
by the presence of the interface are selected for possible conversion. For those, the relative velocity between
311
the interface and the particle is compared to a reference parameter in the form
ED
M
AN US
303
(uint − up ) · nint > δC ,
313
with nint the interface normal pointing in the direction of the phase in which the particle currently is, i i C δC ∈ 0, δ∆tot , ∆C the coarse grid smallest dimension. Only the particles for which equation 43 holds will T
PT
312
(43)
be converted. After the impact, when a particle enters the second phase, it can again be described as a
315
DEM entity in its motion. As pointed out in section 1, this is one of the main advantages of the VOF
316
technique regarding the coupling with DEM. This approach allows switching between a zero-dimensional to
318
319
AC
317
CE
314
a three-dimensional momentum coupling, with a competitive computational cost. In particular, as described in section 5.1, the momentum coupling between DEM entity and fluid is zero-dimensional (based on the drag force formulation), but the interface is perturbed by a three-dimensional entity (the droplet). This technique
320
both allows implicitly taking into account the volume conservation, and modeling a three-dimensional DEM
321
entity-fluid interaction when the higher accuracy is required (i.e when the interface is perturbed by the
322
particles impact). It is worthy to point out that, in the current contribution, there is no density difference 17
ACCEPTED MANUSCRIPT
between the equivalent droplet and the heavy fluid phase. A local modification of the density of the fluid is
324
theoretically possible, but can potentially create numerical instabilities. The density difference between the
325
particle and the fluid phase will, on the other hand, be taken into account in the calculation of the particle
326
trajectory, in accordance with the standard DEM solution. As shown in the example 7.2, this technique is
327
only applicable when the supporting domain provides a fine discretization nearby the surface to properly
328
reconstruct the equivalent droplet in a VOF scheme. Therefore, the additional benefit of treating the impact
329
phenomena in a three-dimensional way is only possible in a multiscale approach.
330
6. Interpolation strategy
CR IP T
323
The dual-grid approach described in section 2 requires the definition of interpolation strategies between
332
the coarse and the fine mesh. This increases the algorithm’s complexity and requires a thoughtful considera-
333
R libraries tion. In this work, we rely on the meshT oM esh interpolation class available in the OpenFOAM
334
that are here used to map fields from the coarse to the fine mesh, and vice versa fine mesh related fields
335
R provides three different interpolation into the coarse domain. At its default implementation, OpenFOAM
336
strategies. However, it is possible to extend the existing interpolation techniques of the meshT oM esh class.
337
R are The three default strategies present in OpenFOAM
M
AN US
331
• meshT oM esh :: M AP that uses a direct mapping of the nearest cell-value, if the source mesh is
339
uniformly discretized (as in most of the case for CFD-DEM couplings) this mapping will preserve the
340
volume fraction and source of momentum. However, the resulting mapped fields will not be smooth
343
PT
342
• meshT oM esh :: IN T ERP OLAT ION uses a (gradient based) inverse distance interpolation. This is not, in general, a linear interpolation, and one can expect it not to have second-order accuracy. • meshT oM esh :: CELL P OIN T IN T ERP OLAT E performs a linear interpolation (preserving second-
CE
341
ED
338
order accuracy) by creating a tetrahedron with the center point of the cell containing the target and
345
three vertexes of the cell, and then uses an inverse distance weights to perform the linear interpolation.
346
347
AC
344
This strategy can nevertheless create problems on boundaries since the boundary value will excessively effect near-wall values.
348
A fourth strategy called CELL P OIN T F ACE IN T ERP OLAT E has been added in the class for this
349
contribution. It is still based on the construction of tetrahedron to perform linear interpolation, but executes
350
an additional check on the boundary. This allows avoiding the above-mentioned problem. With a slight
351
modification of the meshT oM esh class, it is possible to add this interpolation strategy to the options for 18
ACCEPTED MANUSCRIPT
mesh based interpolations. For the sake of simplicity, no other interpolation strategy is considered in this
353
work. In particular, if not differently specified a meshT oM esh :: M AP strategy will be used to map
354
porosity and source of momentum (Ffpi ) field from the coarse to the fine mesh, while a meshT oM esh ::
355
CELL P OIN T F ACE IN T ERP OLAT E will be used to map velocity and pressure fields from the fine
356
mesh to the coarse one. A thoughtful investigation of the optimal mapping technique for the multiscale
357
CFD-DEM method is however of high importance and will be subject for upcoming studies.
358
7. Validation and Results
359
7.1. Three-phase dam break
M
AN US
CR IP T
352
ED
Figure 5: Three-phase dam break. Setup of the simulation.
The dam break is a famous benchmark for two-phase flows simulations. The extension to a gas-liquid-
361
particle configuration represents an interesting challenge. Similar versions of the so-called “three-phase dam
362
break” have been proposed in [20, 21]. In both contributions, the authors propose a dam-break configuration
363
shown in figure 5, without intermediate obstacle, in a box of dimensions 0.2m x 0.1m x 0.3m. A column
364
of water of extension 0.05m x 0.1m x 0.1m with a uniformly distributed bed of spherical particles at the
365
bottom, breaks and stabilizes. The spheres have a diameter of 2.7 mm. We propose the multiscale simulation
366
of the benchmark using a supporting domain discretized with 384k identical cubic cells. The main point
368
CE
AC
367
PT
360
that made us choose this benchmark is that in [20, 21] the authors provide both experimental and numerical results for the test case. This makes it extremely appealing for testing a new method since advantages and
369
drawbacks of the algorithm can be directly tested. The aim is to show how, in this interface problem, the
370
usage of a finer discretization for the fluid domain can add significant advantages in the problem solution.
371
In particular, the dual-grid approach can produce grid convergent results for the continuum phase and
372
correctly capture some additional integral (kinetic energy), and local (shape of the interface) properties of 19
ACCEPTED MANUSCRIPT
373
the fluid. We adopt simulation parameters according to the authors’ suggestions as follows. The liquid
374
(heavy phase) density and viscosity are taken as 1000kg/m3 and 10−3 P a s respectively, while for the gas
375
(light phase) 1 kg/m3 and 10−5 P a s are chosen. The particle density is fixed ad 2500 kg/m3 . A linear
376
dashpot impact model is adopted for both particle-particle and particle-container interactions with spring
CR IP T
constant of 1000 N/m, a restitution coefficient of 0.9 and a friction coefficient of 0.3, as proposed in [20]. In
Liquid and solid front position as a function of Time 4 3.5
Liquid Experiments (Sun et al) Liquid MultiscaleDEM-VOF Solid Experiments (Sun et al) Solid MultiscaleDEM-VOF
AN US
Z* coord
3 2.5 2
1
0
0.5
1
M
1.5
1.5
2
2.5
3
3.5
4
ED
Nondimensional time
PT
Figure 6: Three-phase dam break. Comparison between the solution obtained with the multiscale approach and experimental results relative to the solid and liquid front motion. 377
figure 6, results related to the fluid and solid wave-front positions calculated with the multiscale DEM-VOF
379
method and the experimental results of [21] are proposed. The simulation is performed with a coarse domain
380
discretization of 6k identical cubic cells and a fine domain discretization of 384k cells. Results show a very
381
good agreement for both particles and fluid wave-front.
383
384
A second comparison with the experimental and
AC
382
CE
378
numerical results presented in the literature [20] is shown in figure 7. The discretization of the coarse and fine scale is chosen as 6k and 384k identical cubic cells respectively. It can be observed how the usage of
a finer discretization for the supporting domain (that is easily obtained in a multiscale approach), leads
385
to a better interface reconstruction. In particular, some perturbations of the gas-liquid interface can be
386
captured and the shape of the interface can be reconstructed with a higher resolution. In order to quantify
20
CR IP T
ACCEPTED MANUSCRIPT
ED
M
AN US
Figure 7: Three-phase dam break. Comparison between the solution obtained with the multiscale approach (left), the experimental setup (middle), and the DEM-VOF method [20] at time t = 0.3s , 0.4s.
CE
PT
Source Experiment Sun et al Relative error% Multiscale Relative error%
x∗A 0.887 0.804 9.4 0.882 0.56
x∗B 0.837 0.783 6.7 0.815 1.2
∗ yA 0.162 0.141 13.0 0.160 1.23
∗ yB 0.243 0.194 20.2 0.250 2.9
∗ yC 0.678 0.851 25.5 0.7001 1.3
387
AC
Figure 8: Three-phase dam break, comparison through sampling points. Points of the surface A, B, C, chosen for the comparison (top), and table comparing the nondimensional coordinates of A, B, and C. A reduction of the relative error is observed.
the improvements, we define a set of nondimensional coordinates x∗ =
388
x ∗ y ,y = , l0 l0
(44)
with l0 the length of the computational domain (0.2 m). Both the multiscale DEM-VOF and the DEM-VOF
21
ACCEPTED MANUSCRIPT
method are able to qualitatively reconstruct the shape of the interface at the left border of the domain at
390
time 0.4 s, as shown in figure 8. In order to obtain a quantitative comparison three characteristic points
391
of the curve are chosen. In particular, A and B are respectively the maximum and minimum point of the
392
first and second curve in the x − z plane, with respect to the coordinate z. C is chosen as the first point
393
at which the liquid-gas interface touches the wall. Figure 8 compares the nondimensional coordinates of
394
points A, B, and C. From the comparison it follows that a multiscale approach significantly reduces the
CR IP T
389
relative error with respect to the experimental results. In figure 9, a comparison between the nondimensional
Kinetic energy of the fluid as a function of time 0.35
AN US
48k cells(single scale) 162k cells 384k cells 500k cells
0.3 0.25 0.2
M
0.15 0.1 0.05 0
2
PT
0
ED
Nondimensional kinetic energy
0.4
4
6
8
10
12
14
Nondimensional time
395
CE
Figure 9: Three-phase dam break. Comparison between the nondimensional kinetic energy of the fluid phase obtained with fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of this quantity is observed with discretization finer than 162k cells.
total kinetic energy of the fluid phase obtained with different domain discretization is proposed. We define
397
a nondimensional kinetic energy as the ratio between the kinetic energy and the potential energy of the
398
399
AC
396
column of water at the beginning of the simulation R u · u ρ 1 dV R f f f2 Ekin = D . A ρ gdZ Zcol in h
(45)
In particular, a fine-scale discretization of 48k, 162k, 384k, and 500k identical cubic cells is adopted. In
22
ACCEPTED MANUSCRIPT
400
figure 9, the nondimensional time is defined as ∗
t =t
2g a
21
,
(46)
with a = 0.05m as in [20]. It is worthy to be noted that a structured Cartesian grid is here adopted solely
402
for the sake of simplicity since it is not a requirement for the algorithm proposed. It can be observed how
403
the kinetic energy converges with grid resolutions smaller than the particle size. Therefore a single-scale
404
approach based on the volume averaging technique, for which the finest grid-discretization is of 48k cells,
405
cannot correctly capture this quantity. In particular, the error on the kinetic energy appears higher when
406
the kinetic energy reaches local maxima (corresponding to the impact of the liquid front on the walls). This
407
can be explained considering that when higher velocities occur in the fluid, the interface dynamic becomes
CR IP T
401
0.0006 0.0004 0.0002 0
0
1
2
3
4
5
6
0.0006
500k cells 384k cells 162k cells 48k cells
0.0005 0.0004 0.0003 0.0002 0.0001 0
0
1
2
3
4
5
6
Nondimensional time
ED
Nondimensional time
M
L2 Norm of U [m5/s2]
500k cells 384k cells 162k cells 48k cells
0.0008
L2 Norm of the velocity as a function of time
L2 Norm of U-X-component [m5/s2]
L2 Norm of the velocity as a function of time 0.001
AN US
more complicated and therefore a finer discretization is required to solve them. In figure 10, the L2 norms
Figure 10: Three-phase dam break. Comparison between the L2 Norm for the velocity magnitude (left) and x component (right) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of these quantities is observed with discretization finer than 162k cells.
of the velocity vectors and the velocity x component
CE
409
PT
408
411
412
kuf x kL2 =
Z
ZD
uf · uf dV,
(47)
uf x uf x dV,
(48)
D
AC
410
kuf kL2 =
are proposed. It is interesting to notice how the kinetic energy of the flow is underestimated by the usage of
coarser grids, while the norm of the velocity field is overestimated. This can be explained considering that the two-phases have an high density ratio. Therefore, a difference in the volume fraction can induce very
413
significant changes in the global kinetic energy. A further investigation of this phenomenon is proposed in
414
figure 11, where the L2 norms of the velocity kuliquid kL2 =
Z
D
αuf · uf dV,
kugas kL2 = 23
Z
D
(1 − α)uf · uf dV,
(49)
ACCEPTED MANUSCRIPT
L2 Norm of Uliquid as a function of time
L2 Norm of Ugas as a function of time 0.0009
0.0004
48k cells(single scale) 162k cells 384k cells 500k cells
0.0008
48k cells(single scale) 162k cells 384k cells 500k cells
L2 Norm of Ugas [m5/s2]
L2 Norm of Uliquid [m5/s2]
0.0005
0.0003 0.0002 0.0001
0.0007 0.0006 0.0005 0.0004 0.0003 0.0002
0
1
2
3
4
5
6
0
7
0
1
Nondimensional time
CR IP T
0.0001 0
2
3
4
5
6
7
Nondimensional time
Volume of liquid in interface cells as a function of time 48k cells(single scale) 162k cells 384k cells 500k cells
0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0
0
1
AN US
Volume of liquid [m3]
0.0007
2
3
4
5
6
7
Nondimensional time
M
Figure 11: Three-phase dam break. Comparison between the L2 Norm for the velocity of the liquid phase (top-left), and gas phase (top-right), and volume of liquid occupying an interface cell (bottom) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of these quantities is observed with discretization finer than 162k cells.
are plotted as a function of time. One can observe marked differences in the gas-phase prediction that, in a
416
single-scale approach, has a peak with almost double the value of a grid convergent result. Also, in figure 11
417
the volume of the liquid occupying a cell with α ∈]1, 0[ is proposed. This shows how the adoption of a coarse
418
grid induces almost double of the fluid to be contained within interface cells. This can cause erroneous
419
predictions of the fluid behavior, as the fluid in the interface cells is treated as an equivalent mixture, and
420
therefore the solution is more sensitive to the adopted model.
421
and the volume occupied by the solid phase as functions of time are computed with different fine domain
422
discretization
AC
CE
PT
ED
415
Z
kpkL2 = pp dV, D Z Vsolid = (1 − ) dV.
In figure 12, the L2 norm of the pressure
(50) (51)
D
423
Negligible differences are observed in the pressure norm. No mesh dependency is observed for the volume
424
occupied by the solid, this shows that the projection of the porosity field from the coarse grid to the fine
425
is not sensitive to the discretization of the fine grid, as was expected due to the choice of the mapping 24
ACCEPTED MANUSCRIPT
Volume occupied by the solid as a function of time 42
500k cells 384k cells 162k cells 48k cells
80 70
41.8
60
Volume [cm3]
50 40 30
41.6 41.4
20
500k cells 384k cells 162k cells 48k cells
41.2
10 0
0
1
2
3
4
5
41
6
0
1
Nondimensional time
CR IP T
L2 Norm of the pressure [Pa2 m3]
L2 Norm of the pressure as a function of time 90
2
3
4
5
6
Nondimensional time
Figure 12: Three-phase dam break. Comparison between the L2 Norm of the pressure (left) and the volume occupied by the solid phase (right) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Minor grid dependency is noticed.
strategy. Obtaining grid convergent solutions is a very important characteristic of a numerical method,
427
especially for applications in which a complete experimental validation of the quantities of interest is not
428
possible. Therefore, this makes the proposed method competitive over a single-scale DEM-VOF method,
AN US
426
as it allows this fundamental numerical verification.
In a dual-grid approach, the mesh independence for
Liquid and solid front position as a function of Time
M
4 3.5
ED
Solid Experiments (Sun et al) Coarse 6000 Coarse 26000 Coarse 32000
PT
2.5
*
Z coord
3
CE
2 1.5
AC
1
0
0.5
1
1.5
2
2.5
3
3.5
4
Nondimensional time
Figure 13: Three-phase dam break. Comparison between the position of the particle front obtained with coarse domain discretization of 6k , 26k and 38k cells. Negligible differences are observed. 429
430
the coarse grid must also be addressed. In figure 13, the result described in figure 6 are shown with three
431
different coarse meshes of 6k, 26k, and 38k cells respectively. Negligible effects of the coarse grid choice 25
ACCEPTED MANUSCRIPT
Kinetic energy of the fluid as a function of time 0.35
6k cells(C) 26k cells(C) 38k cells(C)
0.3 0.25 0.2 0.15 0.1 0.05 0
0
2
4
6
8
10
12
14
AN US
Nondimensional time
CR IP T
Nondimensional kinetic energy
0.4
Figure 14: Three-phase dam break. Comparison between the nondimensional kinetic energy of the fluid phase obtained with coarse domain discretization of 6k , 26k and 38k cells. (Coarse) Mesh independence of this quantity is observed.
on the position of the front can be observed. This suggests that the global behavior of the particles can be
433
captured with grid discretization coarser than the one used for the classical single-scale approach regarding
434
the pure CFD-DEM coupling. In figure 14, the influence of the coarse grid choice on the kinetic energy of
435
the fluid is investigated. Three configurations with a fine-scale discretization of 160k, and the previously
436
mentioned coarse grids are confronted. As shown in figure 14, also the total kinetic energy of the fluid
437
expresses mesh independent behavior for a mesh coarser than the mesh used for the single-scale approach.
438
Being able to discretize the coarse grid with a reduced number of cells allows reducing the cost associated
439
with the calculation of the particle-related averaged fields. This can represent a significant part of the
440
cost associated with an Euler-Lagrange coupling, and is an additional important advantage of the dual-grid
441
approach. We want to remark that this dual-grid approach introduces the issue of coarse grid independence.
442
Even if the presented and other cases show coarse grid independence for grids spacing larger than the one
443
used in a single-grid approach, this property must be verified case by case. In conclusion, the ability to
445
ED
PT
CE
AC
444
M
432
separate the solution of the particle-transport scale from the one of the fluid equations is one of the main advantages of the method. In particular, it allows obtaining grid convergent solutions for the continuum
446
phase and at the same time to reduce the computational cost associated with the averaging procedure.
447
7.2. Interface perturbations induced by a single particle
448
A simple three-phase flow configuration is here proposed. It consists of a single particle falling under
449
the effect of gravity into a container partially filled with water. The volume averaging technique is normally 26
AN US
CR IP T
ACCEPTED MANUSCRIPT
(a) Top view Coarse grid (black) and supporting domain grid (red)
(b) Lateral View
M
Figure 15: One Particle falling into Water Simulation setup.
adopted to approach problems characterized by an elevated number of particles. Still, a simple setup like
451
this one can be useful in order to validate the fluid-particle coupling. We use this test case to show how the
452
multiscale approach proposed in this work is able to overcome some of the most important drawbacks of the
453
volume averaging technique concerning the DEM-VOF coupling. The domain dimensions are 0.05 m x 0.05 m
454
x 0.20 m, with the gravity in the -z direction and the initial water level at 0.051m. The spherical particle
455
has a diameter of 4mm, and is initially at rest with its center located at (0.0275m,0.0275m,0.14m). The
456
particle density is fixed at 2000kg/m3 . No-slip boundary conditions are imposed at all boundary except for
457
the top face where a total pressure boundary condition is imposed. The properties of the gas and the liquid
458
phase are chosen as the one of air and water at 25◦ C, respectively. As shown in figure 15a, for the solution
460
PT
CE
AC
459
ED
450
of the fluid-particle interaction a coarse grid discretization of 10 x 10 x 40 cubic cells has been adopted. A volume adding technique described in section 5.2 is adopted. Results are provided for the benchmark using
461
different supporting domain discretization. The aim is to show how the usage of a proper discretization for
462
the supporting domain along with the approach described in section 5.2 is able to capture interface dynamics
463
that are not resolved with a standard DEM-VOF approach. In figure 16, some significant snapshots of the
464
problem solved with a domain discretization of 500k cubic cells are presented. The particle falls in the gas 27
(b)
(c)
ED
M
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(d)
(e)
(f)
PT
Figure 16: One Particle falling into Water. From top to bottom and from left to right significant snapshots of the problem at time t = 0.0s (a), t = 0.13s (b), t = 0.17s(c), t = 0.18s(d), t = 0.21s(e), t = 1.0s(f), solved with a supporting domain discretization of 500k cubic cells.
phase under the effect of the gravitational force, and after t = 0.14s impacts against the liquid-gas interface
466
generating a crater structure. Around t = 0.17 the particle penetrates the interface and continues falling
467
across the liquid column with reduced velocity. The gas-liquid interface stabilizes after forming a central jet
469
470
AC
468
CE
465
at t = 0.21s. After t = 1s the particle is at rest at the bottom of the domain and the interface recovers a smooth shape (except from secondary perturbations). In figure 17a, we propose a comparison between the gas-liquid interface reconstructed at t = 0.18s with four different supporting domain discretization. The first
471
choice for the supporting domain consists of a grid with 10 x 10 x 40 identical cubic cells. This represents the
472
limit case for the multiscale approach in which the coarse domain and the supporting domain are identically
473
discretized. In this condition, the approach is equivalent to a standard DEM-VOF, since no scale separation 28
ACCEPTED MANUSCRIPT
Interface Comparison with different mesh resolution
Interface Comparison with different mesh resolution 0.07
0.054
840k 500k 108k 4k
0.065
0.05
Z-coordinate [m]
0.048 0.046 0.044 0.042
840k 500k 108k 4k
0.04 0.038
0.055 0.05
0.036 0
0.06
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
0.045
0
CR IP T
Z-coordinate [m]
0.052
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
(a) time t = 0.18s
(b) time t = 0.21s
AN US
Figure 17: One Particle falling into water. Interface comparison between a supporting domain discretization of 864k, 500k, 108k and 4k cells respectively.
is introduced. It is possible to observe how this resolution is not able to correctly capture the position of the
475
interface. This is easily explained by considering that the grid spacing is equivalent to 0.005m. As recalled
476
in section 4 the VOF schemes smear the interface in few cells, therefore the uncertainty on the surface
477
collocation will be proportional to the grid size. Using very coarse grids for the phase advection solution
478
naturally leads to high uncertainty on the interface position that, in general, can make the algorithm poorly
479
conservative. Using a finer discretization of 30 x 30 x 120 cells, the error in the interface level reduces
480
significantly, but still, the local reconstruction of the interface is not able to capture the formation of the
481
crater. The discretization of 60 x 60 x 240, on the other hand, succeeds in reproducing the stretch of the
482
interface surface with the formation of a crater structure (time 0.18 s) and a Rayleigh/Worthington internal
483
jet(time 0.21 s). In figure 17b, a comparison between the gas-liquid interface reconstructed at t = 0.21s with
484
the previously described discretization is proposed. Conclusions similar to those of figure 17a one can be
485
drawn for what concerns the first two meshes. Some more significant differences can here be noticed between
486
the 50 x 50 x 200 and the 60 x 60 x 240 discretization. In particular, the proposed discretization is not
487
sufficient to reach convergence for what concerns the simulation of the Rayleigh/Worthington internal jet.
488
This can be explained considering that the characteristic length of the jet and of the droplet detaching from
490
ED
PT
CE
AC
489
M
474
it are smaller than the one of the crater. Therefore, a finer discretization is required from the VOF scheme to properly reconstruct them. Additionally, since the volume fraction field is modified on the coarse coupling
491
scale, the resulting scheme is highly mesh sensitive. The high mesh sensitivity is a major drawback of the
492
volume adding technique. Finally, since the particle influence on the fluid velocity field is taken into account
493
as a volumetric source in the Navier-Stokes equations, it is not always able to model the correct dynamic
29
ACCEPTED MANUSCRIPT
494
at the two-phase interface. In particular, when the magnitude of the relative velocity between particle and
495
fluid is big enough to produce a significant drag force, the volume adding technique can capture surface
496
perturbations. On the other hand, in cases of reduced relative velocity in the proximity of the surface, a
497
more accurate method may be preferred. In order to show the advantages of the last method proposed in section 5.2 (the volume-replacement
499
technique), we adopt a modified version of the benchmark where the particle is positioned at the beginning
500
of the simulation slightly above the interface with an assigned initial velocity of 0.5m/s in the negative
501
z direction (this represents approximately half of the impact velocity of the previous case). The fluid is
502
considered initially at rest. This is to mimic a sudden interface crossing phenomenon that creates a sharp
503
discontinuity in the flow. As described in section 5.2, a phase replacing technique that takes full advantage
504
of the dual-grid approach, can switch between a zero-dimensional and a three-dimensional interaction with
505
the fluid when higher accuracy is required. This also allows studying the influence of the particles on the
506
fluid motion with drag-force based approaches only when the particles are actually within a single phase. In
507
figure 18, the initial configuration, as well as the streamlines of the gas phase during the interface crossing
508
are displayed. The aim is to show the flow field structures characterized by vortexes in the proximity of
AN US
CR IP T
498
In figure 19a, the crater structure formed after
CE
PT
ED
M
the particle-ghost droplet after the conversion happened.
509
510
511
AC
Figure 18: One Particle Injected into Water. Streamlines of the velocity field in the gas phase at time t = 0.0s, t = 0.01s, t = 0.03s resolved with a supporting domain discretization of 500k cubic cells.
0.03s from the particle injection is proposed. In particular, one can observe that by adopting a phase replacing technique, differences in the reconstruction of the crater obtained using a 108k and a 500k fine
512
grid discretization are negligible. Moreover, in figure 19b, it can be observed how the iso-surfaces at 0.1 of
513
the volume fraction, significantly differ when a coarser discretization is used. This is shown to remark once
514
more the introduced modeling error when adopting a VOF scheme with a coarse grid. Even if in figure 19a
515
one can argue that the 4k grid resolution is somehow able to capture the interface position, in figure 19b is 30
ACCEPTED MANUSCRIPT
Interface Comparison with different mesh resolution
Iso Surface Comparison with different mesh resolution
0.052
0.056
0.051 0.054
Z- coord
Z- coord
0.049 0.048 0.047
500k cells 108k cells 32k cells 4k cells
0.046 0.045
0.05
500k cells 108k cells 32k cells 4k cells
0.048
0.044 0.043
0.052
0.046 0
0.01
0.02
0.03
0.04
0.05
0
Y-coordinate[m]
CR IP T
0.05
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
(a) Interface position
(b) Iso-surfaces α = 0.1 of the phase indicator
AN US
Figure 19: One Particle Injected into Water. Interface comparison between a supporting domain discretization of 500k, 108k, 32k, and 4k cells respectively at time t = 0.0.03s.
516
clearly visualized how coarse grid discretization are smearing the interface region over a significant part of
517
the simulation domain. This is in contrast to the very purpose of adopting an interface-capturing method
518
within a CFD-DEM coupling.
519
7.3. Multiple particles injection in quiescent fluid
In order to show the ability of the method to deal with the simulation of particle-laden multiphase flows
521
characterized by an elevated number of particles, we here propose a benchmark originally introduced in [10].
522
The aim of this is to show how the multiscale approach can overcome some of the DEM-VOF drawbacks
523
while still being suitable for handling simulations involving an elevated number of particles. The benchmark
524
consists of the simulation of two fluids at rest like the case from section 7.2, into which multiple particles
525
are injected causing a movement of the gas-liquid interface position. A simple configuration is proposed
PT
ED
M
520
Initial Level
AC
Case
CE
as described in [10] and shown in figure 20. Particles are introduced into the heavy phase and create a
1
0.3
2
0.3
Table 2: Parameters for the Multiple particles injection tests
Injection type
Injection duration
Continuous 1 s (1500p/s) 3000 par- Initial ticles Condition
dp
ρp
Coarse Mesh edge
Container size
0.005m
2600kg/m3 0.02m
0.1 x 0.1 x 0.5 m
0.004m
2600kg/m3 0.02m
0.1 x 0.1 x 0.5 m
526 527
porous zone at the bottom due to the effect of gravity. The inserted particles will occupy part of the volume
528
previously occupied by the pure liquid and will induce a rise of the heavy fluid’s level. The comparison 31
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 20: Multiple particles injection in quiescent fluid. Initial fluid configuration (left) and comparison between the final configuration obtained with a fine discretization of 625 cells(a)), 16875 cells (b)),78125 cells (c)) and 135000 cells(d)). Negligible differences in the predicted level is observed, but significant differences in the local gradients of volume fraction can be noticed.
between the volume fraction distribution obtained with different mesh discretization is shown in figure 20.
530
As already shown in [10], the single-scale approach is capable of predicting the correct level. However, one
531
can observe how the interface region is significantly different from the one obtained with a finer fluid domain
532
discretization, and interests a significant part of the simulation domain. Furthermore, as shown in figure
533
22, the velocity norm shows mesh convergence only with mesh discretization finer than the one adopted for
534
the single-scale approach. This is in accordance with what observed in sections 7.1 and 7.2. In figure 23,
535
a second configuration of the benchmark is proposed in which the injection region is moved to the lighter
536
fluid, making the particles travel across the interface and then forming a porous zone at the bottom under
538
PT
CE
AC
537
ED
529
the action of gravity. The particles are injected as a pre-existing column at the beginning of the simulation, that travels across the interface and deposits at the bottom. The comparison between the averaged water
539
level obtained with different mesh discretization is shown in figure 25. One can notice how, in this case
540
characterized by more complex interface movement, a small error in the level is present for the single-scale
541
approach. This is in accordance to what has already been presented in [10]. However, this problem is solved
542
with the multiscale approach, in which grid convergence and level stability after the column injection are 32
ACCEPTED MANUSCRIPT
Hight of the column of water as a function of time
0.38
Theoretical Level a) b) c) d)
0.34 0.32 0.3
CR IP T
0.36
0
0.2
0.4
0.6
0.8
1
1.2
1.4
AN US
Averaged Water Level[m]
0.4
Time [s]
Figure 21: Multiple particles injection in quiescent fluid. Comparison between the theoretical average level and the ones obtained with the different mesh discretization of figure 20 . Good accordance with the analytic solution can be observed for all the cases.
L Norm of U as a function of time
M
8e-05
a) b) c) d)
ED
6e-05 5e-05
PT
4e-05 3e-05 2e-05
CE
L2 Norm of U [m5/s2]
7e-05
1e-05
AC
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 22: Comparison between the L2 Norm of the fluid velocity field with the domain discretization as proposed in figure 20 . Negligible grid dependence is observed for cases b), c) and d), while significant difference between these cases and the case a) is shown.
543
achieved. Once more, as shown in figure 24, the local approximation of the interface is highly affected by
544
the grid choice. In [10] and [20], the authors already pointed out how the particle size induces an important 33
ACCEPTED MANUSCRIPT
limitation for the standard DEM-VOF in the curvature calculation. In particular, in [10] they compare two
546
cases, one with particles diameter of 5 mm and one with particles diameter of 1 mm, and show how the usage
547
of small particles allows a better curvature calculation which is instead not possible for bigger particles due
548
to the nature of the DEM-VOF coupling. Similarly, in the recent work concerning the DEM-VOF approach
549
to gas-liquid-solid flows [20], the authors pointed out how the inability of the standard DEM-VOF method
550
to capture interface structures smaller than the DEM-VOF grid represents an important “drawback [...]
551
of a volume-averaging simulation”. This problem is solved by the multiscale approach proposed in this
552
contribution. As it was outlined, the key advantage of the approach is that the grid used for the interface
553
reconstruction is independent of the particle dimension. As a matter of fact, being able to numerically
554
simulate those phenomena can be valuable for many engineering processes, as the interfacial area represents
555
a key factor in many mass transport problems.
CE
PT
ED
M
AN US
CR IP T
545
556
AC
Figure 23: Injection of a column of particles inside a quiescent fluid. Relevant steps of the simulation. From left to right: initial configuration, interface perturbation due to the crossing phenomenon, particle forming a porous bed at the bottom of the domain.
In figure 26, a comparison between the L2 norms of the velocity obtained with the different mesh
557
discretization is proposed. An even more marked grid dependence than in figure 10, and 22, can be observed.
558
In particular, a peak in the norm velocity around time 0.1 s can be observed having a value almost ten times
559
higher when a single-scale approach is used. As it can be better observed in figure 27, on the long term only
560
the usage of a properly refined grid allows the L2 norm of the velocity to decay. In particular, the L2 norm 34
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 24: Injection of a column of particles inside a quiescent fluid. Comparison between the shape of the interface obtained with different grid discretization as in figure 20. Significant differences in the shape of the interface can be noted.
Average level as a function of time
M
0.42
ED
0.38 0.36
a) b) c) d)
PT
Average level [m]
0.4
0.34
CE
0.32
AC
0.3
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 25: Injection of a column of particles inside a quiescent fluid: comparison between the averaged levels obtained with the grid discretization as from figure 20 . Negligible differences between case c) and case d) can be observed, while small errors in case a) and b) can be noted.
561
shows comparable behaviors only in cases c) and d). This shows once more how grid independence of the
562
solution cannot in general be granted within a single-scale DEM-VOF coupling.
563
In figure 28, we propose a comparison between the different computational cost of the proposed simula35
ACCEPTED MANUSCRIPT
L2 Norm of the Velocity as a function of time a) b) c) d)
0.8
0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
CR IP T
0.6
1.2
1.4
AN US
L2 Norm of U [m2/s2]
1
1.6
Time[s]
Figure 26: Injection of a column of particles inside a quiescent fluid: Comparison between the L2 norm of the velocity obtained with the grid discretization as from figure 20 . Significant differences can be observed. Negligible differences between case c) and d) can be noted.
L2 Norm of the Velocity as a function of time (Close Up)
0.05
0
AC
M
PT
0.1
a) b) c) d)
ED
0.15
CE
L2 Norm of U [m2/s2]
0.2
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 27: Injection of a column of particles inside a quiescent fluid: Comparison between the L2 norm of the velocity obtained with the grid discretization as from figure 20. Close up of figure 26. Decaying behavior is observed only for case c) and d).
564
tions. This is to show how the increase in computational cost in the multiscale approach, even if extremely
565
relevant for simulations involving a single particle, becomes less important when the number of particles
566
rises. For this reason, results obtained in the cases proposed above are compared to show the dependency 36
ACCEPTED MANUSCRIPT
Single particle 1.5k particles 3k particles 500k particles
CR IP T
100
10
1
100
AN US
Clock Time/ (ClockTime SingleScale)
Normalized Clock time as a function of the number of cells
N cells / (N cells SingleScale)
Figure 28: Comparison between the increase in computational cost as a function of the refinement of the fine-scale grid, logarithmic scale. Strong dependence on the grid refinement is observed in case of Single particle, while negligible dependence is observed when a higher number of particles is adopted.
of the computational cost from the discretization used in the different cases. In order to make the various
568
scenarios comparable, we reported results for the multiscale solutions as a function of the ratio between the
569
cells adopted for the fine discretion and the ones adopted for the coarse one. For the single-scale approach,
570
this ratio will obviously be equal to unity. Similarly, for the y-axis, we defined the characteristic time as the
571
ratio between the elapsed time for the simulation and the time required by the single-scale approach. The
572
results are proposed in a logarithmic scale. One extra case is proposed that is equivalent to the one in figure
573
23, but with 500k particles of reduced diameter. For this case, only 20 fluid iterations have been solved for
574
each grid discretization and the results compared accordingly. One can see how, in cases featuring a large
575
number of particles, the computational cost is not much dependent on the grid discretization adopted. This
576
can be explained considering that in those applications, the computational burden is mainly determined by
ED
PT
CE
AC
577
M
567
the DEM part. This is normally the case for particle-laden three-phase flows, and the standard range of application for methods based on volume averaging technique. In order to better show the computational Table 3: Comparison between the execution times of the single-scale and the multiscale approaches of section 7.1
Case
Coarse Mesh
Fine Mesh
DEM δT
CFD δT
Init
1 2
None 6k
48k 500k
2x10−6 [s] 2x10−4 [s] 1.8[s] 2x10−6 [s] 2x10−4 [s] 8.4[s] 37
DEM Loop
CFD Loop
Total
6346[s] (98%) 6226[s] (73%)
130[s] (2%) 2224[s] (26%)
6476[s] 8450[s]
ACCEPTED MANUSCRIPT
intensity of the DEM and the CFD part for a standard case, we reported in table 3 the detailed data for
579
the simulation of section 7.1. The two cases confronted are a classical single-scale DEM-VOF simulation,
580
where a mesh of 48k cells and no inter-mesh interpolation are used, and the multiscale DEM-VOF operated
581
on a 6k-cells coarse grid and 500k-cells fine grid. One can observe how the increase in overall computational
582
time is rather small and generally in accordance with what is shown in figure 28. It is interesting to observe
583
how the time spent in the DEM computation, that includes the projection of the Lagrangian field on the
584
grid, is smaller in the multiscale approach if compared to the single-scale one. This shows what was antic-
585
ipated in section7.1 and figure 14 . The usage of a coarser grid for the projection of the Lagrangian fields
586
allows reducing the computational burden of this operation. On the other hand, the cost of the CFD part,
587
that includes the mapping between meshes, increases with the multiscale approach. Nevertheless, as this
588
represents a fraction of the total cost, its influence on the total computational burden of the simulation is
589
relatively small.
590
8. Conclusions
AN US
CR IP T
578
A novel multiscale approach for three-phase flows has been presented. The method is based on the
592
volume averaging technique and aims to overcome some of its main drawbacks that limit its extension to
593
complex fluid flows. The numerical scheme is capable of providing solutions for the three-phase problem
594
using an arbitrary discretization for the fluid equations. At the same time, it keeps the computational
595
advantages from coupling the particle-fluid equations with local exchanges of momentum based on drag
596
force calculation. Compared to existing methods, this multiscale approach succeeded in taking into account
597
effects of complex interface dynamics in the description of particle-laden free surface flows. Therefore, it
598
enlarges the computational window to a wide area of real case phenomena. Furthermore, the ability of the
599
multiscale DEM-VOF method to reconstruct finer structures of the flow, and in particular of the gas-liquid
600
interface, makes it suitable to study transport phenomena more accurately. The method has been shown
601
to produce grid convergent solutions. This allows providing an important insight on the correctness of
603
604
ED
PT
CE
AC
602
M
591
the solution. It has been shown how the multiscale DEM-VOF approach is able to agree with benchmark problems presented in the literature, and in particular is able to reproduce experiments with a reduced error if compared to the standard DEM-VOF. Furthermore, the method allows adopting a wider range of volume
605
replacement algorithms, opening the possibility to capture phenomena of bubble and droplet formation, in
606
accordance with the classical VOF theory. Finally, the increase in computational cost of the new method
607
with respect to the classical single-scale approach has been shown to be negligible in cases involving by an 38
ACCEPTED MANUSCRIPT
608
elevated number of particles. The main limitation of the method consists in the higher complexity of the algorithm, that requires a
610
mapping operation between the different grids. This may make its extension to cases involving reacting
611
flows, and heat and mass transfer between particles and fluid, more complex, as the requirements for the
612
mapping operations can become more delicate. Additionally, an efficient parallelization of the algorithm will
613
require a dedicated study due to its enhanced complexity. As a future investigation, the influence of different
614
interpolation strategies on the convergence properties of the algorithm can be studied. Furthermore, the
615
extension of the approach to turbulent configurations is of high interest. In particular, the proposed dual-grid
616
algorithm here proposed brings significant advantages, as the flexibility in the fluid domain discretization
617
can play a fundamental role in the reconstruction of boundary layers and turbulent structures. However,
618
the introduction of sub-grid models in the current algorithm is a topic that must be carefully addressed,
619
with particular reference to the influence of those models on the interface dynamics, and the treatment of
620
sub-grid terms involving the surface tension.
AN US
CR IP T
609
The extension of the multiscale approach presented in this paper to CFD-DEM coupling different from
622
the DEM-VOF is another potential outcome of the present work. It is our opinion that several other
623
couplings between DEM and complex flow can benefit from it.
624
Acknowledgments
ED
625
M
621
The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg [24].
627
fruitful discussions. Finally, we would like to thank Florian Hoffmann, Jack S. Hale, and Amir Houshang
628
Mahmoudi for their invaluable help and suggestions in the last revision of the manuscript.
629
References
633 634 635 636 637 638 639 640 641 642 643 644
CE
632
[1] Afkhami, M., Hassanpour, A., Fairweather, M., Njobuenwu, D., 2015. Fully coupled LES-DEM of particle interaction and agglomeration in a turbulent channel flow. Computers & Chemical Engineering 78, 24–38. URL http://www.sciencedirect.com/science/article/pii/S0098135415001040 [2] Afkhami, S., Zaleski, S., Bussmann, M., 2009. A mesh-dependent model for applying dynamic contact angles to VOF simulations. J. Comput. Physics 228 (15), 5370–5389. URL http://dblp.uni-trier.de/db/journals/jcphy/jcphy228.html#AfkhamiZB09 [3] Anderson, T. B., Jackson, R., 1967. Fluid Mechanical Description of Fluidized Beds. Equations of Motion. Industrial & Engineering Chemistry Fundamentals 6 (4), 527–539. URL http://dx.doi.org/10.1021/i160024a007 [4] Berberovi´ c, E., van Hinsberg, N. P., Jakirli´ c, S., Roisman, I. V., Tropea, C., Mar 2009. Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Phys. Rev. E 79, 036306. URL http://link.aps.org/doi/10.1103/PhysRevE.79.036306 [5] Du, W., Bao, X., Xu, J., Wei, W., 2006. Computational fluid dynamics (CFD) modeling of spouted bed: Assessment of drag coefficient correlations. Chemical Engineering Science 61 (5), 1401–1420. URL http://www.sciencedirect.com/science/article/pii/S0009250905006603
AC
630 631
We wish to thank the anonymous reviewers and the editor for the careful work and
PT
626
39
ACCEPTED MANUSCRIPT
653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
CR IP T
652
AN US
651
M
650
ED
649
PT
648
CE
647
[6] Herrmann, H. J., 1995. Physics of granular media. Chaos, Solitons & Fractals 6, 203–212. URL http://www.sciencedirect.com/science/article/B6TJ4-49WMG9H-X/1/f106bb7c6924a76a8a0c156a99a07718 [7] Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39 (1), 201–225. URL http://www.sciencedirect.com/science/article/pii/0021999181901455 [8] H.Jasak, 1996. Error analysis and estimation for the finite volume method with applications to fluid flows. Ph.D. thesis. [9] Jasak, H., Weller, H. G., Gosman, A. D., Sep. 1999. High resolution NVD differencing scheme for arbitrarily unstructured meshes. International Journal for Numerical Methods in Fluids 31, 431–449. [10] Jing, L., Kwok, C. Y., Leung, Y. F., Sobral, Y. D., 2016. Extended CFD–DEM for free-surface flow with multi-size granules. International Journal for Numerical and Analytical Methods in Geomechanics 40 (1), 62–79, nAG-14-0182.R1. URL http://dx.doi.org/10.1002/nag.2387 [11] Li, Y., Zhang, J., Fan, L.-S., 1999. Numerical simulation of gas–liquid–solid fluidization systems using a combined CFDVOF-DPM method: bubble wake behavior. Chemical Engineering Science 54 (21), 5101–5107. URL http://www.sciencedirect.com/science/article/pii/S0009250999002638 [12] Ling, Y., Zaleski, S., Scardovelli, R., 2015. Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model. International Journal of Multiphase Flow 76, 122–143. URL http://www.sciencedirect.com/science/article/pii/S0301932215001524 [13] Mahmoudi, A. H., Hoffmann, F., Markovic, M., Peters, B., Brem, G., 2016. Numerical modeling of self-heating and selfignition in a packed-bed of biomass using {XDEM}. Combustion and Flame 163, 358–369. URL http://www.sciencedirect.com/science/article/pii/S0010218015003582 [14] Mahmoudi, A. H., Hoffmann, F., Peters, B., 2014. Detailed numerical modeling of pyrolysis in a heterogeneous packed bed using {XDEM}. Journal of Analytical and Applied Pyrolysis 106, 9–20. URL http://www.sciencedirect.com/science/article/pii/S0165237013002660 [15] Mahmoudi, A. H., Hoffmann, F., Peters, B., 2016. Semi-resolved modeling of heat-up, drying and pyrolysis of biomass solid particles as a new feature in {XDEM}. Applied Thermal Engineering 93, 1091–1104. URL http://www.sciencedirect.com/science/article/pii/S1359431115010807 [16] Peters, Bernhard, Pozzetti, Gabriele, 2017. Flow characteristics of metallic powder grains for additive manufacturing. EPJ Web Conf. 140, 13001. URL https://doi.org/10.1051/epjconf/201714013001 [17] Pozzetti, G., Peters, B., 2017. On the choice of a phase interchange strategy for a multiscale DEM-VOF Method. AIP Conference Proceedings 1863. URL http://orbilu.uni.lu/handle/10993/28863 [18] Rui, S., Heng, X., 2015. Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD–DEM. International Journal of Multiphase Flow 72, 233–247. [19] Schiller, L., Naumann, Z., 1935. A Drag Coefficient Corre-lation,. VDI Zeitung Vol. 77,, 318–320. [20] Sun, X., Sakai, M., 2015. Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method. Chemical Engineering Science 134, 531–548. URL http://www.sciencedirect.com/science/article/pii/S0009250915003991 [21] Sun, X., Sakai, M., Yamada, Y., 2013. Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method. Journal of Computational Physics 248, 147–176. URL http://www.sciencedirect.com/science/article/pii/S0021999113002684 [22] Tryggvason, G., Scardovelli, R., Zaleski, S., January 2011. Direct Numerical Simulations of Gas-Liquid Multiphase Flow. Cambridge University Press. [23] Uhlmann, M., 2005. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics 209 (2), 448–476. URL http://www.sciencedirect.com/science/article/pii/S0021999105001385 [24] Varrette, S., Bouvry, P., Cartiaux, H., Georgatos, F., July 2014. Management of an Academic HPC Cluster: The UL Experience. In: Proc. of the 2014 Intl. Conf. on High Performance Computing & Simulation (HPCS 2014). IEEE, Bologna, Italy, pp. 959–967. [25] Weller, H. G., Tabor, G., Jasak, H., Fureby, C., 1998. A tensorial approach to computational continuum mechanics using object-oriented techniques. Computers in Physics 12 (6), 620–631. URL http://scitation.aip.org/content/aip/journal/cip/12/6/10.1063/1.168744 [26] Wu, T.-R., Chu, C.-R., Huang, C.-J., Wang, C.-Y., Chien, S.-Y., Chen, M.-Z., 2014. A two-way coupled simulation of moving solids in free-surface flows. Computers & Fluids 100, 347–355. URL http://www.sciencedirect.com/science/article/pii/S0045793014002023 [27] Xiao, H., Sun, J., 2010. Algorithms in a Robust Hybrid CFD-DEM Solver for Particle-Laden Flows. Communication Computer Physics. [28] Zhu, H., Zhou, Z., Yang, R., Yu, A., 2007. Discrete particle simulation of particulate systems: Theoretical developments. Chemical Engineering Science 62 (13), 3378–3396, frontier of Chemical Engineering - Multi-scale Bridge between Reductionism and Holism. URL http://www.sciencedirect.com/science/article/pii/S000925090700262X
AC
645 646
40
A Multiscale DEM-VOF method for the simulation of three-phase flows Gabriele Pozzetti, Bernhard Peters PII: DOI: Reference:
S0301-9322(16)30732-7 10.1016/j.ijmultiphaseflow.2017.10.008 IJMF 2664
To appear in:
International Journal of Multiphase Flow
Received date: Revised date: Accepted date:
11 December 2016 9 October 2017 9 October 2017
Please cite this article as: Gabriele Pozzetti, Bernhard Peters, A Multiscale DEM-VOF method for the simulation of three-phase flows, International Journal of Multiphase Flow (2017), doi: 10.1016/j.ijmultiphaseflow.2017.10.008
This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.
ACCEPTED MANUSCRIPT
A Multiscale DEM-VOF method for the simulation of three-phase flows Gabriele Pozzettia,∗, Bernhard Petersa a
Campus Kirchberg, Universit du Luxembourg 6, rue Richard Coudenhove-Kalergi L-1359 Luxembourg
CR IP T
Abstract
A novel multiscale approach for three-phase flows is presented. The goal of the proposed method is to solve arbitrary three-phase flow configurations in a computationally efficient way, and in particular taking into account the effects of different length scales while sharply reducing the computational burden. This is particularly important in chemical, environmental, and process engineering, where large-scale quantities are
AN US
normally of interest, but small-scale dynamics cannot be neglected. The method is based on the definition of two different length scales: the bulk scale, and the fluid fine scale. A dual-grid approach is adopted in order to resolve the bulk scale with information from the fluid fine scale. It is shown that the described method succeeds in delivering more accuracy than a standard approach based on the volume averaging technique, still, it remains suitable for the solution of real interest problems. The method is shown to successfully satisfy
M
experimental results presented in the literature. Examples of three-phase flows simulations are provided to show how the proposed numerical approach can describe the physics of particle-laden, free surface flows with
ED
competitive computational cost. It is shown how the proposed method can naturally extend the DEM-VOF method to the presence of complex interface dynamics.
PT
Keywords: Multiphase flows, CFD-DEM coupling, Grid convergence, Multiscale simulation, DEM-VOF
CE
method
Highlights
AC
• A multiscale DEM-VOF coupling is proposed to simulate three-phase flows • It produces grid-convergent solutions for the fluid phase • It can reproduce droplet formation and interface perturbations from a single particle • It is validated against three-phase flows presented in the literature
∗ Principal
corresponding author. Email: [email protected] Email address: [email protected] (Gabriele Pozzetti)
Preprint submitted to Elsevier
October 12, 2017
ACCEPTED MANUSCRIPT
1
1. Introduction Gas-liquid-particles flows are extremely common in industrial and environmental applications, and repre-
3
sent a standard configuration for chemical and process engineering. The engineering community is currently
4
demanding numerical tools able to predict the behavior of such complex flows at reasonable computational
5
cost. In particular, a direct numerical approach to all the phenomena that arise in these flows is unbearable
6
for nowadays computational capacities. Therefore, a considerable modeling effort is required to provide
7
reliable solutions with competitive cost.
CR IP T
2
The dispersed solid phase often plays a key role in the physics of the problem, being responsible for
9
the main chemical and dynamical interactions of the system. A method capable of describing the discrete
10
dynamics of the granular phase could, therefore, be extremely interesting [28]. For this reason, CFD-DEM
11
approaches have become very popular as they are able to combine the accuracy of the Discrete Element
12
Method with the highly developed CFD techniques for the study of complex flows.
AN US
8
Recently, significant efforts have been made for coupling DEM solvers with multiphase fluid solvers using
14
the Volume of Fluid (VOF) technique. This approach has been shown promising to provide reliable numerical
15
solutions of the three-phase flow dynamic and its influence on relevant chemical processes [20, 26, 11, 10, 17].
16
The VOF method is a relatively old approach for solving free surface flows, that after several modifications,
17
is still used today [22]. This gave rise to the development of multiple VOF schemes that try to match
18
a compromise between computational efficiency, and accuracy [22]. Generally, a VOF scheme aims to
19
reconstruct the position of the interface location. The interface reconstruction is particularly critical as it
20
is used to calculate the local curvature which is proportional to the surface tension force that will appear
21
as a contribution in the momentum equation. Furthermore, within a VOF scheme, a volume fraction field
22
is transported by the velocity field in order to define the new interface location.
CE
PT
ED
M
13
Two approaches, commonly referred to as algebraic and geometrical schemes, are widely used today.
24
Geometrical VOF techniques reconstruct the interface using volume fraction at each cell and its neighbors,
25
so to approximate interface information such as area and intersections with cell faces. Many variations of
26
27
AC
23
this strategy have been presented in the literature, and an extended description of various VOF schemes can be found in [22]. Algebraic advection methods perform a discretization of the advection equation by
28
reconstructing volume fluxes from the solution of the velocity field. This is used to compute the mixture
29
volume fraction and its flux over cells. A commonly used algebraic advection scheme is implemented in
30
R [25]. Various modification to the method have been proposed, the algebraic scheme adopted OpenFOAM
31
in this contribution is more extensively described in section 4. 2
ACCEPTED MANUSCRIPT
Depending on the technique chosen to couple the Eulerian and the Lagrangian domains, CFD-DEM
33
couplings can lead to extremely high computational cost [3, 23, 15]. In order to address industrially significant
34
problems, a widely accepted solution is the volume averaging approach[3]. It consists of the projection
35
of Lagrangian fields into Eulerian ones through the usage of volume-averaged variables[3]. CFD-DEM
36
approaches based on the volume averaging technique aim to perform a zero-dimensional coupling between
37
CFD and DEM solution, using cell values or locally interpolated values for the evaluation of the fluid-particle
38
interactions. This significantly reduces the computational burden and represents one of the main advantages
39
of the method. The VOF method, and in general interface-capturing methods, are particularly interesting
40
to be used within these couplings, as they allow solving the well-known single-phase equations within each
41
single fluid phase, smearing the interface region between two phases over a few cells. This can represent a key
42
advantage for a volume-averaged CFD-DEM coupling, as it allows adopting the highly developed empirical
43
models for the coupling between particles and single-phase fluid within the major part of the computational
44
domain.
AN US
CR IP T
32
Due to their attempt to perform inexpensive couplings, volume averaging approaches naturally operate
46
on a length scale at which the exchange of momentum between particles and fluid can be considered as
47
zero-dimensional. We will refer to this scale as the bulk or coupling scale of the simulation. We will instead
48
refer to fluid fine scale as the smallest length scale that needs to be resolved to produce grid independent
49
results for the fluid.
ED
M
45
In standard CFD-DEM couplings based on the volume averaging technique, the bulk scale is automat-
51
ically chosen as the one at which the fluid equations will be written and solved. Therefore, this technique
52
naturally applies to cases in which the fluid fine scale is significantly larger than the particles characteristic
53
dimensions. For this reason, the volume averaging approach was originally proposed for couplings in which
54
the CFD cells can be kept significantly larger than the largest particle diameter. However, it nowadays
55
became a standard to apply this coupling technique with cell-particle diameter ratios close to the unity, by
56
using complex averaging procedures [18].
58
59
CE
AC
57
PT
50
In a VOF based simulation, the fluid fine scale is associated with the interface dynamic, that must be
correctly captured to achieve grid convergence. Generally, the ability of a VOF scheme to reconstruct the interface is bounded to the grid spacing. The particular requirement will depend on the specific scheme
60
adopted, but in all the cases a correct reconstruction of the interface is mandatory as it will affect the cal-
61
culation of the surface tension terms, and therefore the determination of the velocity fields, and finally the
62
advection of the volume fraction. In addition to that, in algebraic VOF schemes, the usage of a sufficiently 3
ACCEPTED MANUSCRIPT
fine discretization is mandatory in order to deal with the problem of spurious currents. This can lead to an
64
erroneous calculation of the fluxes, and therefore a poor conservation of the overall scheme [2, 4]. Further-
65
more, correct grid resolution is particularly important in case of high density ratios. In these configurations,
66
small errors in the volume fraction can easily cause major problems for momentum conservation, leading in
67
the end to a poor accuracy of the solution [4]. For this reason, to extend DEM-VOF couplings to complex
68
interface scenarios, it is important to develop algorithms that can ensure the correct grid resolution in order
69
for the VOF scheme to properly converge.
CR IP T
63
In purely CFD cases, a higher grid resolution is directly linked to an increase in computational cost, that
71
represents a major challenge for the numerical solution of these configurations. In a DEM-VOF coupling
72
based on the volume averaging technique, the scale separation between the simulation domain (or integral
73
scale), the bulk scale for CFD-DEM coupling, and the fluid fine scale, can represent an even bigger problem.
74
This is due to the fact that the bulk scale is the smallest one that a single-scale volume averaging coupling
75
can resolve, but the solution of the fluid fine scale is needed to provide unique solutions.
76
77
AN US
70
The main objective of this contribution is to propose a multiscale method that overcomes this problem. The method we propose uses a dual-grid approach where:
• A coarse grid is used to perform the CFD-DEM coupling at the bulk scale
79
• A fine grid is used to resolve the equations governing the fluid phases at the fluid fine scale.
ED
M
78
Those two scales can have significantly different characteristic lengths, and a correct modeling for both of
81
them is mandatory in order to reach a compromise between accuracy and computational feasibility. This
82
dual-grid multiscale approach aims to capture the most significant phenomena of the two different scales,
83
while limiting as much as possible the computational burden.
PT
80
This method is meant to bring the most advantages when the separation between those two scales is more
85
marked. That happens when the smallest structures of the fluid flow that must be resolved to obtain grid
86
independent results, are significantly smaller than the ones determining the bulk fluid-particle momentum
88
AC
87
CE
84
exchange. In section 2, we refer to the Stokes number as a parameter to identify this scale separability. Compared to the existing methods, our novel multiscale approach is thought to have several advantages.
89
First, it allows achieving grid-convergence of macroscopic flow quantities in cases in which the single-scale
90
DEM-VOF method cannot. Second, it allows more accurate predictions of the interface dynamics that are
91
fundamental to capture the most significant fluid dynamic instabilities. This leads to a highly versatile
92
numerical method that can be applied to several flow configurations. The method can also be naturally 4
ACCEPTED MANUSCRIPT
extended to other CFD-DEM couplings once the length scales are properly individuated and defined. This
94
paper is organized as follows: First, the novel dual-grid coupling procedure is proposed in its more general
95
configuration. Then, the governing equations for both the fluid and the granular phases are presented.
96
The model for the coupling terms between continuum and discrete phases is subsequently proposed with
97
particular reference to the phase interchange problem, and to the interpolation strategies adopted. Some
98
reference cases are finally discussed. To the best of our knowledge the multiscale DEM-VOF method is the
99
first that allows capturing the effect of a single particle on a liquid interface, still keeping the computational
100
cost of a volume-averaged simulation.
101
2. Method description
CR IP T
93
A fundamental nondimensional number in the study of particle-laden flow is the Stokes number. It
103
represents the ratio between the particles and the fluids characteristic time scales and, when the particle
104
Reynolds number is less than unity, can be expressed as: St =
AN US
102
ρp d2p |uf | , dI 18µf
(1)
with ρp the density of the particles, dp their diameter, uf the fluid velocity, dI the hydraulic diameter,
106
µf the viscosity of the fluid. In this paper, we refer to the case of high Stokes number flows, for which a
107
single particle has relatively large mass and diameter and therefore its single action plays an important role
108
in the process dynamic. In this kind of flow, a particle’s motion is mainly influenced by large structures
109
of the fluid due to its inertia. Therefore, the coupling between the CFD and DEM can be operated on a
110
relatively coarse scale, where only the bulk structures of the flow are resolved. Nevertheless, when the fluid
111
phases are characterized by high-density ratios and complex interface dynamics, their bulk behavior cannot
112
be univocally determined without resolving small-scale dynamics. In particular, to correctly predict the
113
phase distribution, the interface must be correctly advected, and this requires detailed information on the
114
velocity field. Furthermore, the very solution of the velocity field is influenced by the local behavior of the
116
117
ED
PT
CE
AC
115
M
105
interface position and curvature. For these reasons, in order to obtain grid independent solutions, one often has to resolve a fluid scale smaller than the one used for the CFD-DEM coupling. In this contribution, we propose a dual-grid algorithm for the CFD-DEM coupling that overcomes this issue, and allows obtaining
118
grid convergent solutions for the fluid while keeping the computational advantages of a volume-averaged
119
simulation.
120
The dual-grid algorithm proposed to achieve this is described in figure 2, and structured as follows. A 5
ACCEPTED MANUSCRIPT
Table 1: List of variables
Symbol Variable
Symbol Variable
uf
fluid velocity
α
µf ρf dI Ffpi
uc Γ ψ C
TΓ
fluid viscosity fluid density hydraulic diameter volumetric source of momentum due to interaction with particles surface tension force
Fb p
body force pressure
σf τf
npc
number of particles in a cell
g
CR IP T
h
Dirac’s delta function describing the interface fluid stress tensor deviatoric part of the stress tensor gravitational acceleration modulus hight
kf
free surface compression factor
Φm Φjm
mass flux partial mass flux (relative to a sub-cycling step) uf F · AF
δTtot
face Area face Area normal vector (pointing outward the owner cell) compression velocity face value = uc · AF simulation time step
δTi
sub-cycle time step
aN
H
non diagonal contribution of the velocity matrix for momentum predictor
ρp dp
particle density particle diameter
AN US
ucF
η
M
AF AF
aP
ED
χ
P porosity field: i Vpi /Vcell with i = 1....npc phase indicator
PT
δΓ
volume fraction of the liquid phase compression velocity Surface Tension particle volume concentration Curvature
m
particle mass
Mext
I Fc Fij
particle moment of inertia contact force acting on a particle contact force acting on the particle i due to the collision with particle j force rising from the particlefluid interaction particle-fluid interface normal vector
up Ap φ
gravitational force torque acting on a particle due to collision events external torque acting on a particle particle velocity particle surface area particle orientation
ω
particle angular velocity
particle based Reynolds number interface local velocity interface normal vector
ρh β CD
density of the heavy phase drag factor drag coefficient
AC
CE
Fg Mc
the matrix coefficient related to the centroid (velocity equation) the coefficients of neighbor centroids (velocity equation)
Ff n
Re uint nint
6
ACCEPTED MANUSCRIPT
Bulk Scale
Fluid 2
CR IP T
Fluid 1
Fluid Fine Scale
AN US
Fluid 2
Fluid 1
M
Figure 1: Different length scales in high-Stokes three-phase flows: bulk (coarse) scale and fluid fine scale.
first coarse mesh is used to construct averaged fields from the granular-phase distribution, according to the
122
standard volume averaging technique. This mesh is used to solve the fluid-particle coupling consisting in
123
the estimation of the drag force, the calculation of the porosity field, and the definition of the particle-
124
related source of momentum in the fluid. Different strategies can be chosen to perform the averaging of
125
particle-related information to the coarse mesh. This operation is also referred to as backward interpolation.
126
Similarly, multiple choices are possible to interpolate the flow fields at the location of the particles. In
127
this contribution, the procedure proposed in [27] is adopted to perform the backward interpolation, while
128
R functions are adopted for the forward interpolation. Also, the semi-implicit the standard OpenFOAM
129
treatment algorithm for the calculation of the fluid-particle interaction terms of [27] is used.
131
PT
CE
AC
130
ED
121
After the fluid-particle interaction problem is solved, we project the updated variables in a second mesh
where the fluid fine scale is resolved. We will refer to this second mesh as the supporting domain of our
132
approach. This procedure is described in section 6. The fluid fine-scale problem is in this way initialized
133
with particle-related information provided by the coarse-scale solution. In the obtained set of equations,
134
the variables derived from the Lagrange fields projection are defined on an averaged scale. The usage of a
135
supporting domain for the fluid equations solution forces us to adopt those averaged fields, but enables us 7
ACCEPTED MANUSCRIPT
-Mapping Lagrangian fields to Eulerian. -Solving fluidparticles interactions
CR IP T
Fluid Solution
Particle related fields
Coarse Mesh
- Solving Fluid equations
AN US
Fine Mesh
Figure 2: Diagram of the solution procedure of the two different length scales for the simulation. The two boxes represent the different models adopted, while the arrows show schematically the communication between the scales. A coarse mesh (top) is used to map the Lagrangian field of the DEM into an Eulerian reference and to solve the fluid-particle interaction. Particle-related fields are mapped to the supporting domain (bottom) where a finer mesh is used to solve the fluid equations.
to sharply reduce the computational burden. The discretization of the supporting domain can be done in
137
accordance with the fluid system requirements, and in general can be totally independent of the particles’
138
characteristics. The usage of sub-grid modeling at the fluid fine scale is also in theory possible, but goes
139
beyond the purpose of this paper, in which we will assume the effect of scales smaller than the fluid fine one
140
as negligible.
PT
ED
M
136
After the set of equations of the supporting domain has been solved, the fluid variables are mapped to the
142
coarse domain where averaged variables are reconstructed from the up to date data. It must be noted that
143
using a finer discretization for the solution of the fluid equations does obviously increase the computational
144
burden with respect to a classical DEM-VOF simulation. However, CFD-DEM couplings based on the
145
volume averaging technique are normally applied in simulation featuring O(103 ) − O(106 ) particles. In this
147
148
AC
146
CE
141
range, the DEM often represents the more expensive part of the simulation. Therefore, as better shown in section 7.3, the additional cost related to the usage of a finer grid can, in many cases, be negligible. 3. Granular phase governing equations
149
The Discrete Element Method (DEM) is a well-established approach for granular materials [6]. It has
150
been proven to be capable of successfully addressing the physics of granular flows in an efficient way. Its 8
ACCEPTED MANUSCRIPT
151
advantages gave rise to multiple attempts to extend DEM solvers by coupling them with flow fields resolved
152
by CFD [14, 10, 13, 16]. In this work, the DEM is used to evolve a set of particles describing physical entities
153
in the presence of a multiphase fluid. Defining with xi the positions, mi the masses, φi the orientation, and
154
Ii the tensorial moment of inertia of the particles one can write d2 xi dt2 d2 Ii 2 φi dt
155
= Mc + Mext ,
where Fc is the collision force Fc =
X
Fij (xj , up j , φj , ωj ),
i6=j
(2)
CR IP T
= Fc + Ff + Fg ,
mi
P
(3)
(4)
with up the particle velocity, and ω the angular velocity. In particular,
157
particles different from particle i, while Fij is the interaction force between particle i and particle j. At the
158
same time, the term Mc indicates the torque acting on the particle due to collisions Mc =
X
AN US
156
i6=j
represents the sum over all
Mij (xj , up j , φj , ωj ),
i6=j
(5)
with Mij the torque acted from particle j to particle i. A detailed description of all the models present in
160
the literature for the equations 4 and 5 is presented in [28], and citations within. In equation 2 the term Fg
161
represents the gravitational force and Ff takes into account the force rising from the interaction with the
162
fluid. An analytic expression for Ff can be obtained as
PT
ED
M
159
Ff =
Z
S
σf · n dγ,
(6)
with σf the fluid stress tensor, n the normal vector to the particle-fluid interface, S the particle surface.
164
The specific treatment of this term will be addressed in section 5.1.
165
4. Governing equations of the fluid scale
167
168
AC
166
CE
163
In this contribution, we refer to an unsteady incompressible, immiscible, two-fluid system, with forcing
terms arising from the third phase in accordance with section 5.1. The equations governing this system are the Navier-Stokes equations in the form ∂ρf uf + ∇ · (ρf uf uf ) = −∇p + ∇ · µf ∇uf + ∇T uf + TΓ + Fb + Ffpi , ∂t 9
(7)
ACCEPTED MANUSCRIPT
169
with TΓ = ΓCnδΓ the surface tension force, Fb a generic body force, Ffpi the fluid-particle interaction force
170
as shown in section 5.1. Within each phase, the incompressibility constraint must be fulfilled in the form ∇ · uf = 0.
(8)
In diluted flow configurations, the presence of the discrete phase is commonly taken into account by the
172
term Ffpi [1], while in configurations characterized by dense granular phase (like for instance, packed beds)
173
porous equations are normally preferred. If one defines the porosity as = Vsolid /Vcell , the equations take
174
the form
CR IP T
171
∂ρf uf + ∇ · (ρf uf uf ) = −∇p + ∇ · µf ∇uf + ∇T uf + TΓ + Fb + Ffpi , ∂t
(9)
where the terms Fb , and Ffpi , must be changed accordingly [5]. Similarly, the incompressibility constraint
176
imposes that, within a control volume, changes in the void fraction are taken into account leading to
AN US
175
∇ · uf =
ρf (x) = ρ1 α(x) + ρ2 (1 − α(x)),
(11)
µf (x) = µ1 α(x) + µ2 (1 − α(x)),
(12)
ED
178
with α the volume fraction defined by
=
χ =
Z 1 χ(x)dx, V V 1 if first fluid, 0
(13)
(14)
if second fluid.
The volume fraction is considered a scalar transported by the fluid flow for which
AC
180
CE
PT
α
179
(10)
Local density and viscosity are dependent on the fluid phase and can be written in the form
M
177
∂ . ∂τ
∂α + ∇ · (αuf ) + ∇ · (α(1 − α)uc ) = 0, ∂t
(15)
must hold, with uc the relative velocity between the two-phases referred to as compression velocity. The
181
third term is introduced in order to avoid an excessive numerical dissipation.
182
4.1. Finite-volume discretization of the fluid equations
183
R [25], an extensively-used openThe computations of the flow fields are performed using OpenFOAM
184
R structure, the control source library for field operation and manipulation. According to the OpenFOAM
10
ACCEPTED MANUSCRIPT
volumes can be of general polyhedral shape, all depended variables share the same control volumes, and
186
equations are treated in a segregated way. We adopt, and here briefly recall, the VOF scheme described in
187
[4]. In particular, the changes proposed by [10] are used for the porous cases. The finite volume method
188
(FVM) is used to discretize the conservation equations. For the source terms as well as the transient term
189
the midpoint rule is adopted and an integration over the cell volume is performed. Convective and diffusive
190
terms are treated converting the volume integrals into surface integrals (Gauss theorem). The integration
191
is then performed over the surface values. Gradients are finally estimated with a linear face interpolation,
192
with a correction term that takes into account the non-orthogonal contributions. For the convective term,
193
a normalized variable diagram (NVD) based face interpolation is chosen, with a high-resolution differencing
194
scheme proposed in [9], in order to ensure boundedness in the discretization. An Euler implicit scheme is
195
used for the temporal discretization [8]. As proposed in [8], the method expresses the face values required
196
by the discretization of diffusive and convective terms with the new time-level cell values. This results in
197
a first-order accurate method that guarantees boundedness of the solution. In particular, the fluid-particle
198
interaction term is treated in a semi-implicit way according to [27]. As proposed in [4] the term uc of equation
199
15 is evaluated as a conservative volume flux derived from the pressure-velocity coupling. In particular, its
200
expression is
M
AN US
CR IP T
185
204
205
ED
PT
203
normal face flux calculated from the gradient of the volume fraction as ξF =
∇α · AF , ||∇α| + ∆|
(17)
with ∆ a stabilization factor introduced for non uniform grids, evaluated on the mesh characteristic by
CE
202
(16)
with η the face volume flux, kf the free surface compression factor here set to 1. Finally, ξF is the face
AC
201
η , max η , uc,F = ξF min kf |AF | |AF |
∆ = P
B
Nc
Vi
Nc
13 ,
(18)
with B a small parameter (here set to 10−7 ) and N the number of computational cells. A temporal sub-
cycling strategy for the advection of the phase is adopted. This is motivated by the fact that for VOF-based
206
methods the stability of the solution procedure is very sensitive to the phase advection problem. In order to
207
address this problem, in addition to the usage of bounded discretization schemes for the divergence terms,
208
the solution of the volume fraction equation is divided into multiple sub-cycles. The time step of the single 11
ACCEPTED MANUSCRIPT
209
sub-cycle is obtained as δTi =
δTtot , Nsc
(19)
with Nsc the number of sub-cycles, δTtot the simulation time step. Within a sub-cycle, the volume fraction
211
α is calculated, and from this the total mass flux Φm is obtained as ΦG m = ρf uf · AF =
Nsc X δTi j Φ , δ tot m j=1 T
212
with Φjm the partial mass flux, Nsc the number of sub-cycles.
213
4.2. PIMPLE loop for the porous equations
CR IP T
210
(20)
R , applied on the We here briefly describe the standard PIMPLE loop as implemented in OpenFOAM
215
porous equation 9. First, the porosity field is mapped from the coarse mesh into the centroids points of the
216
fine mesh, and the porosity values on the cell faces are reconstructed. Then, the matrix for the velocity
217
equation is assembled, and the momentum predictor is solved in the form
AN US
214
au P ufP = H − ∇p − gh∇ρf + ΓC∇α + Ffpi , where au P is the matrix coefficient related to the centroid point, and
M
218
(21)
ED
H = ∇ · µf ∇uf + ∇T uf
−
X
au N ufN ,
(22)
N
with aN the coefficients for the neighbor centroids. The velocity equation is relaxed to ensure stability and
220
solved accordingly. The pressure-correction is then performed as follows. From equation 21, a predicted
221
velocity is obtained as
PT
219
CE
−1 ufP = [au [H − ∇p − gh∇ρf + ΓC∇α + Ffpi ] , P]
(23)
and the face values of the velocity ufF is estimated through face interpolation. The pressure equation is
223
then solved as
AC
222
∂ −1 ∇ · 2 [au ∇p = ∇ · uf + , P] ∂t
(24)
224
that is the discretized version of equation 10, where the time derivative of the porosity fields takes into
225
account the changes of the cell volume due to particle presence. The velocity field is then corrected with
226
the resolved pressure field, so to obtain conservative fluxes, and therefore a divergence-free velocity field.
12
ACCEPTED MANUSCRIPT
227
5. Fluid particle system
228
5.1. Momentum coupling
230
As a standard procedure in fluid dynamics, the fluid stress tensor in equation 6 can be divided into pressure and residual strain contribution, leading to Ff = −
231
Z
S
p I · n dγ +
Z
S
τf · n dγ = Flif t + Fdrag ,
CR IP T
229
with τfij = σfij + pδij ,
(25)
(26)
and p = −1/3 tr(σf ) the fluid pressure, δij the Kronecker delta, I the identity matrix. If the fluid-particle
233
interaction is aimed to be solved on an averaged scale, it is possible to consider the pressure gradient of a
234
region of fluid surrounding the particle as constant. In this scenario, it is convenient to express the first
235
integral as Flif t = −
Z
S
AN US
232
p I · n dγ =
Z
Vp
∇p dV,
(27)
where the Green theorem has been used, Vp is the volume of the particle. Considering the pressure gradient
237
as constant within the integration volume, one obtains
M
236
ED
Flif t = −∇p Vp ,
(28)
with p pressure of the liquid. In a gas-liquid-particles system, the pressure gradient between the liquid and
239
the gas region is normally not negligible due to the difference in density typical of a gas-liquid configuration.
240
This effect must, therefore, be taken into account.
243
244
CE
242
A well-stated method to study the transport of particulate matter into a fluid flow is via the usage of a drag law that takes into account the average effect of the fluid stress on the solid
AC
241
PT
238
Fdrag =
Z
S
τf · n dγ .
(29)
Many authors have focused on the description of the drag force [28, 20] that is normally provided in the form of a semi-empirical law Fdrag = β(uf luid − up ),
(30)
β = β(uf luid − up , ρf , ρp , dp , Ap , µf , ),
(31)
13
ACCEPTED MANUSCRIPT
245
with uf luid , uparticle the fluid and particle velocity respectively, ρf , ρp the respective densities, dp , Ap the
246
particle characteristic length and surface, µf the fluid viscosity. For the sake of generality, we took β as in
247
[19]
Cd =
0.5
d2p 4 π Cd ρf
24(1+0.15Re Re
0.687
)
|up − uf | ,
(32)
if Re < 1000
CR IP T
β=
0.44
(33)
if Re > 1000.
It has to be pointed out that, in the empirical correlation for the drag force of equation 32, the pressure
249
contribution cannot be separated as proposed in equation 25. A contribution caused by the pressure dis-
250
tribution around the particle due to creeping flow is taken into account in those models. However, in a
251
volume-averaged coupling this interaction is not resolved due to the choice of the spatial coarse scale. At
252
this scale, the main assumption is that the particles can be considered as zero-dimensional entities, and
253
therefore the effects of the creeping flow are modeled and not resolved. Therefore, equation 25 holds on the
254
coarse scale since here the fluid variables have to be seen as bulk pressure and bulk velocity. A further as-
255
sumption consists of considering the angular momentum exchange between particles and fluid as negligible.
256
The results of this work focus on spherical particles. However, we want to clarify that in case of particles
257
with a pronounced nonsphericity, the angular momentum exchange between fluid and particle should not
258
be neglected.
259
5.2. Phase interchange strategies
AC
CE
PT
ED
M
AN US
248
Fluid 2
Fluid 2
α=0
α=0
α=1
α=1
Fluid 1
Fluid 1
Figure 3: Particle crossing the interface between two fluids.
260
As previously mentioned, the VOF (and in general the interface capturing methods) can be extremely
261
interesting for a CFD-DEM coupling, as it uses a single velocity and single pressure formulation. This allows 14
CR IP T
ACCEPTED MANUSCRIPT
AN US
Figure 4: Phase replacement technique, particle selected for conversion to equivalent droplet. Coarse grid in red, fine grid in black.
262
simplifying the model adopted for the fluid-particle coupling since for the most part of the computational
263
domain (the one not interested by the presence of the interface) it will be reduced to a single-phase coupling.
264
Nevertheless, the interaction between particle and interface remains a complex phenomenon that needs to
265
be addressed.
First, as shown in figure 3, when a particle crosses the interface between a Fluid 1 and a Fluid 2 it will
267
not occupy volume in the Fluid 2 but instead a volume in Fluid 1. This phase interchange phenomenon
268
must be taken into account and is particularly important in dense granular flows. Then, the particle motion
269
will induce a perturbation of the interface. This phenomenon is particularly interesting to be studied
270
since a local perturbation of the interface can have a major influence on the global two-phase configuration.
271
Finally, capillary forces arise at a gas-liquid interface that can influence the particle dynamic. At a gas-liquid
272
interface, a difference of pressure will be established in the form
274
ED
PT
CE
∆p = ΓC,
(34)
with Γ the surface tension, C the interface local curvature. The capillary force will act in the normal
AC
273
M
266
direction to the surface, namely Fcap = ΓC n.
(35)
275
As mentioned in section 2, this work focuses on high Stokes number flows for which the effects of capillary
276
forces are assumed being negligible with respect to the inertia of the particles. This is reasonable when
277
operating with common values of surface tension (indicatively below 100 (10−3 N )/m). 15
ACCEPTED MANUSCRIPT
278
We here propose three different strategies for handling the phase interchange with reference to the
279
multiscale DEM-VOF method. In [10] the authors propose a porous approach for the phase interchange
280
problem in a single-scale DEM-VOF coupling. This approach consists of modifying the volume fraction
281
transport equation in order to take into account the presence of the granular phase as a porosity field,
282
leading to
CR IP T
∂α + ∇ · (αuf ) + ∇ · (α(1 − α)uc ) = 0. ∂t
(36)
The authors showed that the method can provide good results in case of packed beds of solid particles
284
interacting with a steady flow under the effect of gravity. This approach can be seen as a natural extension
285
to multiphase environments of the classic volume averaging technique for dense granular phases. As shown
286
in [10] for a sedimentation problem, this approach leads to errors in the volume conservation in the order of
287
the classical algebraic two-phase VOF on a coarse grid. This method has an obvious limitation in case of
288
high porosity (dilute) configurations since the interface motion will tend to be not affected by the particle
289
presence. In a multiscale DEM-VOF coupling, this phase interchange strategy can be used as well. The
290
porosity field defined on the coarse scale will be used as well in the modified advection equation 36. In
291
addition, one can expect the VOF scheme to be less dissipative when operated on a finer grid [7].
M
AN US
283
Another possible phase interchange strategy purely based on the definition of the color function consists
293
of re-defining the standard two-phase color function while keep using equation 15 for the phase advection.
294
In a two-phase environment, the color function represents the ratio
ED
292
Vliquid , Vtot
(37)
Vliquid . Vliquid + Vgas
(38)
297
CE
296
or
αtwo phase =
When particles enter in a two-phase environment, the total volume can be expressed as
AC
295
PT
αtwo phase =
Vtot = Vliquid + Vgas + Vparticles .
(39)
If one keeps operating on the coarse scale one can define a particle volumetric concentration as ψ=
Vparticles . Vtot
16
(40)
ACCEPTED MANUSCRIPT
298
For gas-liquid-particles mixtures, it can be useful to define a heavy-phase volume concentration as αthree phase =
299
Vliquid + Vparticles . Vliquid + Vgas + Vparticles
(41)
From the previous definitions, one can easily obtain the relation αthree phase = αtwo phase (1 − ψ) + ψ.
CR IP T
(42)
300
In this way, the volume occupied by the particle is substituted by a ghost fluid that enters directly in the
301
phase advection problem. The phase advection equation 15 will remain unchanged, while the variable itself
302
is substituted. As shown in section 7.2, this approach is extremely sensitive to the grid resolution.
Finally, using a dual-grid multiscale algorithm for the DEM-VOF coupling allows adopting a modified
304
version of the previous strategy that will be referred to as phase replacement technique. This method is
305
based on the replacement between DEM particle and an equivalent droplet. In [12] the authors, study
306
atomization processes with a VOF technique. In order to save computational cost, the authors propose to
307
substitute small droplets with Lagrangian point-particle entities, when they are far from the interface. The
308
idea we propose here is to substitute the DEM particle with a liquid droplet of the same volume when it
309
impacts the interface. In particular, as shown in figure 4, particles that are in a part of the domain interested
310
by the presence of the interface are selected for possible conversion. For those, the relative velocity between
311
the interface and the particle is compared to a reference parameter in the form
ED
M
AN US
303
(uint − up ) · nint > δC ,
313
with nint the interface normal pointing in the direction of the phase in which the particle currently is, i i C δC ∈ 0, δ∆tot , ∆C the coarse grid smallest dimension. Only the particles for which equation 43 holds will T
PT
312
(43)
be converted. After the impact, when a particle enters the second phase, it can again be described as a
315
DEM entity in its motion. As pointed out in section 1, this is one of the main advantages of the VOF
316
technique regarding the coupling with DEM. This approach allows switching between a zero-dimensional to
318
319
AC
317
CE
314
a three-dimensional momentum coupling, with a competitive computational cost. In particular, as described in section 5.1, the momentum coupling between DEM entity and fluid is zero-dimensional (based on the drag force formulation), but the interface is perturbed by a three-dimensional entity (the droplet). This technique
320
both allows implicitly taking into account the volume conservation, and modeling a three-dimensional DEM
321
entity-fluid interaction when the higher accuracy is required (i.e when the interface is perturbed by the
322
particles impact). It is worthy to point out that, in the current contribution, there is no density difference 17
ACCEPTED MANUSCRIPT
between the equivalent droplet and the heavy fluid phase. A local modification of the density of the fluid is
324
theoretically possible, but can potentially create numerical instabilities. The density difference between the
325
particle and the fluid phase will, on the other hand, be taken into account in the calculation of the particle
326
trajectory, in accordance with the standard DEM solution. As shown in the example 7.2, this technique is
327
only applicable when the supporting domain provides a fine discretization nearby the surface to properly
328
reconstruct the equivalent droplet in a VOF scheme. Therefore, the additional benefit of treating the impact
329
phenomena in a three-dimensional way is only possible in a multiscale approach.
330
6. Interpolation strategy
CR IP T
323
The dual-grid approach described in section 2 requires the definition of interpolation strategies between
332
the coarse and the fine mesh. This increases the algorithm’s complexity and requires a thoughtful considera-
333
R libraries tion. In this work, we rely on the meshT oM esh interpolation class available in the OpenFOAM
334
that are here used to map fields from the coarse to the fine mesh, and vice versa fine mesh related fields
335
R provides three different interpolation into the coarse domain. At its default implementation, OpenFOAM
336
strategies. However, it is possible to extend the existing interpolation techniques of the meshT oM esh class.
337
R are The three default strategies present in OpenFOAM
M
AN US
331
• meshT oM esh :: M AP that uses a direct mapping of the nearest cell-value, if the source mesh is
339
uniformly discretized (as in most of the case for CFD-DEM couplings) this mapping will preserve the
340
volume fraction and source of momentum. However, the resulting mapped fields will not be smooth
343
PT
342
• meshT oM esh :: IN T ERP OLAT ION uses a (gradient based) inverse distance interpolation. This is not, in general, a linear interpolation, and one can expect it not to have second-order accuracy. • meshT oM esh :: CELL P OIN T IN T ERP OLAT E performs a linear interpolation (preserving second-
CE
341
ED
338
order accuracy) by creating a tetrahedron with the center point of the cell containing the target and
345
three vertexes of the cell, and then uses an inverse distance weights to perform the linear interpolation.
346
347
AC
344
This strategy can nevertheless create problems on boundaries since the boundary value will excessively effect near-wall values.
348
A fourth strategy called CELL P OIN T F ACE IN T ERP OLAT E has been added in the class for this
349
contribution. It is still based on the construction of tetrahedron to perform linear interpolation, but executes
350
an additional check on the boundary. This allows avoiding the above-mentioned problem. With a slight
351
modification of the meshT oM esh class, it is possible to add this interpolation strategy to the options for 18
ACCEPTED MANUSCRIPT
mesh based interpolations. For the sake of simplicity, no other interpolation strategy is considered in this
353
work. In particular, if not differently specified a meshT oM esh :: M AP strategy will be used to map
354
porosity and source of momentum (Ffpi ) field from the coarse to the fine mesh, while a meshT oM esh ::
355
CELL P OIN T F ACE IN T ERP OLAT E will be used to map velocity and pressure fields from the fine
356
mesh to the coarse one. A thoughtful investigation of the optimal mapping technique for the multiscale
357
CFD-DEM method is however of high importance and will be subject for upcoming studies.
358
7. Validation and Results
359
7.1. Three-phase dam break
M
AN US
CR IP T
352
ED
Figure 5: Three-phase dam break. Setup of the simulation.
The dam break is a famous benchmark for two-phase flows simulations. The extension to a gas-liquid-
361
particle configuration represents an interesting challenge. Similar versions of the so-called “three-phase dam
362
break” have been proposed in [20, 21]. In both contributions, the authors propose a dam-break configuration
363
shown in figure 5, without intermediate obstacle, in a box of dimensions 0.2m x 0.1m x 0.3m. A column
364
of water of extension 0.05m x 0.1m x 0.1m with a uniformly distributed bed of spherical particles at the
365
bottom, breaks and stabilizes. The spheres have a diameter of 2.7 mm. We propose the multiscale simulation
366
of the benchmark using a supporting domain discretized with 384k identical cubic cells. The main point
368
CE
AC
367
PT
360
that made us choose this benchmark is that in [20, 21] the authors provide both experimental and numerical results for the test case. This makes it extremely appealing for testing a new method since advantages and
369
drawbacks of the algorithm can be directly tested. The aim is to show how, in this interface problem, the
370
usage of a finer discretization for the fluid domain can add significant advantages in the problem solution.
371
In particular, the dual-grid approach can produce grid convergent results for the continuum phase and
372
correctly capture some additional integral (kinetic energy), and local (shape of the interface) properties of 19
ACCEPTED MANUSCRIPT
373
the fluid. We adopt simulation parameters according to the authors’ suggestions as follows. The liquid
374
(heavy phase) density and viscosity are taken as 1000kg/m3 and 10−3 P a s respectively, while for the gas
375
(light phase) 1 kg/m3 and 10−5 P a s are chosen. The particle density is fixed ad 2500 kg/m3 . A linear
376
dashpot impact model is adopted for both particle-particle and particle-container interactions with spring
CR IP T
constant of 1000 N/m, a restitution coefficient of 0.9 and a friction coefficient of 0.3, as proposed in [20]. In
Liquid and solid front position as a function of Time 4 3.5
Liquid Experiments (Sun et al) Liquid MultiscaleDEM-VOF Solid Experiments (Sun et al) Solid MultiscaleDEM-VOF
AN US
Z* coord
3 2.5 2
1
0
0.5
1
M
1.5
1.5
2
2.5
3
3.5
4
ED
Nondimensional time
PT
Figure 6: Three-phase dam break. Comparison between the solution obtained with the multiscale approach and experimental results relative to the solid and liquid front motion. 377
figure 6, results related to the fluid and solid wave-front positions calculated with the multiscale DEM-VOF
379
method and the experimental results of [21] are proposed. The simulation is performed with a coarse domain
380
discretization of 6k identical cubic cells and a fine domain discretization of 384k cells. Results show a very
381
good agreement for both particles and fluid wave-front.
383
384
A second comparison with the experimental and
AC
382
CE
378
numerical results presented in the literature [20] is shown in figure 7. The discretization of the coarse and fine scale is chosen as 6k and 384k identical cubic cells respectively. It can be observed how the usage of
a finer discretization for the supporting domain (that is easily obtained in a multiscale approach), leads
385
to a better interface reconstruction. In particular, some perturbations of the gas-liquid interface can be
386
captured and the shape of the interface can be reconstructed with a higher resolution. In order to quantify
20
CR IP T
ACCEPTED MANUSCRIPT
ED
M
AN US
Figure 7: Three-phase dam break. Comparison between the solution obtained with the multiscale approach (left), the experimental setup (middle), and the DEM-VOF method [20] at time t = 0.3s , 0.4s.
CE
PT
Source Experiment Sun et al Relative error% Multiscale Relative error%
x∗A 0.887 0.804 9.4 0.882 0.56
x∗B 0.837 0.783 6.7 0.815 1.2
∗ yA 0.162 0.141 13.0 0.160 1.23
∗ yB 0.243 0.194 20.2 0.250 2.9
∗ yC 0.678 0.851 25.5 0.7001 1.3
387
AC
Figure 8: Three-phase dam break, comparison through sampling points. Points of the surface A, B, C, chosen for the comparison (top), and table comparing the nondimensional coordinates of A, B, and C. A reduction of the relative error is observed.
the improvements, we define a set of nondimensional coordinates x∗ =
388
x ∗ y ,y = , l0 l0
(44)
with l0 the length of the computational domain (0.2 m). Both the multiscale DEM-VOF and the DEM-VOF
21
ACCEPTED MANUSCRIPT
method are able to qualitatively reconstruct the shape of the interface at the left border of the domain at
390
time 0.4 s, as shown in figure 8. In order to obtain a quantitative comparison three characteristic points
391
of the curve are chosen. In particular, A and B are respectively the maximum and minimum point of the
392
first and second curve in the x − z plane, with respect to the coordinate z. C is chosen as the first point
393
at which the liquid-gas interface touches the wall. Figure 8 compares the nondimensional coordinates of
394
points A, B, and C. From the comparison it follows that a multiscale approach significantly reduces the
CR IP T
389
relative error with respect to the experimental results. In figure 9, a comparison between the nondimensional
Kinetic energy of the fluid as a function of time 0.35
AN US
48k cells(single scale) 162k cells 384k cells 500k cells
0.3 0.25 0.2
M
0.15 0.1 0.05 0
2
PT
0
ED
Nondimensional kinetic energy
0.4
4
6
8
10
12
14
Nondimensional time
395
CE
Figure 9: Three-phase dam break. Comparison between the nondimensional kinetic energy of the fluid phase obtained with fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of this quantity is observed with discretization finer than 162k cells.
total kinetic energy of the fluid phase obtained with different domain discretization is proposed. We define
397
a nondimensional kinetic energy as the ratio between the kinetic energy and the potential energy of the
398
399
AC
396
column of water at the beginning of the simulation R u · u ρ 1 dV R f f f2 Ekin = D . A ρ gdZ Zcol in h
(45)
In particular, a fine-scale discretization of 48k, 162k, 384k, and 500k identical cubic cells is adopted. In
22
ACCEPTED MANUSCRIPT
400
figure 9, the nondimensional time is defined as ∗
t =t
2g a
21
,
(46)
with a = 0.05m as in [20]. It is worthy to be noted that a structured Cartesian grid is here adopted solely
402
for the sake of simplicity since it is not a requirement for the algorithm proposed. It can be observed how
403
the kinetic energy converges with grid resolutions smaller than the particle size. Therefore a single-scale
404
approach based on the volume averaging technique, for which the finest grid-discretization is of 48k cells,
405
cannot correctly capture this quantity. In particular, the error on the kinetic energy appears higher when
406
the kinetic energy reaches local maxima (corresponding to the impact of the liquid front on the walls). This
407
can be explained considering that when higher velocities occur in the fluid, the interface dynamic becomes
CR IP T
401
0.0006 0.0004 0.0002 0
0
1
2
3
4
5
6
0.0006
500k cells 384k cells 162k cells 48k cells
0.0005 0.0004 0.0003 0.0002 0.0001 0
0
1
2
3
4
5
6
Nondimensional time
ED
Nondimensional time
M
L2 Norm of U [m5/s2]
500k cells 384k cells 162k cells 48k cells
0.0008
L2 Norm of the velocity as a function of time
L2 Norm of U-X-component [m5/s2]
L2 Norm of the velocity as a function of time 0.001
AN US
more complicated and therefore a finer discretization is required to solve them. In figure 10, the L2 norms
Figure 10: Three-phase dam break. Comparison between the L2 Norm for the velocity magnitude (left) and x component (right) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of these quantities is observed with discretization finer than 162k cells.
of the velocity vectors and the velocity x component
CE
409
PT
408
411
412
kuf x kL2 =
Z
ZD
uf · uf dV,
(47)
uf x uf x dV,
(48)
D
AC
410
kuf kL2 =
are proposed. It is interesting to notice how the kinetic energy of the flow is underestimated by the usage of
coarser grids, while the norm of the velocity field is overestimated. This can be explained considering that the two-phases have an high density ratio. Therefore, a difference in the volume fraction can induce very
413
significant changes in the global kinetic energy. A further investigation of this phenomenon is proposed in
414
figure 11, where the L2 norms of the velocity kuliquid kL2 =
Z
D
αuf · uf dV,
kugas kL2 = 23
Z
D
(1 − α)uf · uf dV,
(49)
ACCEPTED MANUSCRIPT
L2 Norm of Uliquid as a function of time
L2 Norm of Ugas as a function of time 0.0009
0.0004
48k cells(single scale) 162k cells 384k cells 500k cells
0.0008
48k cells(single scale) 162k cells 384k cells 500k cells
L2 Norm of Ugas [m5/s2]
L2 Norm of Uliquid [m5/s2]
0.0005
0.0003 0.0002 0.0001
0.0007 0.0006 0.0005 0.0004 0.0003 0.0002
0
1
2
3
4
5
6
0
7
0
1
Nondimensional time
CR IP T
0.0001 0
2
3
4
5
6
7
Nondimensional time
Volume of liquid in interface cells as a function of time 48k cells(single scale) 162k cells 384k cells 500k cells
0.0006 0.0005 0.0004 0.0003 0.0002 0.0001 0
0
1
AN US
Volume of liquid [m3]
0.0007
2
3
4
5
6
7
Nondimensional time
M
Figure 11: Three-phase dam break. Comparison between the L2 Norm for the velocity of the liquid phase (top-left), and gas phase (top-right), and volume of liquid occupying an interface cell (bottom) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Mesh independence of these quantities is observed with discretization finer than 162k cells.
are plotted as a function of time. One can observe marked differences in the gas-phase prediction that, in a
416
single-scale approach, has a peak with almost double the value of a grid convergent result. Also, in figure 11
417
the volume of the liquid occupying a cell with α ∈]1, 0[ is proposed. This shows how the adoption of a coarse
418
grid induces almost double of the fluid to be contained within interface cells. This can cause erroneous
419
predictions of the fluid behavior, as the fluid in the interface cells is treated as an equivalent mixture, and
420
therefore the solution is more sensitive to the adopted model.
421
and the volume occupied by the solid phase as functions of time are computed with different fine domain
422
discretization
AC
CE
PT
ED
415
Z
kpkL2 = pp dV, D Z Vsolid = (1 − ) dV.
In figure 12, the L2 norm of the pressure
(50) (51)
D
423
Negligible differences are observed in the pressure norm. No mesh dependency is observed for the volume
424
occupied by the solid, this shows that the projection of the porosity field from the coarse grid to the fine
425
is not sensitive to the discretization of the fine grid, as was expected due to the choice of the mapping 24
ACCEPTED MANUSCRIPT
Volume occupied by the solid as a function of time 42
500k cells 384k cells 162k cells 48k cells
80 70
41.8
60
Volume [cm3]
50 40 30
41.6 41.4
20
500k cells 384k cells 162k cells 48k cells
41.2
10 0
0
1
2
3
4
5
41
6
0
1
Nondimensional time
CR IP T
L2 Norm of the pressure [Pa2 m3]
L2 Norm of the pressure as a function of time 90
2
3
4
5
6
Nondimensional time
Figure 12: Three-phase dam break. Comparison between the L2 Norm of the pressure (left) and the volume occupied by the solid phase (right) with a fine domain discretization of 48k (single-scale), 162k, 384k and 500k cells. Minor grid dependency is noticed.
strategy. Obtaining grid convergent solutions is a very important characteristic of a numerical method,
427
especially for applications in which a complete experimental validation of the quantities of interest is not
428
possible. Therefore, this makes the proposed method competitive over a single-scale DEM-VOF method,
AN US
426
as it allows this fundamental numerical verification.
In a dual-grid approach, the mesh independence for
Liquid and solid front position as a function of Time
M
4 3.5
ED
Solid Experiments (Sun et al) Coarse 6000 Coarse 26000 Coarse 32000
PT
2.5
*
Z coord
3
CE
2 1.5
AC
1
0
0.5
1
1.5
2
2.5
3
3.5
4
Nondimensional time
Figure 13: Three-phase dam break. Comparison between the position of the particle front obtained with coarse domain discretization of 6k , 26k and 38k cells. Negligible differences are observed. 429
430
the coarse grid must also be addressed. In figure 13, the result described in figure 6 are shown with three
431
different coarse meshes of 6k, 26k, and 38k cells respectively. Negligible effects of the coarse grid choice 25
ACCEPTED MANUSCRIPT
Kinetic energy of the fluid as a function of time 0.35
6k cells(C) 26k cells(C) 38k cells(C)
0.3 0.25 0.2 0.15 0.1 0.05 0
0
2
4
6
8
10
12
14
AN US
Nondimensional time
CR IP T
Nondimensional kinetic energy
0.4
Figure 14: Three-phase dam break. Comparison between the nondimensional kinetic energy of the fluid phase obtained with coarse domain discretization of 6k , 26k and 38k cells. (Coarse) Mesh independence of this quantity is observed.
on the position of the front can be observed. This suggests that the global behavior of the particles can be
433
captured with grid discretization coarser than the one used for the classical single-scale approach regarding
434
the pure CFD-DEM coupling. In figure 14, the influence of the coarse grid choice on the kinetic energy of
435
the fluid is investigated. Three configurations with a fine-scale discretization of 160k, and the previously
436
mentioned coarse grids are confronted. As shown in figure 14, also the total kinetic energy of the fluid
437
expresses mesh independent behavior for a mesh coarser than the mesh used for the single-scale approach.
438
Being able to discretize the coarse grid with a reduced number of cells allows reducing the cost associated
439
with the calculation of the particle-related averaged fields. This can represent a significant part of the
440
cost associated with an Euler-Lagrange coupling, and is an additional important advantage of the dual-grid
441
approach. We want to remark that this dual-grid approach introduces the issue of coarse grid independence.
442
Even if the presented and other cases show coarse grid independence for grids spacing larger than the one
443
used in a single-grid approach, this property must be verified case by case. In conclusion, the ability to
445
ED
PT
CE
AC
444
M
432
separate the solution of the particle-transport scale from the one of the fluid equations is one of the main advantages of the method. In particular, it allows obtaining grid convergent solutions for the continuum
446
phase and at the same time to reduce the computational cost associated with the averaging procedure.
447
7.2. Interface perturbations induced by a single particle
448
A simple three-phase flow configuration is here proposed. It consists of a single particle falling under
449
the effect of gravity into a container partially filled with water. The volume averaging technique is normally 26
AN US
CR IP T
ACCEPTED MANUSCRIPT
(a) Top view Coarse grid (black) and supporting domain grid (red)
(b) Lateral View
M
Figure 15: One Particle falling into Water Simulation setup.
adopted to approach problems characterized by an elevated number of particles. Still, a simple setup like
451
this one can be useful in order to validate the fluid-particle coupling. We use this test case to show how the
452
multiscale approach proposed in this work is able to overcome some of the most important drawbacks of the
453
volume averaging technique concerning the DEM-VOF coupling. The domain dimensions are 0.05 m x 0.05 m
454
x 0.20 m, with the gravity in the -z direction and the initial water level at 0.051m. The spherical particle
455
has a diameter of 4mm, and is initially at rest with its center located at (0.0275m,0.0275m,0.14m). The
456
particle density is fixed at 2000kg/m3 . No-slip boundary conditions are imposed at all boundary except for
457
the top face where a total pressure boundary condition is imposed. The properties of the gas and the liquid
458
phase are chosen as the one of air and water at 25◦ C, respectively. As shown in figure 15a, for the solution
460
PT
CE
AC
459
ED
450
of the fluid-particle interaction a coarse grid discretization of 10 x 10 x 40 cubic cells has been adopted. A volume adding technique described in section 5.2 is adopted. Results are provided for the benchmark using
461
different supporting domain discretization. The aim is to show how the usage of a proper discretization for
462
the supporting domain along with the approach described in section 5.2 is able to capture interface dynamics
463
that are not resolved with a standard DEM-VOF approach. In figure 16, some significant snapshots of the
464
problem solved with a domain discretization of 500k cubic cells are presented. The particle falls in the gas 27
(b)
(c)
ED
M
AN US
(a)
CR IP T
ACCEPTED MANUSCRIPT
(d)
(e)
(f)
PT
Figure 16: One Particle falling into Water. From top to bottom and from left to right significant snapshots of the problem at time t = 0.0s (a), t = 0.13s (b), t = 0.17s(c), t = 0.18s(d), t = 0.21s(e), t = 1.0s(f), solved with a supporting domain discretization of 500k cubic cells.
phase under the effect of the gravitational force, and after t = 0.14s impacts against the liquid-gas interface
466
generating a crater structure. Around t = 0.17 the particle penetrates the interface and continues falling
467
across the liquid column with reduced velocity. The gas-liquid interface stabilizes after forming a central jet
469
470
AC
468
CE
465
at t = 0.21s. After t = 1s the particle is at rest at the bottom of the domain and the interface recovers a smooth shape (except from secondary perturbations). In figure 17a, we propose a comparison between the gas-liquid interface reconstructed at t = 0.18s with four different supporting domain discretization. The first
471
choice for the supporting domain consists of a grid with 10 x 10 x 40 identical cubic cells. This represents the
472
limit case for the multiscale approach in which the coarse domain and the supporting domain are identically
473
discretized. In this condition, the approach is equivalent to a standard DEM-VOF, since no scale separation 28
ACCEPTED MANUSCRIPT
Interface Comparison with different mesh resolution
Interface Comparison with different mesh resolution 0.07
0.054
840k 500k 108k 4k
0.065
0.05
Z-coordinate [m]
0.048 0.046 0.044 0.042
840k 500k 108k 4k
0.04 0.038
0.055 0.05
0.036 0
0.06
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
0.045
0
CR IP T
Z-coordinate [m]
0.052
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
(a) time t = 0.18s
(b) time t = 0.21s
AN US
Figure 17: One Particle falling into water. Interface comparison between a supporting domain discretization of 864k, 500k, 108k and 4k cells respectively.
is introduced. It is possible to observe how this resolution is not able to correctly capture the position of the
475
interface. This is easily explained by considering that the grid spacing is equivalent to 0.005m. As recalled
476
in section 4 the VOF schemes smear the interface in few cells, therefore the uncertainty on the surface
477
collocation will be proportional to the grid size. Using very coarse grids for the phase advection solution
478
naturally leads to high uncertainty on the interface position that, in general, can make the algorithm poorly
479
conservative. Using a finer discretization of 30 x 30 x 120 cells, the error in the interface level reduces
480
significantly, but still, the local reconstruction of the interface is not able to capture the formation of the
481
crater. The discretization of 60 x 60 x 240, on the other hand, succeeds in reproducing the stretch of the
482
interface surface with the formation of a crater structure (time 0.18 s) and a Rayleigh/Worthington internal
483
jet(time 0.21 s). In figure 17b, a comparison between the gas-liquid interface reconstructed at t = 0.21s with
484
the previously described discretization is proposed. Conclusions similar to those of figure 17a one can be
485
drawn for what concerns the first two meshes. Some more significant differences can here be noticed between
486
the 50 x 50 x 200 and the 60 x 60 x 240 discretization. In particular, the proposed discretization is not
487
sufficient to reach convergence for what concerns the simulation of the Rayleigh/Worthington internal jet.
488
This can be explained considering that the characteristic length of the jet and of the droplet detaching from
490
ED
PT
CE
AC
489
M
474
it are smaller than the one of the crater. Therefore, a finer discretization is required from the VOF scheme to properly reconstruct them. Additionally, since the volume fraction field is modified on the coarse coupling
491
scale, the resulting scheme is highly mesh sensitive. The high mesh sensitivity is a major drawback of the
492
volume adding technique. Finally, since the particle influence on the fluid velocity field is taken into account
493
as a volumetric source in the Navier-Stokes equations, it is not always able to model the correct dynamic
29
ACCEPTED MANUSCRIPT
494
at the two-phase interface. In particular, when the magnitude of the relative velocity between particle and
495
fluid is big enough to produce a significant drag force, the volume adding technique can capture surface
496
perturbations. On the other hand, in cases of reduced relative velocity in the proximity of the surface, a
497
more accurate method may be preferred. In order to show the advantages of the last method proposed in section 5.2 (the volume-replacement
499
technique), we adopt a modified version of the benchmark where the particle is positioned at the beginning
500
of the simulation slightly above the interface with an assigned initial velocity of 0.5m/s in the negative
501
z direction (this represents approximately half of the impact velocity of the previous case). The fluid is
502
considered initially at rest. This is to mimic a sudden interface crossing phenomenon that creates a sharp
503
discontinuity in the flow. As described in section 5.2, a phase replacing technique that takes full advantage
504
of the dual-grid approach, can switch between a zero-dimensional and a three-dimensional interaction with
505
the fluid when higher accuracy is required. This also allows studying the influence of the particles on the
506
fluid motion with drag-force based approaches only when the particles are actually within a single phase. In
507
figure 18, the initial configuration, as well as the streamlines of the gas phase during the interface crossing
508
are displayed. The aim is to show the flow field structures characterized by vortexes in the proximity of
AN US
CR IP T
498
In figure 19a, the crater structure formed after
CE
PT
ED
M
the particle-ghost droplet after the conversion happened.
509
510
511
AC
Figure 18: One Particle Injected into Water. Streamlines of the velocity field in the gas phase at time t = 0.0s, t = 0.01s, t = 0.03s resolved with a supporting domain discretization of 500k cubic cells.
0.03s from the particle injection is proposed. In particular, one can observe that by adopting a phase replacing technique, differences in the reconstruction of the crater obtained using a 108k and a 500k fine
512
grid discretization are negligible. Moreover, in figure 19b, it can be observed how the iso-surfaces at 0.1 of
513
the volume fraction, significantly differ when a coarser discretization is used. This is shown to remark once
514
more the introduced modeling error when adopting a VOF scheme with a coarse grid. Even if in figure 19a
515
one can argue that the 4k grid resolution is somehow able to capture the interface position, in figure 19b is 30
ACCEPTED MANUSCRIPT
Interface Comparison with different mesh resolution
Iso Surface Comparison with different mesh resolution
0.052
0.056
0.051 0.054
Z- coord
Z- coord
0.049 0.048 0.047
500k cells 108k cells 32k cells 4k cells
0.046 0.045
0.05
500k cells 108k cells 32k cells 4k cells
0.048
0.044 0.043
0.052
0.046 0
0.01
0.02
0.03
0.04
0.05
0
Y-coordinate[m]
CR IP T
0.05
0.01
0.02
0.03
0.04
0.05
Y-coordinate[m]
(a) Interface position
(b) Iso-surfaces α = 0.1 of the phase indicator
AN US
Figure 19: One Particle Injected into Water. Interface comparison between a supporting domain discretization of 500k, 108k, 32k, and 4k cells respectively at time t = 0.0.03s.
516
clearly visualized how coarse grid discretization are smearing the interface region over a significant part of
517
the simulation domain. This is in contrast to the very purpose of adopting an interface-capturing method
518
within a CFD-DEM coupling.
519
7.3. Multiple particles injection in quiescent fluid
In order to show the ability of the method to deal with the simulation of particle-laden multiphase flows
521
characterized by an elevated number of particles, we here propose a benchmark originally introduced in [10].
522
The aim of this is to show how the multiscale approach can overcome some of the DEM-VOF drawbacks
523
while still being suitable for handling simulations involving an elevated number of particles. The benchmark
524
consists of the simulation of two fluids at rest like the case from section 7.2, into which multiple particles
525
are injected causing a movement of the gas-liquid interface position. A simple configuration is proposed
PT
ED
M
520
Initial Level
AC
Case
CE
as described in [10] and shown in figure 20. Particles are introduced into the heavy phase and create a
1
0.3
2
0.3
Table 2: Parameters for the Multiple particles injection tests
Injection type
Injection duration
Continuous 1 s (1500p/s) 3000 par- Initial ticles Condition
dp
ρp
Coarse Mesh edge
Container size
0.005m
2600kg/m3 0.02m
0.1 x 0.1 x 0.5 m
0.004m
2600kg/m3 0.02m
0.1 x 0.1 x 0.5 m
526 527
porous zone at the bottom due to the effect of gravity. The inserted particles will occupy part of the volume
528
previously occupied by the pure liquid and will induce a rise of the heavy fluid’s level. The comparison 31
AN US
CR IP T
ACCEPTED MANUSCRIPT
M
Figure 20: Multiple particles injection in quiescent fluid. Initial fluid configuration (left) and comparison between the final configuration obtained with a fine discretization of 625 cells(a)), 16875 cells (b)),78125 cells (c)) and 135000 cells(d)). Negligible differences in the predicted level is observed, but significant differences in the local gradients of volume fraction can be noticed.
between the volume fraction distribution obtained with different mesh discretization is shown in figure 20.
530
As already shown in [10], the single-scale approach is capable of predicting the correct level. However, one
531
can observe how the interface region is significantly different from the one obtained with a finer fluid domain
532
discretization, and interests a significant part of the simulation domain. Furthermore, as shown in figure
533
22, the velocity norm shows mesh convergence only with mesh discretization finer than the one adopted for
534
the single-scale approach. This is in accordance with what observed in sections 7.1 and 7.2. In figure 23,
535
a second configuration of the benchmark is proposed in which the injection region is moved to the lighter
536
fluid, making the particles travel across the interface and then forming a porous zone at the bottom under
538
PT
CE
AC
537
ED
529
the action of gravity. The particles are injected as a pre-existing column at the beginning of the simulation, that travels across the interface and deposits at the bottom. The comparison between the averaged water
539
level obtained with different mesh discretization is shown in figure 25. One can notice how, in this case
540
characterized by more complex interface movement, a small error in the level is present for the single-scale
541
approach. This is in accordance to what has already been presented in [10]. However, this problem is solved
542
with the multiscale approach, in which grid convergence and level stability after the column injection are 32
ACCEPTED MANUSCRIPT
Hight of the column of water as a function of time
0.38
Theoretical Level a) b) c) d)
0.34 0.32 0.3
CR IP T
0.36
0
0.2
0.4
0.6
0.8
1
1.2
1.4
AN US
Averaged Water Level[m]
0.4
Time [s]
Figure 21: Multiple particles injection in quiescent fluid. Comparison between the theoretical average level and the ones obtained with the different mesh discretization of figure 20 . Good accordance with the analytic solution can be observed for all the cases.
L Norm of U as a function of time
M
8e-05
a) b) c) d)
ED
6e-05 5e-05
PT
4e-05 3e-05 2e-05
CE
L2 Norm of U [m5/s2]
7e-05
1e-05
AC
0
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 22: Comparison between the L2 Norm of the fluid velocity field with the domain discretization as proposed in figure 20 . Negligible grid dependence is observed for cases b), c) and d), while significant difference between these cases and the case a) is shown.
543
achieved. Once more, as shown in figure 24, the local approximation of the interface is highly affected by
544
the grid choice. In [10] and [20], the authors already pointed out how the particle size induces an important 33
ACCEPTED MANUSCRIPT
limitation for the standard DEM-VOF in the curvature calculation. In particular, in [10] they compare two
546
cases, one with particles diameter of 5 mm and one with particles diameter of 1 mm, and show how the usage
547
of small particles allows a better curvature calculation which is instead not possible for bigger particles due
548
to the nature of the DEM-VOF coupling. Similarly, in the recent work concerning the DEM-VOF approach
549
to gas-liquid-solid flows [20], the authors pointed out how the inability of the standard DEM-VOF method
550
to capture interface structures smaller than the DEM-VOF grid represents an important “drawback [...]
551
of a volume-averaging simulation”. This problem is solved by the multiscale approach proposed in this
552
contribution. As it was outlined, the key advantage of the approach is that the grid used for the interface
553
reconstruction is independent of the particle dimension. As a matter of fact, being able to numerically
554
simulate those phenomena can be valuable for many engineering processes, as the interfacial area represents
555
a key factor in many mass transport problems.
CE
PT
ED
M
AN US
CR IP T
545
556
AC
Figure 23: Injection of a column of particles inside a quiescent fluid. Relevant steps of the simulation. From left to right: initial configuration, interface perturbation due to the crossing phenomenon, particle forming a porous bed at the bottom of the domain.
In figure 26, a comparison between the L2 norms of the velocity obtained with the different mesh
557
discretization is proposed. An even more marked grid dependence than in figure 10, and 22, can be observed.
558
In particular, a peak in the norm velocity around time 0.1 s can be observed having a value almost ten times
559
higher when a single-scale approach is used. As it can be better observed in figure 27, on the long term only
560
the usage of a properly refined grid allows the L2 norm of the velocity to decay. In particular, the L2 norm 34
AN US
CR IP T
ACCEPTED MANUSCRIPT
Figure 24: Injection of a column of particles inside a quiescent fluid. Comparison between the shape of the interface obtained with different grid discretization as in figure 20. Significant differences in the shape of the interface can be noted.
Average level as a function of time
M
0.42
ED
0.38 0.36
a) b) c) d)
PT
Average level [m]
0.4
0.34
CE
0.32
AC
0.3
0
0.2
0.4
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 25: Injection of a column of particles inside a quiescent fluid: comparison between the averaged levels obtained with the grid discretization as from figure 20 . Negligible differences between case c) and case d) can be observed, while small errors in case a) and b) can be noted.
561
shows comparable behaviors only in cases c) and d). This shows once more how grid independence of the
562
solution cannot in general be granted within a single-scale DEM-VOF coupling.
563
In figure 28, we propose a comparison between the different computational cost of the proposed simula35
ACCEPTED MANUSCRIPT
L2 Norm of the Velocity as a function of time a) b) c) d)
0.8
0.4 0.2 0
0
0.2
0.4
0.6
0.8
1
CR IP T
0.6
1.2
1.4
AN US
L2 Norm of U [m2/s2]
1
1.6
Time[s]
Figure 26: Injection of a column of particles inside a quiescent fluid: Comparison between the L2 norm of the velocity obtained with the grid discretization as from figure 20 . Significant differences can be observed. Negligible differences between case c) and d) can be noted.
L2 Norm of the Velocity as a function of time (Close Up)
0.05
0
AC
M
PT
0.1
a) b) c) d)
ED
0.15
CE
L2 Norm of U [m2/s2]
0.2
0.6
0.8
1
1.2
1.4
1.6
Time[s]
Figure 27: Injection of a column of particles inside a quiescent fluid: Comparison between the L2 norm of the velocity obtained with the grid discretization as from figure 20. Close up of figure 26. Decaying behavior is observed only for case c) and d).
564
tions. This is to show how the increase in computational cost in the multiscale approach, even if extremely
565
relevant for simulations involving a single particle, becomes less important when the number of particles
566
rises. For this reason, results obtained in the cases proposed above are compared to show the dependency 36
ACCEPTED MANUSCRIPT
Single particle 1.5k particles 3k particles 500k particles
CR IP T
100
10
1
100
AN US
Clock Time/ (ClockTime SingleScale)
Normalized Clock time as a function of the number of cells
N cells / (N cells SingleScale)
Figure 28: Comparison between the increase in computational cost as a function of the refinement of the fine-scale grid, logarithmic scale. Strong dependence on the grid refinement is observed in case of Single particle, while negligible dependence is observed when a higher number of particles is adopted.
of the computational cost from the discretization used in the different cases. In order to make the various
568
scenarios comparable, we reported results for the multiscale solutions as a function of the ratio between the
569
cells adopted for the fine discretion and the ones adopted for the coarse one. For the single-scale approach,
570
this ratio will obviously be equal to unity. Similarly, for the y-axis, we defined the characteristic time as the
571
ratio between the elapsed time for the simulation and the time required by the single-scale approach. The
572
results are proposed in a logarithmic scale. One extra case is proposed that is equivalent to the one in figure
573
23, but with 500k particles of reduced diameter. For this case, only 20 fluid iterations have been solved for
574
each grid discretization and the results compared accordingly. One can see how, in cases featuring a large
575
number of particles, the computational cost is not much dependent on the grid discretization adopted. This
576
can be explained considering that in those applications, the computational burden is mainly determined by
ED
PT
CE
AC
577
M
567
the DEM part. This is normally the case for particle-laden three-phase flows, and the standard range of application for methods based on volume averaging technique. In order to better show the computational Table 3: Comparison between the execution times of the single-scale and the multiscale approaches of section 7.1
Case
Coarse Mesh
Fine Mesh
DEM δT
CFD δT
Init
1 2
None 6k
48k 500k
2x10−6 [s] 2x10−4 [s] 1.8[s] 2x10−6 [s] 2x10−4 [s] 8.4[s] 37
DEM Loop
CFD Loop
Total
6346[s] (98%) 6226[s] (73%)
130[s] (2%) 2224[s] (26%)
6476[s] 8450[s]
ACCEPTED MANUSCRIPT
intensity of the DEM and the CFD part for a standard case, we reported in table 3 the detailed data for
579
the simulation of section 7.1. The two cases confronted are a classical single-scale DEM-VOF simulation,
580
where a mesh of 48k cells and no inter-mesh interpolation are used, and the multiscale DEM-VOF operated
581
on a 6k-cells coarse grid and 500k-cells fine grid. One can observe how the increase in overall computational
582
time is rather small and generally in accordance with what is shown in figure 28. It is interesting to observe
583
how the time spent in the DEM computation, that includes the projection of the Lagrangian field on the
584
grid, is smaller in the multiscale approach if compared to the single-scale one. This shows what was antic-
585
ipated in section7.1 and figure 14 . The usage of a coarser grid for the projection of the Lagrangian fields
586
allows reducing the computational burden of this operation. On the other hand, the cost of the CFD part,
587
that includes the mapping between meshes, increases with the multiscale approach. Nevertheless, as this
588
represents a fraction of the total cost, its influence on the total computational burden of the simulation is
589
relatively small.
590
8. Conclusions
AN US
CR IP T
578
A novel multiscale approach for three-phase flows has been presented. The method is based on the
592
volume averaging technique and aims to overcome some of its main drawbacks that limit its extension to
593
complex fluid flows. The numerical scheme is capable of providing solutions for the three-phase problem
594
using an arbitrary discretization for the fluid equations. At the same time, it keeps the computational
595
advantages from coupling the particle-fluid equations with local exchanges of momentum based on drag
596
force calculation. Compared to existing methods, this multiscale approach succeeded in taking into account
597
effects of complex interface dynamics in the description of particle-laden free surface flows. Therefore, it
598
enlarges the computational window to a wide area of real case phenomena. Furthermore, the ability of the
599
multiscale DEM-VOF method to reconstruct finer structures of the flow, and in particular of the gas-liquid
600
interface, makes it suitable to study transport phenomena more accurately. The method has been shown
601
to produce grid convergent solutions. This allows providing an important insight on the correctness of
603
604
ED
PT
CE
AC
602
M
591
the solution. It has been shown how the multiscale DEM-VOF approach is able to agree with benchmark problems presented in the literature, and in particular is able to reproduce experiments with a reduced error if compared to the standard DEM-VOF. Furthermore, the method allows adopting a wider range of volume
605
replacement algorithms, opening the possibility to capture phenomena of bubble and droplet formation, in
606
accordance with the classical VOF theory. Finally, the increase in computational cost of the new method
607
with respect to the classical single-scale approach has been shown to be negligible in cases involving by an 38
ACCEPTED MANUSCRIPT
608
elevated number of particles. The main limitation of the method consists in the higher complexity of the algorithm, that requires a
610
mapping operation between the different grids. This may make its extension to cases involving reacting
611
flows, and heat and mass transfer between particles and fluid, more complex, as the requirements for the
612
mapping operations can become more delicate. Additionally, an efficient parallelization of the algorithm will
613
require a dedicated study due to its enhanced complexity. As a future investigation, the influence of different
614
interpolation strategies on the convergence properties of the algorithm can be studied. Furthermore, the
615
extension of the approach to turbulent configurations is of high interest. In particular, the proposed dual-grid
616
algorithm here proposed brings significant advantages, as the flexibility in the fluid domain discretization
617
can play a fundamental role in the reconstruction of boundary layers and turbulent structures. However,
618
the introduction of sub-grid models in the current algorithm is a topic that must be carefully addressed,
619
with particular reference to the influence of those models on the interface dynamics, and the treatment of
620
sub-grid terms involving the surface tension.
AN US
CR IP T
609
The extension of the multiscale approach presented in this paper to CFD-DEM coupling different from
622
the DEM-VOF is another potential outcome of the present work. It is our opinion that several other
623
couplings between DEM and complex flow can benefit from it.
624
Acknowledgments
ED
625
M
621
The experiments presented in this paper were carried out using the HPC facilities of the University of Luxembourg [24].
627
fruitful discussions. Finally, we would like to thank Florian Hoffmann, Jack S. Hale, and Amir Houshang
628
Mahmoudi for their invaluable help and suggestions in the last revision of the manuscript.
629
References
633 634 635 636 637 638 639 640 641 642 643 644
CE
632
[1] Afkhami, M., Hassanpour, A., Fairweather, M., Njobuenwu, D., 2015. Fully coupled LES-DEM of particle interaction and agglomeration in a turbulent channel flow. Computers & Chemical Engineering 78, 24–38. URL http://www.sciencedirect.com/science/article/pii/S0098135415001040 [2] Afkhami, S., Zaleski, S., Bussmann, M., 2009. A mesh-dependent model for applying dynamic contact angles to VOF simulations. J. Comput. Physics 228 (15), 5370–5389. URL http://dblp.uni-trier.de/db/journals/jcphy/jcphy228.html#AfkhamiZB09 [3] Anderson, T. B., Jackson, R., 1967. Fluid Mechanical Description of Fluidized Beds. Equations of Motion. Industrial & Engineering Chemistry Fundamentals 6 (4), 527–539. URL http://dx.doi.org/10.1021/i160024a007 [4] Berberovi´ c, E., van Hinsberg, N. P., Jakirli´ c, S., Roisman, I. V., Tropea, C., Mar 2009. Drop impact onto a liquid layer of finite thickness: Dynamics of the cavity evolution. Phys. Rev. E 79, 036306. URL http://link.aps.org/doi/10.1103/PhysRevE.79.036306 [5] Du, W., Bao, X., Xu, J., Wei, W., 2006. Computational fluid dynamics (CFD) modeling of spouted bed: Assessment of drag coefficient correlations. Chemical Engineering Science 61 (5), 1401–1420. URL http://www.sciencedirect.com/science/article/pii/S0009250905006603
AC
630 631
We wish to thank the anonymous reviewers and the editor for the careful work and
PT
626
39
ACCEPTED MANUSCRIPT
653 654 655 656 657 658 659 660 661 662 663 664 665 666 667 668 669 670 671 672 673 674 675 676 677 678 679 680 681 682 683 684 685 686 687 688 689 690 691 692 693 694 695 696 697 698 699 700 701 702 703 704 705
CR IP T
652
AN US
651
M
650
ED
649
PT
648
CE
647
[6] Herrmann, H. J., 1995. Physics of granular media. Chaos, Solitons & Fractals 6, 203–212. URL http://www.sciencedirect.com/science/article/B6TJ4-49WMG9H-X/1/f106bb7c6924a76a8a0c156a99a07718 [7] Hirt, C., Nichols, B., 1981. Volume of fluid (VOF) method for the dynamics of free boundaries. Journal of Computational Physics 39 (1), 201–225. URL http://www.sciencedirect.com/science/article/pii/0021999181901455 [8] H.Jasak, 1996. Error analysis and estimation for the finite volume method with applications to fluid flows. Ph.D. thesis. [9] Jasak, H., Weller, H. G., Gosman, A. D., Sep. 1999. High resolution NVD differencing scheme for arbitrarily unstructured meshes. International Journal for Numerical Methods in Fluids 31, 431–449. [10] Jing, L., Kwok, C. Y., Leung, Y. F., Sobral, Y. D., 2016. Extended CFD–DEM for free-surface flow with multi-size granules. International Journal for Numerical and Analytical Methods in Geomechanics 40 (1), 62–79, nAG-14-0182.R1. URL http://dx.doi.org/10.1002/nag.2387 [11] Li, Y., Zhang, J., Fan, L.-S., 1999. Numerical simulation of gas–liquid–solid fluidization systems using a combined CFDVOF-DPM method: bubble wake behavior. Chemical Engineering Science 54 (21), 5101–5107. URL http://www.sciencedirect.com/science/article/pii/S0009250999002638 [12] Ling, Y., Zaleski, S., Scardovelli, R., 2015. Multiscale simulation of atomization with small droplets represented by a Lagrangian point-particle model. International Journal of Multiphase Flow 76, 122–143. URL http://www.sciencedirect.com/science/article/pii/S0301932215001524 [13] Mahmoudi, A. H., Hoffmann, F., Markovic, M., Peters, B., Brem, G., 2016. Numerical modeling of self-heating and selfignition in a packed-bed of biomass using {XDEM}. Combustion and Flame 163, 358–369. URL http://www.sciencedirect.com/science/article/pii/S0010218015003582 [14] Mahmoudi, A. H., Hoffmann, F., Peters, B., 2014. Detailed numerical modeling of pyrolysis in a heterogeneous packed bed using {XDEM}. Journal of Analytical and Applied Pyrolysis 106, 9–20. URL http://www.sciencedirect.com/science/article/pii/S0165237013002660 [15] Mahmoudi, A. H., Hoffmann, F., Peters, B., 2016. Semi-resolved modeling of heat-up, drying and pyrolysis of biomass solid particles as a new feature in {XDEM}. Applied Thermal Engineering 93, 1091–1104. URL http://www.sciencedirect.com/science/article/pii/S1359431115010807 [16] Peters, Bernhard, Pozzetti, Gabriele, 2017. Flow characteristics of metallic powder grains for additive manufacturing. EPJ Web Conf. 140, 13001. URL https://doi.org/10.1051/epjconf/201714013001 [17] Pozzetti, G., Peters, B., 2017. On the choice of a phase interchange strategy for a multiscale DEM-VOF Method. AIP Conference Proceedings 1863. URL http://orbilu.uni.lu/handle/10993/28863 [18] Rui, S., Heng, X., 2015. Diffusion-based coarse graining in hybrid continuum–discrete solvers: Applications in CFD–DEM. International Journal of Multiphase Flow 72, 233–247. [19] Schiller, L., Naumann, Z., 1935. A Drag Coefficient Corre-lation,. VDI Zeitung Vol. 77,, 318–320. [20] Sun, X., Sakai, M., 2015. Three-dimensional simulation of gas–solid–liquid flows using the DEM–VOF method. Chemical Engineering Science 134, 531–548. URL http://www.sciencedirect.com/science/article/pii/S0009250915003991 [21] Sun, X., Sakai, M., Yamada, Y., 2013. Three-dimensional simulation of a solid–liquid flow by the DEM–SPH method. Journal of Computational Physics 248, 147–176. URL http://www.sciencedirect.com/science/article/pii/S0021999113002684 [22] Tryggvason, G., Scardovelli, R., Zaleski, S., January 2011. Direct Numerical Simulations of Gas-Liquid Multiphase Flow. Cambridge University Press. [23] Uhlmann, M., 2005. An immersed boundary method with direct forcing for the simulation of particulate flows. Journal of Computational Physics 209 (2), 448–476. URL http://www.sciencedirect.com/science/article/pii/S0021999105001385 [24] Varrette, S., Bouvry, P., Cartiaux, H., Georgatos, F., July 2014. Management of an Academic HPC Cluster: The UL Experience. In: Proc. of the 2014 Intl. Conf. on High Performance Computing & Simulation (HPCS 2014). IEEE, Bologna, Italy, pp. 959–967. [25] Weller, H. G., Tabor, G., Jasak, H., Fureby, C., 1998. A tensorial approach to computational continuum mechanics using object-oriented techniques. Computers in Physics 12 (6), 620–631. URL http://scitation.aip.org/content/aip/journal/cip/12/6/10.1063/1.168744 [26] Wu, T.-R., Chu, C.-R., Huang, C.-J., Wang, C.-Y., Chien, S.-Y., Chen, M.-Z., 2014. A two-way coupled simulation of moving solids in free-surface flows. Computers & Fluids 100, 347–355. URL http://www.sciencedirect.com/science/article/pii/S0045793014002023 [27] Xiao, H., Sun, J., 2010. Algorithms in a Robust Hybrid CFD-DEM Solver for Particle-Laden Flows. Communication Computer Physics. [28] Zhu, H., Zhou, Z., Yang, R., Yu, A., 2007. Discrete particle simulation of particulate systems: Theoretical developments. Chemical Engineering Science 62 (13), 3378–3396, frontier of Chemical Engineering - Multi-scale Bridge between Reductionism and Holism. URL http://www.sciencedirect.com/science/article/pii/S000925090700262X
AC
645 646
40