∎

11institutetext: 1Laboratory of Fluid Dynamics and Technical Flows, University
of Magdeburg “Otto von Guericke”, Magdeburg, 39106, Germany.
2Department of Mechanical and Process Engineering, ETH Zürich, 8092, Zürich,
Switzerland.
3 Interdisciplinary Centre for Advanced Materials Simulation (ICAMS), Ruhr-
University Bochum, Bochum, 44801, Germany.

# Modeling gas flows in packed beds with the lattice Boltzmann method:
validation against experiments

Tanya Neeraj*1 Corresponding author: Tanya Neeraj<EMAIL_ADDRESS>Christin Velten1 Gabor Janiga1 Katharina Zähringer1 Reza Namdar3 Fathollah
Varnik3 Dominique Thévenin1 Seyed Ali Hosseini1,2

(Received: date / Accepted: date)

###### Abstract

This study aims to validate the lattice Boltzmann method and assess its
ability to accurately describe the behavior of gaseous flows in packed beds.
To that end, simulations of a model packed bed reactor, corresponding to an
experimental bench, are conducted, and the results are directly compared with
experimental data obtained by Particle Image Velocimetry measurements. It is
found that the lattice Boltzmann solver exhibits very good agreement with
experimental measurements. Then, the numerical solver is further used to
analyze the effect of the number of packing layers on the flow structure and
to determine the minimum bed height above which the changes in flow structure
become insignificant. Finally, flow fluctuations in time are discussed. The
findings of this study provide valuable insights into the behavior of the gas
flow in packed bed reactors, opening the door for further investigations
involving additionally chemical reactions, as found in many practical
applications.

###### Keywords:

Lattice Boltzmann Packed bed Body-Centered-Cubic packing Computational Fluid
Dynamics (CFD) Particle Image Velocimetry (PIV)

##  1 Introduction

Packed bed reactors are an essential component of many engineering, chemical
and petrochemical processes. Their relevance lies in their ability to provide
a large surface area for mass transfer, high pressure drop, and high heat
transfer rates. Investigating these packed beds provides valuable insights
into their performance and efficiency, crucial for optimizing their design and
operation. The flow behavior, pressure drop, heat transfer, and mass transfer
characteristics of packed beds are being analyzed and this information then
aids in predicting their performance under different operating conditions and
facilitates the development of more efficient and cost-effective industrial
processes. These reactors consist of a fixed bed involving many solid
particles, often reacting, in contact with the surrounding liquid or gas
phase. In the present study, only gas flows will be considered, since they
correspond to the high-temperature applications (drying, heat storage,
roasting, calcination, pyrolysis…) driving the present investigations. Packed
bed reactors allow for efficient heat and mass transfer and – as a consequence
– reactions between gas and solid phase. The combustion and reaction processes
in these packed bed reactors are intricate and involve a number of coupled
physical and chemical phenomena – most prominently mass transfer, heat
transfer, and fluid dynamics (Bao and Meskas, 2011; Dixon and Partopour, 2020;
Young-Il, 2021; Eppinger et al, 2014).

Unfortunately, the full experimental characterization of these packed bed
geometries faces several challenges including the complex and three-
dimensional nature of the beds, tracking individual particles for moving beds,
capturing non-uniformities and heterogeneities of the flow within the bed,
scaling up experimental findings, and ensuring reproducibility. Additionally,
high-temperature processes in packed beds cannot be investigated
experimentally with most existing diagnostics, either because of the very
difficult access to the inter-particle spaces and/or of the harsh conditions
found there. Getting accurate experimental data regarding the exchange
processes at this level is essential to optimize further all underlying
applications, improving energy efficiency and reducing environmental impact
(Bao and Meskas, 2011; Alobaid et al, 2022).

As a complementary source of information, sophisticated mathematical and
computational models should be used to accurately simulate the behavior of
packed bed reactors. With the advent of computational fluid dynamics (CFD) and
constant progress in computational power, there is a growing interest in
modeling numerically all kinds of multiphase flows (Sommerfeld et al, 2008;
Wan et al, 2008; Wang et al, 2019), supporting also design, operation, and
optimization of packed bed reactors (Illana Mahiques et al, 2023; Alobaid et
al, 2022; Freund et al, 2003; Augier et al, 2010). Numerical methods, such as
the lattice Boltzmann method (LBM), can help in the experimental
characterization of these packed bed geometries by providing a complementary
approach to experimental techniques. LBM is an approach that is rapidly
gaining popularity due to its ability to simulate a wide variety of complex
flows. It is extremely well-suited for describing fluid flows in porous media
(Bao and Meskas, 2011; Chen and Doolen, 1998; Eshghinejadfard et al, 2016). It
simulates the movement of particles on a discrete lattice (Sudhakar and Das,
2020; Yang and Boek, 2013) and can provide information about the flow
behavior, velocity distribution, and mass transfer characteristics within the
packed bed at the pore-scale level. LBM simulations excel in capturing the
complexities of three-dimensional packed bed structures and provide a means to
visualize internal flow patterns that are inaccessible experimentally. It can
also be used to optimize the design and operation of packed beds by predicting
their performance under different operating conditions. This method is
particularly efficient for complex geometries due to a very easy grid
generation and treatment of boundary conditions. Additionally, LBM simulations
can be scaled up to larger bed sizes (as done in the present study), aiding in
the prediction of packed bed performance in industrial settings. In summary,
LBM is in principle very simple, shows excellent scalability on parallel
supercomputers, and extensions have been proposed to represent accurately all
kinds of flows (Sudhakar and Das, 2020; Krüger et al, 2011; Ya-Ling et al,
2019; Boivin et al, 2021; Zhao et al, 2020). Therefore, by integrating
numerical methods like LBM with experimental data, we can gain a more
comprehensive understanding of packed bed behavior, optimize their design, and
enhance the efficiency of the industrial processes.

However, to fully leverage the benefits of simulations, it is imperative to
ensure the accuracy and reliability of the results obtained. Even with
sophisticated models and computational techniques, there is a need to
establish trust in numerical outcomes. This necessitates the combination of
experimental and numerical studies for benchmarking, as experimental
validation serves as a critical reference for assessing the accuracy of
numerical simulations. Here also, many challenges must be solved to ensure
sufficiently accurate predictions for such complex two-phase flows while
keeping computational requirements at an acceptable level. Due to non-linear
coupling processes between momentum, heat, and mass transfer, very slight
changes in geometry at small scale or minute modifications of the operating
conditions may have a significant impact on the behavior of the whole reactor
(Prokopová and Prokop, 2009; Gernaey et al, 2015).

In recent decades, the lattice Boltzmann method (LBM) has proved to be an
interesting numerical approach to address this type of phenomena. Rong et al
(2014), for example, used LBM to explore the effects of particle size
distribution on flow behavior and drag forces within binary mixtures of
particles. They identified flow characteristics related to pore geometries and
proposed a new equation for calculating mean individual drag force. Another
study by the same authors Rong et al (2015) used LBM to simulate fluid flow
through packed beds of uniform ellipsoids and investigated the impact of
particle shape, orientation, and volume fraction on flow behavior and pressure
drop. Sullivan et al (2005) on the other hand, explored the application of
lattice Boltzmann methods for simulating hydrodynamics, reaction, and mass
transfer in a packed bed reactor of catalyst particles at pore-scale.
Spatially varying values of diffusivity have been considered to differentiate
between intra- and inter-particle diffusion. A similar study used LBM to
simulate concentration profiles in a fixed bed filled with spherical porous
adsorbents (Verma et al, 2007). The central objective was to investigate how
concentration profiles change due to inter- and intra-particle mass transport
in an adsorber with small tube-to-sphere diameter ratios, down to less than
10. Various packing arrangements and different pore diffusivities were
considered, validated with experimental data obtained by near-infrared optical
tomography in a tubular adsorber packed with zeolite particles. Tariq and Liu
(2022) investigated flow structures at the level of the inter-particle space,
without discussions of the overall flow features within and above a full bed
of particles.

It is noteworthy that, despite the large variety of lattice Boltzmann studies
of the flow through complex geometries including packed beds, studies, which
also include an experimental validation of numerical results are quite sparse.
Combined experimental-numerical investigations proposed in the present study
thus aims at filling this gap. Considering the complexity of this problem,
combining experimental and numerical investigations of exactly the same
configuration appears to be the best solution in order to push further our
understanding. This is the central objective of the present study. It was
already stressed that experimental measurements are extremely challenging due
to limited access and high temperatures. Avoiding in this first investigation
the latter issue by considering only a two-phase cold flow, a solution must
still be found to measure accurately flow properties in-between particles
within a packed bed.

Again, at the difference of most published studies, the present investigation
concentrates only on gaseous flows through packed beds. To the best of our
knowledge, no publication using LBM considered such gas flows in large-scale
packed beds up to now. Therefore, and even if the in-house LBM solver employed
in this study (ALBORZ, described in particular in Hosseini (2020); Hosseini et
al (2020a, 2021)) has been previously validated for a variety of
configurations, it is particularly important to quantify again its accuracy
and reliability by direct comparisons with experimental data for the same
packed-bed configuration. In what follows, the air flow within a body-
centered-cubic (BCC) packing will be compared to results from PIV measurements
regarding flow structures and velocities, both directly above as well as
within the packed bed. This will deepen our understanding of the flow behavior
inside the reactor in order to push further process design, ultimately
enhancing performance and efficiency. Two different particle Reynolds numbers
will be examined in a range relevant for practical applications. In addition,
the numerical results will be analyzed to better understand the overall
process as well as the flow behavior at selected positions. The occurrence of
flow fluctuations and instabilities will finally be discussed. The LBM results
show good agreement with experimental data, providing a strong basis for the
numerical model developed in this study.

##  2 Numerical solver ALBORZ: the lattice Boltzmann method

###  2.1 Modified central Hermite multiple-relaxation time lattice Boltzmann
solver

The lattice Boltzmann method is a solver for the Boltzmann equation in the
limit of the hydrodynamic regime (He and Luo, 1997), tailored to recover the
Navier-Stokes equations. It can be derived by first discretizing the Boltzmann
equation in phase space, the space of particle degrees of freedom, via
projection onto Hermite polynomials and the use of the Gauss-Hermite
quadrature. The resulting system of coupled hyperbolic equations is then
integrated along characteristics (Hosseini and Karlin, 2023). This in turn
yields the now-famous ”stream-collide” equation:

$f_{i}\left(\bm{r}+\bm{c}_{i}\delta t,t+\delta
t\right)-f_{i}\left(\bm{r},t\right)=\Omega_{i},$ (1)

for discrete distribution functions $f_{i}$, where $\bm{c}_{i}$ are discrete
particle velocities, $\Omega_{i}$ is the discrete collision operator, modeling
equilibration of the distribution function due to molecular collisions, and
$\delta t$ is the time-step. Close to local thermodynamic equilibrium (LTE)
conditions (which is valid for small Knudsen numbers, as found in all flows
considered here) a linear relaxation operator (Bhatnagar-Gross-Krook, BGK) is
usually employed, which in its discrete form reads

$\Omega_{i}=\omega\left(f^{\rm eq}_{i}-f_{i}\right)+\Xi_{i},$ (2)

where $\omega$ is the relaxation frequency coefficient and $f^{\rm eq}_{i}$ is
the discrete form of the continuous attractor obtained via minimization of
entropy, i.e., the Maxwell-Boltzmann distribution (Chapman and Cowling, 1990).
It must be noted that here the collision operator is supplemented with a term
$\Xi_{i}$ accounting for deviations in the diagonal components of the third-
order moments of the equilibrium term (for models relying on third-order
quadratures). Details on this term can be found in Feng et al (2015);
Prasianakis and Karlin (2007); Hosseini et al (2020b), among other sources. In
the context of the lattice Boltzmann method the discrete equilibrium
distribution function is expressed as a finite-order expansion of the Maxwell-
Boltzmann distribution in terms of Hermite polynomials by means of the Gauss-
Hermite quadrature:

$f^{\rm
eq}_{i}=w_{i}\sum_{n=0}^{N}\frac{1}{n!c_{s}^{2n}}\mathcal{H}_{n}(\bm{c}_{i}):a_{n}^{\rm
eq}(\rho,\bm{u},p/\rho),$ (3)

where $\mathcal{H}_{n}$ is the Hermite polynomial tensor of order/rank $n$,
$a_{n}^{\rm eq}$ the corresponding equilibrium coefficient, $w_{i}$ are
weights associated with the quadrature and $c_{s}$ the lattice speed of sound,
also tied to the quadrature. Further quantities appearing in the system are
$\rho$ the fluid density, $\bm{u}$ the velocity, $p$ the pressure and $N$ the
order of expansion. For the remainder of the manuscript we will restrict
ourselves to the D2Q9 and D3Q27 stencils, i.e. third-order quadratures, where
$D$ is the physical dimension (2D or 3D), and Q the number of discrete
velocities. In two dimensions, $N$ is set to four, while in 3D it is set to
six. Both are dictated by the maximum number of terms supported by the
stencils. Contrary to classical second-order polynomial expansions, the
present model allows the discrete form to preserve Galilean invariance of the
dissipation rate of shear modes at the Navier-Stokes level, which is an
essential aspect. The term $\Xi_{i}$ restores Galilean invariance to the
dissipation rate of normal modes and is defined as (Hosseini et al, 2020b;
Feng et al, 2015):

$\Xi_{i}=\left(1-\frac{\omega}{2}\right)\frac{w_{i}}{2c_{s}^{4}}\left(\bm{\nabla}\mathcal{H}_{2}\right):\Delta^{\rm
eq}_{3},$ (4)

where $\Delta^{\rm eq}_{3}$ is a diagonal tensor of rank 3 defined as:

$\Delta^{\rm eq}_{\alpha\alpha\alpha}=\rho
u_{\alpha}\left(u_{\alpha}^{2}+3p/\rho\right)-\sum_{i}c_{i,\alpha}^{3}f^{\rm
eq}_{i}.$ (5)

More advanced versions of the BGK operator are usually employed on the basis
of either physical, for instance variable Prandtl number, or numerical
arguments. Among those, the multiple relaxation time (MRT) formulation allows
to relax the distribution function in a space of linearly independent base
functions, typically moments of the distribution function. Owing to its
improved stability and accuracy it has become quite popular in the lattice
Boltzmann community. Here, we make use of a MRT in a space of modified Hermite
central moments allowing for independent control over the bulk viscosity, as
first described in Hosseini et al (2022c). At the difference of the classical
Hermite polynomial space (Hosseini et al, 2021; Huang et al, 2022), the trace
and trace-free contributions to the second-order moments relax independently
in this formulation:

$\Omega_{i}=\mathcal{T}^{-1}W\mathcal{T}(f^{\rm eq}_{i}-f_{i})+\Xi_{i}$ (6)

where $\mathcal{T}$ and $\mathcal{T}^{-1}$ are the moments transform tensor
and its inverse, and $W$ is the diagonal tensor of relaxation rates. The
moments transform tensor is defined via a set of modified central Hermite
polynomials $\widetilde{\mathcal{H}}_{n}(\bm{c}_{i})$ , which for the D3Q27
stencil are:

$\widetilde{\mathcal{H}}_{n}\in\\{\widetilde{\mathcal{H}}_{0},\widetilde{\mathcal{H}}_{x},\widetilde{\mathcal{H}}_{y},\widetilde{\mathcal{H}}_{z},\widetilde{\mathcal{H}}_{xy},\widetilde{\mathcal{H}}_{xz},\widetilde{\mathcal{H}}_{yz},\widetilde{\mathcal{H}}_{xx}-\widetilde{\mathcal{H}}_{yy},\widetilde{\mathcal{H}}_{xx}-\widetilde{\mathcal{H}}_{zz},\widetilde{\mathcal{H}}_{xx}+\widetilde{\mathcal{H}}_{yy}+\widetilde{\mathcal{H}}_{zz},\\\
\widetilde{\mathcal{H}}_{xyz},\widetilde{\mathcal{H}}_{xxy},\widetilde{\mathcal{H}}_{xxz},\widetilde{\mathcal{H}}_{xyy},\widetilde{\mathcal{H}}_{yyz},\widetilde{\mathcal{H}}_{xzz},\widetilde{\mathcal{H}}_{yzz},\widetilde{\mathcal{H}}_{xxyy},\widetilde{\mathcal{H}}_{xxzz},\widetilde{\mathcal{H}}_{yyzz},\widetilde{\mathcal{H}}_{xxyz},\\\
\widetilde{\mathcal{H}}_{xyyz},\widetilde{\mathcal{H}}_{xyzz},\widetilde{\mathcal{H}}_{xyyzz},\widetilde{\mathcal{H}}_{xxyzz},\widetilde{\mathcal{H}}_{xyzz},\widetilde{\mathcal{H}}_{xxyyzz}\\}.$
(7)

For detailed expression of central Hermite polynomials we refer interested
readers to Hosseini (2020). Equilibrium moments based on these polynomials are
the central Hermite coefficients $\widetilde{a}^{\rm eq}_{n}$, which for the
equilibrium distribution function using D3Q27 are:
$\widetilde{a}^{\rm eq}_{n}\in\\{\rho,0,0,0,0,0,0,0,0,3(p-\rho
c_{s}^{2}),0,0,0,0,0,0,0,{(p-\rho c_{s}^{2})}^{2},\\\ {(p-\rho
c_{s}^{2})}^{2},{(p-\rho c_{s}^{2})}^{2},0,0,0,0,0,0,{(p-\rho
c_{s}^{2})}^{3}\\}.$ (8) Setting $p=\rho c_{s}^{2}$, all equilibrium
coefficients reduce to zero except at order zero. The relaxation frequency
tensor $W$ for the D3Q27 stencil is defined as:
$W={\rm
diag}(1,1,1,1,\omega_{s},\omega_{s},\omega_{s},\omega_{s},\omega_{s},\omega_{b},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\\\
\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g},\omega_{g}),$
(9) where the operator ${\rm diag}$ is defined as:

${\rm diag}(\bm{A})=(\bm{A}\otimes\bm{1})\circ\bm{I},$ (10)

for a given vector $\bm{A}$, with $\bm{1}$ a vector with elements 1, $\bm{I}$
the unitary tensor, and $\circ$ the Hadamard product, while $\omega_{s}$,
$\omega_{b}$ and $\omega_{g}$ are the shear, bulk, and ghost-modes relaxation
frequency, respectively. The quantity $\omega_{b}$ is related to the bulk
viscosity, $\eta$, as Hosseini et al (2022b):

$\omega_{b}=\frac{\delta t}{(\frac{2+D}{D}-\frac{\partial\ln
p}{\partial\ln\rho})\frac{\eta}{p}+\delta t/2},$ (11)

while the shear modes relax with:

$\omega_{s}=\frac{\delta t}{\frac{\mu}{p}+\delta t/2},$ (12)

and $\omega_{g}=1$. The correction term for the diagonal components of the
third-order moments appearing in Eq. 6 changes in the context of the multiple
relaxation time collision operator into:

$\Xi_{i}=\left(1-\frac{\omega_{b}\omega_{s}}{\omega_{b}+\omega_{s}}\right)\frac{w_{i}}{2c_{s}^{4}}\left(\bm{\nabla}\mathcal{H}_{2}\right):\Delta^{\rm
eq}_{3}.$ (13)

The lattice Boltzmann model as introduced here allows for an independent bulk
viscosity and Galilean-invariant dynamic and bulk viscosities. It is worth
noting that all the LBM equations mentioned in this article are expressed in
S.I. units and no non-dimensionalization has been used.

###  2.2 Boundary conditions

While a variety of different boundary conditions have been developed for LBM,
in the context of the present study all solid boundaries can be simply modeled
using the half-way bounce-back scheme, where missing distribution functions
after the collision-streaming steps are computed as (Krüger et al, 2017):

$f_{i}\left(\bm{x},t+\delta t\right)=f^{*}_{\bar{i}}\left(\bm{x},t\right),$
(14)

where $f^{*}$ is the post-collision population (prior to streaming) and
$\bar{i}$ is the index of the particle velocity opposite that of $i$. To get
rid of any staircase approximation, the half-way bounce-back is supplemented
with the advanced treatment for curved boundaries proposed in Bouzidi et al
(2001).
At a given boundary node $\bm{x}_{f}$, the missing incoming populations are
computed as:

$\displaystyle f_{i}(\bm{x}_{f},t+\delta t)$
$\displaystyle=2qf_{\bar{i}}(\bm{x}_{f}+\bm{c}_{\bar{i}},t+\delta t)$
$\displaystyle+\left(1-2q\right)f_{\bar{i}}(\bm{x}_{f},t+\delta t),\forall
q<\frac{1}{2},$ (15a) $\displaystyle f_{i}(\bm{x}_{f},t+\delta t)$
$\displaystyle=\frac{1}{2q}f_{\bar{i}}(\bm{x}_{f}+\bm{c}_{\bar{i}},t+\delta
t)$ $\displaystyle+\frac{2q-1}{2q}f_{\bar{i}}(\bm{x}_{f},t+\delta t),\forall
q\geq\frac{1}{2},$ (15b)

where $\bar{i}$ designates the direction opposite $i$ and $q$:

$q=\frac{\lvert\lvert\bm{x}_{f}-\bm{x}_{s}\lvert\lvert}{\lvert\lvert\bm{c}_{i}\lvert\lvert},$
(16)

with $\bm{x}_{s}$ denoting the wall position in direction $i$.
Other boundary conditions applied at the inlets and outlets, i.e., constant
velocity and/or constant pressure, are realized using the non-equilibrium
extrapolation approach proposed in Zhao-Li et al (2002). In this method the
missing populations at the boundary nodes are reconstructed as:

$f_{i}\left(\bm{x}_{w},t+\delta t\right)=f^{\rm
eq}_{i}\left(\bm{x}_{w},t+\delta t\right)+f^{\rm
neq}_{i}\left(\bm{x}_{w},t+\delta t\right),$ (17)

where $\bm{x}_{w}$ is the position of the boundary node and $f^{\rm neq}$ is
the non-equilibrium part of the distribution function. To be able to use this
equation one needs both, density and velocity at the boundary node. One of
these two variables is known via the boundary condition, and the other one is
computed by interpolation from neighboring nodes. The non-equilibrium part of
the distribution function is also computed via interpolation in space.
In addition to velocity and pressure boundary conditions, the zero-gradient
outlet boundary condition was implemented following Krüger et al (2017):

$f_{i}\left(\bm{x},t+\delta t\right)=f_{i}\left(\bm{x}-\bm{n}\delta x,t+\delta
t\right),$ (18)

where $\bm{n}$ is the outward-pointing unit vector normal to the boundary
surface.

One attractive feature of LBM is that, using such boundary conditions, there
is no need for any complex grid adaption to resolve the structure of the
packed bed. A simple, regular, equidistant grid is used in all simulations.
The employed resolution in time and space is discussed later on.

##  3 Configuration considered in this study

As explained previously, the same geometry, same size, same operating
conditions have been considered for numerical simulations and experimental
measurements, enabling direct comparisons discussed later. However, real
experiments and diagnostics involve of course many more aspects needed for a
practical realization. This is why the next sections describe separately
particular features of the experimental bench and of the LBM simulations.

Practical applications of interest for our group involve typical particle
Reynolds numbers $\mathrm{Re}_{P}$ between 100 and 1000, in most cases between
300 and 500. This is why these two particular values have been considered in
the present study. Here the Reynolds number is defined as:

$\mathrm{\mathrm{Re}_{P}}=\frac{\rho d_{p}v_{int}}{\mu},$ (19)

where $d_{p}$ is particle diameter, $\mu$ is the dynamic viscosity of air for
ambient pressure and temperature ($p=1~{}\mathrm{bar}$ and $T$ = 20 °C),
$\mu=1.822\cdot 10^{-5}\,\mathrm{Pa\cdot s}$ and $\rho$ is the air density
under the same conditions, $\rho=1.204\,\mathrm{kg/m^{3}}$. The interstitial
velocity $v_{int}$ is defined as:

$v_{\mathrm{int}}=\frac{v_{s}}{\varphi}=\frac{Q}{A\cdot\varphi},$ (20)

where $v_{s}$ is the superficial velocity based on flow-rate and cross-section
of the reactor $A=0.132\times 0.132~{}\mathrm{m^{2}}$ and $\varphi$ the
porosity of the packed bed, here $\varphi=0.32$. $Q$ is the volume flow rate,
here taken to be $Q=38.2$ and $63.6~{}\mathrm{l/min}$ for Reynolds 300 and
500, respectively.

###  3.1 Experimental setup

Advanced imaging technologies, such as Radar Tomography or Magnetic Resonance
Imaging are considered in companion projects of this study. They appear very
promising, but are extremely complex and do not provide a very high resolution
in space and/or time (Tropea et al, 2016; Piro et al, 2017; Poelma, 2020). The
experimental experience of our own group is focused on optical flow
measurements. Regarding velocity and fluctuations, Particle Image Velocimetry
(PIV) is currently the gold standard. However, PIV is only possible with
sufficient optical access. For almost all published studies using optical
measurement techniques in packed beds, Refractive Index Matching (RIM) has
been used, employing a liquid with the same refraction index as the
(transparent) bed particles to enable PIV (Hassan and Dominguez-Ontiveros,
2008; Larsson et al, 2018; Wood et al, 2015) or adjusting the solid material
to the refractive index of the working fluid (Harshani et al, 2016). Hassan
and Dominguez-Ontiveros (2008) applied Particle Tracking Velocimetry (PTV) and
high-speed PIV to investigate the flow in a vertical packed bed of randomly
arranged PMMA spheres with a diameter of 4.7 mm for particle Reynolds numbers
of 80 to 500 and characterised the vortex formation at the pore level scale.
Two dimensional flow fields of laminar flows with a particle Reynolds number
in the range of 5 to 600 through a porous bed of glass spheres with a diameter
of 15 mm where obtained using PIV by Wood et al (2015). The imaging system was
used to identify the bead centre location allowing to generate the exact
geometry of the randomly assembled bed, used for direct comparison to DNS
simulations. Larsson et al (2018) applied tomographic PIV to the horizontal
flow through a thin porous bed consisting of vertical cylinders, which allowed
a visualization of three-dimensional vortex structures and a determination of
velocity magnitude for particle Reynolds numbers from 45 to 950.
However, the applications considered in the present study all involve gas
flows, and it is almost impossible to keep identical all non-dimensional
parameters characterizing momentum, heat, and mass transfer when switching to
a liquid using flow similarity. As a consequence, Refractive Index Matching
cannot be used, and alternatives must be found. In 2D porous media, innovative
methods have been introduced for visualizing interstitial pore fluid flow
using fluorescence imaging and high-speed particle tracking (Li and Iskander,
2022). However, they cannot be directly applied to real, 3D configurations. To
the best of our knowledge, no systematic PIV measurements have been carried
out for gaseous flows in packed beds yet. A corresponding procedure is
currently being developed in our group to avoid Refractive Index Matching.
After first using a simple calibration to correct resulting optical
distortions, Martins et al (2018) applied a ray-tracing based correction to
measure by PIV a jet flow behind one to three transparent spheres. A current
project of our group considers further improvements of this procedure,
enabling PIV within beds constituted of optically-transparent spheres with
gaseous flow. This leads to many challenges discussed in separate publications
(Velten et al, 2022; Ebert et al, 2022).

####  3.1.1 Description and instrumentation

The gas flow field through a packed bed with 21 layers of spheres (the choice
of this particular number will be discussed later) assembled in a body-
centered-cubic (BCC) packing is investigated by standard PIV. The main
components of the experiment, consisting of the bulk reactor mounted on a
3D-traversing unit, the laser light sheet (in green), and the PIV camera are
shown in Fig. 1(left). Liquid Di-Ethyl-Sebacat (DEHS) tracer particles
provided by a liquid nebulizer (Type AGF 10.0, Palas GmbH) enter together with
the air adjusted by a mass-flow controller (Bronkhorst D-6371; max. 500 l/min
air) through the gas inlet at the bottom (orange arrow in Fig. 1, left). The
light sheet is created by a Nd:YAG PIV laser (Quantel Q-smart Twins 850) with
a maximum energy of 380 mJ/pulse operated at a wavelength of 532 nm. In what
follows, PIV measurements will be discussed both on the freeboard, just above
the packed bed, as well as in the interstices between the particles. For the
measurement above the packing, an Imager LX 8M camera equipped with a Nikkor
35 mm f/2D Nikon lens has been used. Within the packed bed, an Imager Intense
camera equipped with a Micro-Nikkor 105 mm f/2.8D Nikon lens has been
employed. Image acquisition took place at a frequency of 2.5 Hz. Different
recording times were applied as function of the particle Reynolds number,
acquiring 300 double-frame images for $\mathrm{Re}_{P}=300$ and 1000 for
$\mathrm{Re}_{P}=500$. All PIV evaluations were carried out with the
commercial software Davis 8.4 from LaVision GmbH. A standard cross-correlation
method (multi-pass, decreasing size from 64 x 64 to 32 x 32 pixels with 50%
overlap) has been employed to calculate the two velocity components within the
laser sheet, leading to two-dimensional vector fields. The 3D traversing unit
allowed for an accurate adjustment of the measurement position without moving
any optical component.

Figure 1: (a) Scheme of the experimental setup with the camera (C) imaging the
flow within or above the bulk reactor (R) containing the BCC packed bed, the
whole system being installed on a 3D traversing unit (T) and illuminated by a
laser light sheet (L). The orange arrow shows the gas inlet. (b) Top view of
the reactor baseplate, showing with grey circles the nine bearings holding the
spheres of the first layer at their target positions. (c) Numerical model and
dimensions of the reactor used for the simulations. The packing is not shown
here for increased readability (see later Fig. 2 for a representation with BCC
packing).

####  3.1.2 Details concerning the experimental packing geometry

The flow through the reactor shown in Fig. 1(left) involves three distinct
parts: first, the inlet zone, until the baseplate; then, the BCC packing
itself; finally, an extended outlet zone with no particles but side walls, to
minimize any uncontrolled influence from the surroundings. The reactor has
been designed to yield as much as possible symmetric, stable flow conditions,
since it involves a fully symmetric geometry and constant inlet flow-rate.
Additionally, a diffuser with a combination of honeycombs and
$4~{}\mathrm{mm}$ glass spheres is used to homogenise the flow in the inlet
zone and distribute it equally over the whole cross-section of the baseplate
shown in Fig. 1(b). A pattern of 312 regularly-arranged holes (with a diameter
of $d_{h}=4~{}\mathrm{mm}$ each) controls the flow entering the packing. No
holes are placed near the vertical walls in order to reduce wall channeling.
The nine spheres building the first particle layer are fixed in proper
position by small cavities (maximal depth of $0.8$ mm), shown as gray circles
in Fig. 1(b).

The BCC packing itself consists of 139 high-precision polypropylene spheres
with a diameter of 40 mm (tolerance $150~{}\mu$m, from Ningyang Xinxin
Stainless Steel Ball Manufacture Co., Ltd.), involving a maximum of 11 so-
called ”full layers” (3 x 3 particles, touching the vertical walls of the
container) and 10 ”weak layers” in-between (only 2 x 2 particles, without
contact to the walls). The resulting height of the BCC packing consisting of
21 layers is $0.51$ m starting from the baseplate. The distance between two
opposite vertical walls is $0.132$ m (square cross-section). Including
additionally the outlet zone, the vertical walls of the container extend up to
$0.92$ m above the ground plate.

Note that, in order to avoid any micro-movement or vibration of the spheres
when activating the traversing unit, all spheres have been glued together by a
melting adhesive, creating a local bridge at the contact points between the
particles. This leads to minute deviations from the ideal BCC packing.
Considering that these slight modifications cannot be characterized with
sufficient accuracy – and have, therefore, not been taken into account in the
simulations – all comparisons shown later on will assume that the
experimentally-employed BCC packing corresponds to an unperturbed, ideal
configuration.

###  3.2 Simulation setup

####  3.2.1 Numerical domain

In the LBM simulations, a slightly simplified version of the geometry is used,
as shown in Fig. 1(c). Note that the BCC packing is not shown in this figure
for a better readability. All important features controlling the flow are kept
identical to the experimental conditions. The geometry consists of a 3D domain
with a total height of $0.97$ m, and a constant cross-section of $0.132\times
0.132$ m2 (identical to the inner dimensions of the experimentally-employed
container). The same baseplate shown in Fig. 1(b) is placed 0.14 m downstream
of the inlet plane in the simulations. The boundary condition at this inlet
plane corresponds to a fully-homogeneous flow, mimicking the action of
diffuser, honeycombs, and glass beads used to homogenize the flow in the
experiments. After the baseplate (10 mm thickness), an ideal BCC packing is
implemented (not shown in Fig. 1c), involving up to 21 layers – as in the
experiments. After reaching the top-end of the packed bed, the numerical
domain is extended by $0.31\,\rm{m}$ in streamwise direction until the outlet
boundary condition (BC). Preliminary tests have shown that this is sufficient
to avoid any undue influence of the outlet BC on the comparisons discussed
later on (though it is somewhat shorter than in the experiments). When
implementing the highest number of layers (21), the total height of the
numerical domain is $0.97\,\rm{m}$.

The inlet boundary condition is set to a constant volume flow-rate boundary
condition implemented via a modified form of the half-way bounce-back method
(Krüger et al, 2017). At the outlet, a constant-pressure boundary condition
enforced through the non-equilibrium extrapolation method (Zhao-Li et al,
2002) is used, and no-slip boundary conditions are implemented along the
vertical walls. The bounce-back method with curved treatment extension,
detailed in previous sections, is used to enforce the no-slip condition.

The dimensions of the 3D computational domain are $0.132\times 0.132\times
0.97\,\rm{m^{3}}$. After preliminary tests, the grid resolution used for all
simulations corresponds to $\delta x=5\times 10^{-4}\rm{m}$ (same in all
directions) so that each hole in the baseplate is resolved by 8 grid points.
This leads finally to a grid with $\approx 135$ Mio. grid points. The time-
step sizes are set to $\delta t=6.88\times 10^{-5}\,\rm{s}$ for
$\mathrm{Re}_{P}=300$, and $\delta t=4.13\times 10^{-5}\,\rm{s}$ for
$\mathrm{Re}_{P}=500$. This choice of time-step size guarantees that the
maximum non-dimensional velocity in the domain (corresponding to the Courant-
Friedrichs-Lewy condition, CFL) below $0.05$ guaranteeing minimal
compressibility effects in the obtained solution. Note that in the context of
the present study, we consistently use the convective CFL number based on
maximal fluid velocity, defined as:

${\rm CFL}=\frac{\|\bm{u}\|_{\rm max}\delta t}{\delta x}$ (21)

All simulations have been carried out for a total of $12\ t_{c}$ in order to
reach steady-state conditions (possibly with remaining, intrinsic
fluctuations, as discussed later on), where $t_{c}$ is the flow-through time
defined as:

$t_{c}=\frac{L_{c}}{v_{s}},$ (22)

where the characteristic speed is the superficial velocity and characteristic
length $L_{c}$ is the length of the entire flow domain in streamwise
direction. Time-averaging is only performed on data from the second half of
the simulation (6–12 $t_{c}$) to compute the average flow fields (and
associated fluctuations).

###  3.3 Locations used for comparisons and analysis

PIV measurements for $\mathrm{Re}_{P}=300$ or $500$ have been carried out for
two different configurations using the same experimental setup: Configuration
1, within selected inter-particle spaces for layer #18, 19, 20, or 21,
delivering useful information within the packed bed for locations that are
optically accessible from outside, and Configuration 2, just above the BCC
packing involving 21 layers of spheres (sometimes also called freeboard
measurements in what follows). All these measurement zones are illustrated and
described in Fig. 2. Since the corresponding information is not easy to
convey, but important for the later analysis, the reader is invited to spend
sufficient time on this figure, and come back to it later if necessary.

While all experimental results discussed in this article have been obtained
for a total number of 21 layers, the simulations have been repeated as well in
a similar procedure for 15, 17, 19, and 21 total layers in the bed. When
discussing corresponding positions in an absolute manner, the same Cartesian
coordinate system is always used, as shown in Fig. 2. The $x$ and $y$ axes
define horizontal cross-sections, with the origin at the geometrical center of
the reactor, $x$ called horizontal direction and $y$ called depth. The
$z$-axis corresponds to the vertical, streamwise direction (also called
height).

Figure 2: Front view of the packed bed reactor with 21 layers of spheres in a
BCC packing (left). A zoom of the region of interest for the later comparisons
(layers #18 to #21, and freeboard above the bed) is marked with a blue
rectangle in the left figure, and shown in the central image. In this central
view, the regions where the flow fields were measured by PIV and will be
compared to simulations are marked: the red rectangle corresponds to the
freeboard measurements above the bed (Config. 2), while the small regions
delimited by dark blue lines show the interstitial measurements (Config. 1)
within the full layers (#19 and #21), and those in light blue corresponding to
the weak layers (#18 and #20). Finally, to clarify the corresponding positions
in cross-direction (along the $y$-axis), a top view is added on the right: red
lines show measurement planes above the bed, while short blue lines correspond
to measurements within the packed bed (B1 and A2).

##  4 Comparisons with experimental data and cross-validation

The numerical results obtained with lattice Boltzmann are now compared to the
PIV data based on the two configurations described previously (within the bed,
above the bed). This procedure allows for a cross-validation, a good agreement
demonstrating the accuracy of the simulation procedure, while using at the
same time the additional details provided by LBM to better understand
experimental observations. The employed scales and domain dimensions are
always the same for experimental and numerical data, enabling a direct
comparison.

###  4.1 Configuration 1: within the bed

First comparisons concern the flow within the packed bed. Only small optically
accessible regions of around 7 x 7 mm can be used for PIV between the spheres,
as shown previously in Fig. 2. They are now considered for direct comparisons
with simulation results. Two different positions (B1 and A2, as plotted in
Fig. 2) are involved: B1, situated in the same $y$ position as the back-plane
discussed in the next section, is farther from the back wall of the container,
whereas A2 is close (only half a sphere diameter) to the front wall.

Figure 3: Comparison of the averaged in-plane velocity magnitude from LBM
(left) and PIV (right) with an overlay showing streamlines for position B1 in
layer #21 for $\mathrm{Re}_{P}=300$ (top row) and $\mathrm{Re}_{P}=500$
(bottom row). For a better visibility, only every 9th vector (PIV) and every
2nd grid value (LBM) are displayed. Figure 4: Comparison of the averaged in-
plane velocity magnitude from LBM (left) and PIV (right) with an overlay
showing streamlines for position B1 in layer #19 for $\mathrm{Re}_{P}=300$
(top row) and $\mathrm{Re}_{P}=500$ (bottom row). Please note that different
color bars are used for the two different Reynolds numbers – but are of course
the same for numerical and experimental results. For a better visibility, only
every 9th vector (PIV) and every 2nd grid value (LBM) are displayed.

Figure 3 shows the flow fields in the interstice B1 for layer #21. This is the
last inter-particle space before the flow gets released to the surface of the
bed. In Fig. 4 the same interstice, but now in layer #19 (deeper within the
bed) is presented. The same dimensions and color scales are used for LBM and
PIV, so that these visualisations allow a direct comparison between numerical
and experimental results for both Reynolds numbers, $\mathrm{Re}_{P}=300$ and
$\mathrm{Re}_{P}=500$. All plots show the in-plane average velocity magnitudes
involving the two velocity components in $x$ and $z$ directions (since the
$y$-component could not be measured by PIV), with an overlay of streamlines
showing also flow directions.

Because of the flow getting released to the free surface, the measured fields
in the top layer (#21, Fig. 3) are very different from those deeper inside the
bed (#19, Fig. 4); even the structures are completely different. While layer
#21 is dominated by a back-flow at $\mathrm{Re}_{P}=300$, or a large
recirculation zone at $\mathrm{Re}_{P}=500$, the results within layer #19
reveal a strongly horizontal flow from right to left. Since in-plane velocity
magnitudes increase only slightly with particle Reynolds numbers near the free
surface, the same color scale has been used. On the other hand, the peak
velocity becomes much higher when increasing ReP within the bed (layer #19) so
that different color scales had to be employed in Fig. 4.
For both positions and both values of $\mathrm{Re}_{P}$ a very good agreement
between numerical predictions and experimental measurements is observed. Only
very slight shifts in position are sometimes observed. The differences are a
bit larger for $\mathrm{Re}_{P}=500$, perhaps due to larger fluctuations (as
discussed later in Sec. 5.2).
The in-plane velocity magnitudes at position A2 (see again Fig. 2) situated
close to the reactor walls are now compared. Results from LBM and PIV for
$\mathrm{Re}_{P}=300$ and $\mathrm{Re}_{P}=500$ are shown in Figs. 5 and 6 for
layers #20 (closer to the bed surface) and #18 (deeper in the bed),
respectively. The position of A2 in the horizontal symmetry plane of the bed
($x=0$) should yield symmetrical flow conditions, which is overall captured by
the results, in particular for the numerical data. The flow field in the
inter-particle space at position A2 is characterized by a strong vortex pair,
which is found both in LBM and in PIV. The velocity magnitudes are always
higher within layer #20, which is the last weak layer before the flow gets
released to the surface, in comparison to the deeper layer inside the bed
(layer #18); note the different color scales used in Figs. 5 and 6. The most
pronounced difference between PIV and LBM is visible for $\mathrm{Re}_{P}=500$
in Fig. 5 (bottom row). Here, the experimental data does not show fully
symmetrical conditions and the flow structures are shifted to the upper left
and the lower right corners, revealing perhaps a minute misalignment in the
placement of the spheres for the experiments.
Overall, slightly larger differences are observed for position A2 (closer to
the container wall) compared to position B1. Close to the A2 location the
local porosity of the bed is higher because of the missing half-sphere at the
wall in each weak layer. This, in connection with very slight modifications
compared to an ideal BCC packing (as discussed in Section 4.3) might explain
the non-symmetrical flow features measured by PIV, as well as differences
compared to numerical results. Obviously, LBM simulations assume a perfect BCC
packing and only contact points between the spheres/with the walls.

Figure 5: Comparison of the averaged in-plane velocity magnitude from LBM
(left) and PIV (right) with an overlay showing streamlines for position A2 in
layer #20 for $\mathrm{Re}_{P}=300$ (top row) and $\mathrm{Re}_{P}=500$
(bottom row). Please note that different color bars are used for the two
different Reynolds numbers – but are of course the same for numerical and
experimental results. For a better visibility, only every 9th vector (PIV) and
every 2nd grid value (LBM) are displayed. Figure 6: Comparison of the averaged
in-plane velocity magnitude from LBM (left) and PIV (right) with an overlay
showing streamlines for position A2 in layer #18 for $\mathrm{Re}_{P}=300$
(top row) and $\mathrm{Re}_{P}=500$ (bottom row). Please note that different
color bars are used for the two different Reynolds numbers – but are of course
the same for numerical and experimental results. For a better visibility, only
every 9th vector (PIV) and every 2nd grid value (LBM) are displayed.

###  4.2 Configuration 2: just above the bed

The flow field downstream of the surface of the packed bed was measured within
the red rectangle shown in Fig. 2 (center). The resulting flow fields obtained
experimentally and numerically for $\mathrm{Re}_{P}=300$ and $500$ are
compared for two different measurement positions: back (Figs. 7 and Fig. 8)
and middle (Figs. 9 and Fig. 10). The corresponding positions have been
illustrated in Fig. 2. All the results shown start $10~{}\mathrm{mm}$ above
the bed. Assuming a stable flow, symmetric results are expected due to the
symmetric geometry. Indeed, symmetric flow features are mostly observed in
Figs. 7 to 10.

Figure 7: Comparison of numerical data (LBM, top) and experimental results
(PIV, bottom) for $\mathrm{Re}_{P}=300$ at the back position above the bed
concerning the averaged in-plane velocity magnitudes, with the corresponding
vector field shown as an overlay to illustrate the flow direction. For a
better visibility, only every 10th vector (PIV) and 15th grid value (LBM)
vector are displayed – the full resolution is much higher. Figure 8:
Comparison of numerical data (LBM, top) and experimental results (PIV, bottom)
for $\mathrm{Re}_{P}=500$ at the back position above the bed concerning the
averaged in-plane velocity magnitudes, with the corresponding vector field
shown as an overlay to illustrate the flow direction. For a better visibility,
only every 10th vector (PIV) and 15th grid value (LBM) vector are displayed.

The back position (Figs. 7 and 8) yields a flow field that is highly
influenced by wall channeling effects, due to the missing half-sphere close to
the wall in all weak layers. Therefore, two high-speed jets appear close to
the walls. They are found for both Reynolds numbers, and are captured
similarly by simulations and experiments. Comparing the numerical results with
the experimental ones, a very good agreement is found for
$\mathrm{Re}_{P}=300$. However, the LBM simulations overestimate slightly the
height and the width of the jets, inducing somewhat smaller recirculation
zones in the central part of the domain. The flow at $\mathrm{Re}_{P}=500$ is
again dominated by the high-speed jets near the side walls. Experiments and
simulations reveal two additional, but less pronounced oblique jets at
positions where, in this investigation plane, the flow is not blocked by the
underlying layer of spheres. The agreement between PIV and LBM is again good
at this location.

Figure 9: Comparison of numerical data (LBM, top) and experimental results
(PIV, bottom) for $\mathrm{Re}_{P}=300$ at the middle position above the bed
concerning the averaged in-plane velocity magnitudes, with the corresponding
vector field shown as an overlay to illustrate the flow direction. For a
better visibility, only every 10th vector (PIV) and 15th grid value (LBM)
vector are displayed. Figure 10: Comparison of numerical data (LBM, top) and
experimental results (PIV, bottom) for $\mathrm{Re}_{P}=500$ at the middle
position above the bed concerning the averaged in-plane velocity magnitudes,
with the corresponding vector field shown as an overlay to illustrate the flow
direction. For a better visibility, only every 10th vector (PIV) and 15th grid
value (LBM) vector are displayed.

The next comparisons involve the middle position (Figs. 9 and 10). Here also,
the dominating flow structure is well captured, mainly characterised by two
jets appearing above the free inter-particle spaces of layer #21 and leading
to two vortices close to the walls. The recirculation zone with $z$-negative
velocity direction in the center is also reproduced by the simulation, but is
noticeably stronger in the experiment at $\mathrm{Re}_{P}=300$. Overall, Fig.
9 shows the poorest comparison between simulation and experiment, which might
indicate some unexpected perturbation coming from the outflow. The
fluctuations discussed later in Sec. 5.2 might also possibly explain at least
part of this visible difference. At higher inflow velocity
($\mathrm{Re}_{P}=500$, Fig. 10) the high-speed jets dominate even more the
resulting flow fields. The position of the left jet appears somewhat shifted
toward the top-left corner in the simulation compared to PIV data, but the
agreement is much better than at $\mathrm{Re}_{P}=300$.

To provide a more quantitative comparison, line plots for both,
$\mathrm{Re}_{P}=300$ and $\mathrm{Re}_{P}=500$, of the vertical ($v$) and
horizontal ($u$) velocity components are now compared. These are the two
velocity components measured by planar PIV within the light sheet. The
comparisons are shown at both depth positions (back and middle) 1 cm above the
particle bed (Fig. 11). Additionally, the confidence intervals corresponding
to twice the standard deviation ($2\sigma$) of the PIV data used for
averaging, corresponding to a 95% confidence level, are plotted. This allows a
better quantification of the fluctuations observed experimentally from image
to image.
Based on Fig. 11 it can be concluded that most of the flow structures are
captured well. Keeping in mind that the vertical scales (i.e., the peak
velocity) differ from image to image, small dynamic ranges obviously magnify
relative differences. In Fig. 11, the observed differences correspond mostly
to a slight over- or underestimation of peak magnitude, and/or to a shift in
positions, as already observed previously. This will be part of the discussion
in the next section.

Figure 11: Comparison of experimental (blue line) and numerical data (red
line) for $\mathrm{Re}_{P}=300$ (top half) and $\mathrm{Re}_{P}=500$ (bottom
half) 1 cm above the bed. For each case the horizontal velocity $u$ (top part
of each subplot) and the vertical velocity $v$ (bottom part of each subplot)
are shown, together with the $2\sigma$ confidence interval from PIV as a
light-blue region. Note the different vertical scales.

###  4.3 Discussion of the comparisons and validation

Overall, a very good agreement has been obtained by comparing directly results
of LBM simulations and experimental PIV measurements. Taking additionally into
account the numerous previous verification and validation steps involving
ALBORZ (discussed, among others, in Hosseini et al (2020a)) the numerical
solver is, therefore, considered as being properly validated, and can now be
used for further investigations, as described in the next section. The
observed differences depend on Reynolds number, distance from the side walls,
and position within or outside of the bed, justifying a more detailed
discussion regarding the observed differences.

First, going from a particle Reynolds number of $\mathrm{Re}_{P}=300$ to
$\mathrm{Re}_{P}=500$, stronger discrepancies are observed between LBM and PIV
in most cases. As will be discussed separately in Sec. 5.2, this is possibly
due to low-frequency flow fluctuations connected to the onset of flow
instabilities. Additionally, it must be kept in mind that different averaging
times are used for PIV and LBM. The experimental data was recorded with a
frequency of 2.5 Hz for a total of 300 images (for $\mathrm{Re}_{P}=300$) or
1000 images (for $\mathrm{Re}_{P}=500$), leading respectively to an averaging
time of 120 s or even 400 s. On the other hand, the averaging time for the
numerical simulation is only 32 s due to the very high computational costs
coming with this large set-up. As a consequence, it is expected that the PIV
results show the real averaged fields, while the rather short averaging time
of LBM might suggest that the numerical results are not always fully converged
to steady-state. This would also explain why some numerical results are still
not symmetrical, while perfect symmetrical set-up and boundary conditions are
considered in the simulation.

Flow fluctuations are probably also the reason why the comparisons are better
within the bed than above it. Obviously, the flow direction is strongly
constrained within the BCC packing, enforcing specific flow structures.
Approaching the end of the particle bed (within the very last layers), and
even more above it, this is not true anymore. The resulting flow becomes far
more sensitive to intrinsic instabilities, as well as to possible
perturbations induced at the outlet of the reactor/of the computational
domain. The video showing the freeboard submitted as supplementary material to
this article shows clearly the flow fluctuations appearing after leaving the
packed bed, while the second video demonstrates that fluctuations in time are
visually not apparent within the particle bed itself.

Finally, the better agreement observed farther from the side walls is expected
to come from strong channelling effects, resulting first from the small number
of spheres in cross-section direction (only $3\times 3$ in strong layers).
Later studies will have to increase this number, and/or use 3D-printing
techniques to include the missing half-spheres along the vertical walls.
Additionally, all the geometrical uncertainties coming with the experimental
set-up may play a role to magnify this problem. While the LBM simulations
consider a perfect BCC packing with point-contact between the spheres/with the
walls, the reactor used for PIV measurements involves unavoidable
uncertainties, in particular: manufacturing tolerance of the spheres, leading
to very small variations in diameter or imperfect sphericity, insufficient
precision of construction or placement of the side walls, slightly varying
wall thickness, walls being locally non-flat (bulging toward or away from the
packing). As a consequence, small gaps could appear locally between packing
and container walls, impacting noticeably channeling effects. Finally, the
manual positioning of the spheres and the gluing process might induce slight
deviations from an ideal BCC packing. The small glue bridges between the
spheres are definitely larger than the contact points assumed in the
simulations.

Keeping all those points in mind, the overall agreement observed between LBM
simulations and PIV measurements is deemed very good, successfully terminating
the validation step regarding lattice Boltzmann simulations of the flow
through packed beds. The rest of this article deals with a deeper flow
analysis using only numerical simulations.

##  5 Further numerical studies regarding gas flows in BCC packing

###  5.1 Effect of the number of particle layers

This section discusses the impact of the number of layers on the resulting
flow structure at the end of and above the bed. Preliminary experimental
studies revealed that an odd number of layers around 20 is necessary to get
results that do not change when adding further particle layers; 17 was a
minimum, 21 layers have been finally used to have a safety margin.

To check that point by simulations, the flow conditions obtained with various
bed heights are now assessed, comparing the results obtained with 15, 17, 19,
and 21 layers. In a packed bed reactor, it is well-known from the literature
that the bed height affects pressure drop, channeling effects, flow
distribution, and residence time of the fluid. If the bed height is too small,
the reactor may not provide sufficient residence time for the desired reaction
to occur, while a bed height that is too large will lead to excessive pressure
drop and possibly poor flow distribution. Therefore, determining the minimum
number of layers required for the flow to become independent of the bed height
is crucial in ensuring optimal reactor performance.

Figure 12: Comparison of the average in-plane velocity magnitude in the
central vertical plane ($y=0$) just above the BCC packing for
$\mathrm{Re}_{P}=500$ with decreasing number of layers (from top to bottom):
21, 19, 17, 15.

To this end, four distinct simulations were performed for
$\mathrm{Re}_{P}=500$, using 15, 17, 19, or 21 layers, respectively. The
highest bed with 21 layers corresponds to the configuration used up to now in
this article. All simulations were conducted using identical boundary
conditions, the same grid spacing and time-step, in order to ensure a fair
comparison between the results. After careful analysis of the results, as
confirmed visually by looking at Fig. 12, it becomes evident that the flow
pattern does not change noticeably in the simulations when using an odd number
of layers exceeding at least 17. While clear differences are observed in the
flow above 15 or above 17 layers in the BCC packing, increasing further the
total number of layers to 19 or 21 leads only to very small changes. This
confirms the validity of our previous choice, using 21 layers for comparisons
between PIV measurements and simulation data.

In Figure 12 the average in-plane velocity magnitudes at $\mathrm{Re}_{P}=500$
in the middle plane of the reactor containing 15, 17, 19 or 21 layers are
presented. Note that, in this purely numerical study, all three velocity
components are of course known and could be used to compute the velocity
magnitude. In order to facilitate comparisons with the previous figures, a
representation of the in-plane velocity magnitude (computed from the $x$ and
$z$-components only) has been kept.

To gain a more comprehensive understanding of these changes in flow velocity,
Figure 13 provides a quantitative analysis of the in-plane velocity magnitude
within and above the last particle layer, considering again 15, 17, 19, or 21
layers in total. In Fig. 13, the results have been plotted in the center of
the reactor ($y=0$) at three different heights: 1) in the middle of the last
particle layer, at the level of sphere center (left); 2) directly after
exiting the particle bed (center); 3) 3 cm above particle bed (right figure).
No noticeable differences are observed in the left figure, still located
within the bed (middle of last particle layer). Just after leaving the
particle bed (central part of Fig. 13), the results obtained with 17, 19, and
21 layers are fully identical, while the blue curve (corresponding to 15
layers) already shows small but visible deviations. This confirms our previous
statement – at least 17 layers should be used to get exactly the same flow
conditions up to leaving the particle bed. Finally, as the flow proceeds
further in the free space above the bed, as shown in the right part of Fig.
13, the differences amplify strongly. The results with a 15-layer bed look
completely different there. Qualitatively, the velocity profiles with 17, 19,
or 21 layers look the same, but visible differences start to appear as well,
and the profiles are not fully symmetric any more. This will be the subject of
the following last section.

Figure 13: Comparison of in-plane velocity magnitudes in the center of the
reactor ($y=0$) for $\mathrm{Re}_{P}=500$ at 3 different heights when
increasing the total number of layers from 15, 17, 19, up to 21 layers: (left)
in the center plane of the last particle layer; (middle) just after leaving
the particle bed; (right) 3 cm above the particle bed. Note the different
vertical scales.

###  5.2 Velocity fluctuations

In the course of this study, the experiments and a video of the numerical
results (provided as supplementary material to this article) revealed that the
gas flow, after exiting the packed bed, shows clear fluctuations in the
freeboard region above the packing. Due to the relatively low particle
Reynolds numbers involved in this study, the corresponding frequencies are
quite low. They are, therefore, not easy to analyze accurately.
Experimentally, PIV data was acquired at a low frequency, with no attempt to
get full resolution in time. Regarding LBM simulations, the physical time that
can be solved was up to now limited to 32 s due to the high cost of these
computations.

In classical CFD, unsteady phenomena can be first identified by looking at the
evolution of the residuals. Such residuals can of course also be computed in
LBM simulations. Doing this, it was found that the residuals for velocity
magnitude drop very sharply during the first 13 s (for $\mathrm{Re}_{P}=300$)
or 8 s (for $\mathrm{Re}_{P}=500$) of physical time before reaching a low-
level, nearly constant value. This plateau, again, depends on the particle
Reynolds number. While the normalized residuals oscillate slightly around
$5\cdot 10^{-6}$ at $\mathrm{Re}_{P}=300$, they cannot drop below $4\cdot
10^{-5}$ at $\mathrm{Re}_{P}=500$, revealing larger variations of velocity
with time at the higher particle Reynolds number. The analysis of flow
fluctuations over time provides valuable insights regarding the dynamics of
the flow, and is an important step for a better understanding of fluid flows
in porous media (Boutt et al, 2007).

In the previous comparisons and discussion, it was mentioned several times
that low-frequency fluctuations have been observed in the numerical
simulations, possibly explaining non-symmetrical flow features and differences
from the experiments. For this reason, it was decided to pursue the LBM
simulations for a noticeably longer physical time by accessing a High-
Performance Computer. However, saving all the corresponding results on disk
was not possible due to storage limitations. This is why only selected results
can be discussed in what follows.

As a revealing example of the observed, low-frequency fluctuations, Figure 14
provides a visual representation of how the changes in flow structures take
place with time. This figure shows the instantaneous in-plane velocity
magnitude in the center plane of the reactor ($y=0$) within and above the last
particle layer (layer #21), in the region bounded by the red rectangle plotted
on the left side. Results separated by 6 s of physical time during somewhat
more than one minute in total (from 40 to 106 s) are shown from top left to
bottom right. A very clear jet flapping phenomenon is observed here for
$\mathrm{Re}_{P}=300$. Interestingly, jet flapping is not symmetric, the right
jet changing completely horizontal direction, while the left jet keeps its
dominating orientation towards the wall in the simulated time windows of 66
seconds. It is expected that, if it would be possible to pursue this very
costly simulation even further, symmetry would re-establish later on with a
left jet switching as well its horizontal direction. Similar low-frequency
fluctuations have been observed also for the higher particle Reynolds number.

Analyzing individual images, it appears possible that these flow movements are
the results of weak vortex shedding in the wake of the top-layer particles. An
oscillating flow pattern is created as a result of such vortices being shed
alternatively on both sides of the obstruction. Due to the interactions
between the shedding processes behind the different particles and with the
side walls symmetry is broken, inducing complex flapping movements. The exact
geometry of packing and reactor, fluid characteristics, flow rate and Reynolds
number will impact the flapping behavior (Shi et al, 2011). For all
applications where the behavior of the flow after leaving the packed bed is
important for the final process outcome, a better understanding of these low-
frequency fluctuations – or possibly finding a way how they could be
suppressed or controlled – would be important. This will be the subject of
future work.

Figure 14: Representation of the instantaneous in-plane velocity magnitude in
the central symmetry plane ($y=0$) within the last particle layer and above
the bed (see red rectangle on the left side) for a particle Reynolds number of
300. From top left to bottom right, 12 results are shown separated by 6 s of
physical time each, leading to a total duration of 66 s. Figure 15: Root Mean
Square (RMS) of velocity components $u$ (horizontal direction, left part of
each subfigure) and $v$ (vertical direction, right part of each subfigure) for
$\mathrm{Re}_{P}=300$ (left) or $\mathrm{Re}_{P}=500$ (right). The values
inside the packed bed are represented with a logarithmic color scale in order
to improve readability, while the much larger RMS values above the bed are
shown with a linear scale, but with different scales for the two different
Reynolds numbers.

Root Mean Square (RMS) is a widely used statistical measure that quantifies
the fluctuations in a signal. Here, RMS is used to analyze the magnitude of
velocity fluctuations, both within the BCC packing and in the freeboard above
it. Low RMS values indicate negligible fluctuations, while high values will be
associated with noticeable unsteady features (Kawamura et al, 2002). The
resulting values, presented in Fig. 15, have been computed separately for the
horizontal ($u$) and vertical ($v$) in-plane velocity components. Again, the
two particle Reynolds numbers $\mathrm{Re}_{P}=300$ and $\mathrm{Re}_{P}=500$
have been considered. The results reveal that the flow above the packed bed,
with RMS values plotted using a linear color scale, shows fluctuations orders
of magnitude larger than within the BCC packing, in which a logarithmic scale
must be used to be readable. This shows that, within the packed bed,
fluctuations are extremely small except in the last particle layer. They do
not increase much when going from $\mathrm{Re}_{P}=300$ to
$\mathrm{Re}_{P}=500$, but the regions showing noticeable fluctuations grow in
size.

The situation at the exit of and above the bed is completely different, with
fluctuations showing amplitudes exceeding locally 6 cm/s for
$\mathrm{Re}_{P}=500$. Looking at the different color scales, it can be seen
that the peak RMS values found above the bed are roughly three times higher
for $\mathrm{Re}_{P}=500$ than those for $\mathrm{Re}_{P}=300$. This indicates
that velocity variations with time become considerable after exiting the
particle bed. A more detailed analysis of these velocity fluctuations by
Proper Orthogonal Decomposition (Lehwald et al, 2012; Janiga, 2019; Chi et al,
2022) would be very interesting; however, POD necessitates acquisition (by
PIV) and storage (for LBM simulations) of more snapshots over longer times.

##  6 Conclusions

The focus of this study is to develop and validate a numerical model that can
accurately describe the behavior of gas flows through a packed bed reactor. To
achieve this, an existing hybrid solver combining the lattice Boltzmann
formulation with finite differences is improved further to examine the flow
behavior of a model packed bed, particularly above and in-between the last
layers of the BCC packing. For validation purposes, the results obtained from
the numerical model have been directly compared with experimental data
obtained by Particle Image Velocimetry for exactly the same configuration.
Detailed comparisons within the last 4 layers of the packed bed reveal very
good agreement with the measurement data. The remaining, small discrepancies
can be mostly attributed to minute differences regarding packing arrangement
and reactor geometry in the real set-up, compared to the ideal geometry
considered in the simulations, as well as to low-frequency fluctuations in
time. Then, the validated numerical model has been used to investigate the
effect of the number of particle layers on the flow structure at the outlet of
the bed. This analysis revealed that a minimum number of 17 layers must be
used to get always the same flow structures at the exit of the BCC packing,
confirming preliminary experimental observations.
Interestingly, the study also revealed considerable velocity fluctuations in
the freeboard just above the bed, even at low particle Reynolds numbers. At
the same time, the magnitude of the velocity fluctuations within the particle
bed is extremely low, except within the very last layer. Overall, these
findings indicate that the lattice Boltzmann simulation is an effective tool
for modeling and analyzing the fluid dynamics of such systems. To summarize:

  * •

A combined experimental (PIV) and numerical (LBM) study of the same
configuration is a particularly attractive solution for cross-validation and
to maximize our understanding of the main flow features in packed beds.

  * •

The lattice Boltzmann solver developed in this study is suitable for getting
accurate predictions of the resulting flow.

  * •

A minimum of 17 particle layers must be used to ensure that the flow
structures obtained at the outlet of the particle bed do not change any more.

  * •

Even at such low particle Reynolds numbers as 300, very noticeable velocity
fluctuations are observed at the outlet of and above the particle bed,
associated to low-frequency oscillations of the flow structure.

The results of this study have the potential to enhance the understanding of
gas dynamics in packed bed reactors. A future perspective would be to extend
the methodology by taking into account heat transfer around non-spherical
particles (Namdar et al, 2023) and chemical reactions (Hosseini et al, 2022a),
with the ultimate objective of optimizing corresponding high-temperature
conversion processes in packed beds.

## Ethics declarations

  * •

Funding:  This work has been funded by the Deutsche Forschungsgemeinschaft
(DFG, German Research Foundation) within CRC/TRR 287 “BULK-REACTION” under
number 422037413.

  * •

Conflict of interest/Competing interests:  The authors declare that they have
no conflict of interest.

  * •

Ethics approval:  Not applicable. The research did not involve any human or
animal participants.

  * •

Informed consent:  Not applicable. The research did not involve any human
participants.

  * •

Authors’ contributions:  Tanya Neeraj is in charge of all numerical
simulations and writing of the article. Chirstin Velten is in charge of all
experimental measurements and writing of the article. Gabor Janiga provided
support for post-processing and image analysis. Katharina Zähringer provided
scientific supervision of experimental work. Reza Namdar provided support for
the article correction. Fathollah Varnik provided support for the article
correction, project management. Dominique Thévenin provided support for the
article correction, project management. Seyed Ali Hosseini provided scientific
supervision of numerical work.

## Acknowledgments

This work has been funded by the Deutsche Forschungsgemeinschaft (DFG, German
Research Foundation) within CRC/TRR 287 “BULK-REACTION” under number
422037413. Regarding final simulations on High-Performance Computers, the
support of the Leibniz Supercomputing Center in Munich for accessing SuperMUC-
NG is gratefully acknowledged.

## References

  * Alobaid et al (2022) Alobaid F, Almohammed N, Massoudi F M, May J, Rößger P, Richter A, Epple B (2022) Progress in CFD simulations of fluidized beds for chemical and energy process engineering. Prog Energ Combust Sci 91:100930, DOI https://doi.org/10.1016/j.pecs.2021.1009
  * Augier et al (2010) Augier F, Idoux F, Delenne JY (2010) Numerical simulations of transfer and transport properties inside packed beds of spherical particles. Chem Eng Sci 65:1055–1064, DOI https://doi.org/10.1016/j.ces.2009.09.059
  * Bao and Meskas (2011) Bao Y, Meskas J (2011) Lattice Boltzmann Method for Fluid Simulations. Inst. Math. Sci., New York
  * Boivin et al (2021) Boivin P, Tayyab M, Zhao S (2021) Benchmarking a lattice-Boltzmann solver for reactive flows: Is the method worth the effort for combustion? Phys Fluids 33:071703
  * Boutt et al (2007) Boutt D, Jerolmack D, Hobley D, Brodu N, Roy S (2007) Direct simulation of fluid-solid mechanics in porous media using the discrete element and lattice-Boltzmann methods. J Geophys Res 112(B10)
  * Bouzidi et al (2001) Bouzidi M, Firdaouss M, Lallemand P (2001) Momentum transfer of a Boltzmann-lattice fluid with boundaries. Phys Fluids 13(11):3452–3459
  * Chapman and Cowling (1990) Chapman S, Cowling T (1990) The mathematical theory of non-uniform gases: an account of the kinetic theory of viscosity, thermal conduction and diffusion in gases. Cambridge University Press
  * Chen and Doolen (1998) Chen S, Doolen G (1998) Lattice Boltzmann method for fluid flows. Ann Rev Fluid Mech 30(1):329–364
  * Chi et al (2022) Chi C, Thévenin D, Smits A, Wolligandt S, Theisel H (2022) Identification and analysis of very-large-scale turbulent motions using multi-scale proper orthogonal decomposition. Phys Rev Fluids 7:084603/1–19
  * Dixon and Partopour (2020) Dixon A, Partopour B (2020) Computational fluid dynamics for fixed bed reactor design. Ann Rev Chem Biomol Eng 11:109–130
  * Ebert et al (2022) Ebert M, Velten C, Zähringer K, Lessig C (2022) Efficient PIV measurements in the interior of complex, transparent geometries. Conference on Modelling Fluid Flow (CMFF’22) pp 78–86
  * Eppinger et al (2014) Eppinger T, Jurtz N, Aglave R (2014) Automated workflow for spatially resolved packed bed reactors with spherical and non-spherical particles. In: Proceedings of the 10th International Conference on CFD in Oil & Gas, Metallurgical and process Industries, pp 17–19
  * Eshghinejadfard et al (2016) Eshghinejadfard A, Daróczy L, Janiga G, Thévenin D (2016) Calculation of the permeability in porous media using the lattice Boltzmann method. Int J Heat Fluid Flow 62:93–103
  * Feng et al (2015) Feng Y, Sagaut P, Tao W (2015) A three dimensional lattice model for thermal compressible flow on standard lattices. J Comput Phys 303:514–529
  * Freund et al (2003) Freund H, Zeiser T, Huber F, Klemm E, Brenner G, Durst F, Emig G (2003) Numerical simulations of single phase reacting flows in randomly packed fixed-bed reactors and experimental validation. Chem Eng Sci 58(3-6):903–910
  * Gernaey et al (2015) Gernaey K, Huusom J, Gani R (2015) Modeling the fixed-bed Fischer-Tropsch reactor in different reaction media. Comput Aid Chem Eng 37:143–148
  * Harshani et al (2016) Harshani H, Galindo-Torres S, Scheuermann A, Muhlhaus H (2016) Experimental study of porous media flow using hydro-gel beads and LED based PIV. Meas Sci Technol 28(1):015902, DOI https://doi.org/10.1088/1361-6501/28/1/015902
  * Hassan and Dominguez-Ontiveros (2008) Hassan Y, Dominguez-Ontiveros E (2008) Flow visualization in a pebble bed reactor experiment using PIV and refractive index matching techniques. Nucl Eng Design 238(11):3080–3085
  * He and Luo (1997) He X, Luo L (1997) Theory of the lattice Boltzmann method: From the Boltzmann equation to the lattice Boltzmann equation. Phys Rev E 56(6):6811–6817, DOI https://doi.org/10.1103/PhysRevE.56.6811
  * Hosseini and Karlin (2023) Hosseini S, Karlin I (2023) Lattice boltzmann for non-ideal fluids: Fundamentals and practice. arXiv preprint arXiv:230102011
  * Hosseini et al (2020a) Hosseini S, Abdelsamie A, Darabiha N, Thévenin D (2020a) Low-mach hybrid lattice boltzmann–finite differences solver for combustion in complex flows. Phys Fluids 32:077105/1–14
  * Hosseini et al (2020b) Hosseini S, Darabiha N, Thévenin D (2020b) Compressibility in lattice Boltzmann on standard stencils: effects of deviation from reference temperature. Phil Trans Royal Soc A 378(2175):20190399, DOI https://doi.org/10.1098/rsta.2019.0399
  * Hosseini et al (2021) Hosseini S, Berg P, Huang F, Roloff C, Janiga G, Thevenin D (2021) Central moments multiple relaxation time LBM for hemodynamic simulations in intracranial aneurysms: An in-vitro validation study using PIV and PC-MRI. Comput Biol Med 131:104251
  * Hosseini et al (2022a) Hosseini S, Darabiha N, Thévenin D (2022a) Low Mach number lattice Boltzmann model for turbulent combustion: flow in confined geometries. Proc Combust Inst in press
  * Hosseini (2020) Hosseini SA (2020) Development of a lattice Boltzmann-based numerical method for the simulation of reacting flows. PhD thesis, Université Paris-Saclay & Otto von Guericke University
  * Hosseini et al (2022b) Hosseini SA, Dorschner B, Karlin IV (2022b) Towards a consistent lattice boltzmann model for two-phase fluids. Journal of Fluid Mechanics 953:A4
  * Hosseini et al (2022c) Hosseini SA, Huang F, Thevenin D (2022c) Lattice Boltzmann model for simulation of flow in intracranial aneurysms considering non-newtonian effects. Phys Fluids 34(7):073105
  * Huang et al (2022) Huang F, Noel R, Berg P, Hosseini S (2022) Simulation of the FDA nozzle benchmark: A lattice Boltzmann study. Comput Meth Prog Biomed 221:106863
  * Illana Mahiques et al (2023) Illana Mahiques E, Brömmer M, Wirtz S, van Wachem B, Scherer V (2023) Simulation of reacting, moving granular assemblies of thermally thick particles by Discrete Element Method/Computational Fluid Dynamics. Chem Eng Technol in press
  * Janiga (2019) Janiga G (2019) Quantitative assessment of 4D hemodynamics in cerebral aneurysms using proper orthogonal decomposition. J Biomech 82:80–86
  * Kawamura et al (2002) Kawamura H, Abe H, Shingai K (2002) Dns of turbulence and heat transport in a channel flow with different Reynolds and Prandtl numbers and boundary conditions. In: Proceedings of the 3rd International Symposium on Turbulence, Heat and Mass Transfer
  * Krüger et al (2011) Krüger T, Varnik F, Raabe D (2011) Efficient and accurate simulations of deformable particles immersed in a fluid using a combined immersed boundary lattice oltzmann finite element method. Comput Math Applic 61:3485–3505, DOI https://doi.org/10.1016/j.camwa.2010.03.057
  * Krüger et al (2017) Krüger T, Kusumaatmaja H, Kuzmin A, Shardt O, Silva G, Viggen E (2017) The Lattice Boltzmann Method: Principles and Practice. Springer
  * Larsson et al (2018) Larsson IAS, Lundström TS, Lycksam H (2018) Tomographic PIV of flow through ordered thin porous media. Exp Fluids 59(6), DOI https://doi.org/10.1007/s00348-018-2548-6
  * Lehwald et al (2012) Lehwald A, Janiga G, Thévenin D, Zähringer K (2012) Simultaneous investigation of macro- and micro-mixing in a static mixer. Chem Eng Sci 79:8–18
  * Li and Iskander (2022) Li L, Iskander M (2022) Visualization of interstitial pore fluid flow. J Imaging 8(2):32, DOI https://doi.org/10.3390/jimaging8020032
  * Martins et al (2018) Martins F, Da Silva CC, Lessig C, Zähringer K (2018) Ray-tracing based image correction of optical distortion for PIV measurements in packed beds. J Adv Opt Photon 1(2):71–94, DOI https://doi.org/10.32604/jaop.2018.03870
  * Namdar et al (2023) Namdar R, Khodsiani M, Safari H, Neeraj T, Hosseini S, Beyrau F, Fond B, Thévenin D, Varnik F (2023) Numerical study of convective heat transfer in static arrangements of particles with arbitrary shapes: A monolithic hybrid lattice Boltzmann-finite difference-phase field solver. Particuol in press
  * Piro et al (2017) Piro M, Wassermann F, Grundmann S, Tensuda B, Kim S, Christon M, Berndt M, Nishimura M, Tropea C (2017) Fluid flow investigations within a 37 element CANDU fuel bundle supported by magnetic resonance velocimetry and computational fluid dynamics. Int J Heat Fluid Flow 66:27–42
  * Poelma (2020) Poelma C (2020) Measurement in opaque flows: a review of measurement techniques for dispersed multiphase flows. Acta Mech 231(6):2089–2111, DOI https://doi.org/10.1007/s00707-020-02683-x
  * Prasianakis and Karlin (2007) Prasianakis N, Karlin I (2007) Lattice Boltzmann method for thermal flow simulation on standard lattices. Phys Rev E 76(1):016702
  * Prokopová and Prokop (2009) Prokopová Z, Prokop R (2009) Modelling and simulation of chemical industrial reactors. In: 23rd European Conference on Modelling and Simulation, pp 378–383
  * Rong et al (2014) Rong L, Dong K, Yu A (2014) Lattice-boltzmann simulation of fluid flow through packed beds of spheres: Effect of particle size distribution. Chemical Engineering Science 116:508–523, DOI https://doi.org/10.1016/j.ces.2014.05.025, URL https://www.sciencedirect.com/science/article/pii/S0009250914002358
  * Rong et al (2015) Rong L, Zhou Z, Yu A (2015) Lattice–boltzmann simulation of fluid flow through packed beds of uniform ellipsoids. Powder Technology 285:146–156, DOI https://doi.org/10.1016/j.powtec.2015.06.047, URL https://www.sciencedirect.com/science/article/pii/S0032591015005057, pharmaceutical Particle Technology
  * Shi et al (2011) Shi L, Yu Z, Jaworski A (2011) Investigation into the strouhal numbers associated with vortex shedding from parallel-plate thermoacoustic stacks in oscillatory flow conditions. Europ J Mech B/Fluids 30(2):206–217
  * Sommerfeld et al (2008) Sommerfeld M, van Wachem B, Oliemans R (2008) Computational Fluid Dynamics of dispersed multi-phase flows. ERCOFTAC
  * Sudhakar and Das (2020) Sudhakar T, Das AK (2020) Evolution of multiphase lattice Boltzmann method: A review. J Inst Eng (India): Ser C 101(4):711–719
  * Sullivan et al (2005) Sullivan SP, Sani FM, Johns ML, Gladden LF (2005) Simulation of packed bed reactors using lattice Boltzmann methods. Chem Eng Sci 60:3405–3418, DOI https://doi.org/10.1016/j.ces.2005.01.038
  * Tariq and Liu (2022) Tariq A, Liu Z (2022) Heat transfer and friction factor correlations for slip gaseous fluid flow in confined porous medium with pore-scale LBM modelling. Int J Therm Sci 173:107382, DOI https://doi.org/10.1016/j.ijthermalsci.2021.107382
  * Tropea et al (2016) Tropea C, Yarin A, Foss J (2016) Experimental Fluid Mechanics. Springer
  * Velten et al (2022) Velten C, Ebert M, Lessig C, Zähringer K (2022) Ray tracing based reconstruction of PIV measurements in the outlet zone of gaseous flow in packed beds. In: 20th International Symposium on the Application of Laser and Imaging Techniques to Fluid Mechanics, Lisbon, Portugal, pp 51/1–11
  * Verma et al (2007) Verma N, Salem K, Mewes D (2007) Simulation of micro- and macro-transport in a packed bed of porous adsorbents by lattice Boltzmann methods. Chem Eng Sci 62(14):3685–3698, DOI https://doi.org/10.1016/j.ces.2007.04.005
  * Wan et al (2008) Wan Z, Fan F, Luo K (2008) Combined multi-direct forcing and immersed boundary method for simulating flows with moving particles. Int J Multiphase Flow 34(3):283–302
  * Wang et al (2019) Wang S, Luo K, Hu C, Lin J, Fan J (2019) CFD-DEM simulation of heat transfer in fluidized beds: Model verification, validation, and application. Chem Eng Sci 197:280–295
  * Wood et al (2015) Wood B, Apte S, Liburdy J, Ziazi R, He X, Finn J, Patil V (2015) A comparison of measured and modeled velocity fields for a laminar flow in a porous medium. Adv Water Res 85:45–63, DOI https://doi.org/10.1016/j.advwatres.2015.08.013
  * Ya-Ling et al (2019) Ya-Ling H, Qing L, Qing L, Wen-Quan T (2019) Lattice Boltzmann methods for single-phase and solid-liquid phase-change heat transfer in porous media: A review. Int J Heat Mass Transf 129:160–197, DOI https://doi.org/10.1016/j.ijheatmasstransfer.2018.08.135
  * Yang and Boek (2013) Yang J, Boek E (2013) A comparison study of multi-component lattice Boltzmann models for flow in porous media applications. Comput Math Applic 65(6):882–890
  * Young-Il (2021) Young-Il L (2021) Computational Fluid Dynamics (CFD) of Chemical Processes. MDPI
  * Zhao et al (2020) Zhao S, Farag G, Boivin P, Sagaut P (2020) Toward fully conservative hybrid lattice Boltzmann methods for compressible flows. Phys Fluids 32:126118
  * Zhao-Li et al (2002) Zhao-Li G, Chu-Guang Z, Bao-Chang S (2002) Non-equilibrium extrapolation method for velocity and pressure boundary conditions in the lattice Boltzmann method. Chin Phys 11(4):366