# Revisiting the Influence of Confinement on the Pressure Drop in Fixed Beds

Florencia Falkinhoff¹,²,³ Ⓒ · Jean-Lou Pierson¹ · Lionel Gamet¹ · Mickäel Bourgoin² · Romain Volk²

Received: 8 May 2023 / Accepted: 12 August 2023 / Published online: 31 August 2023
© The Author(s) 2023

## Abstract

We study the non-turbulent pipe flow of a Newtonian fluid through a confined porous medium made of randomly arranged spherical particles in the situation where the ratio between the pipe diameter over the particle diameter (D/d) is less than 10. Using experiments and numerical simulations, we examine the relation between the flow rate and the mean pressure gradient as a function of the Reynolds number and particle size, and how it is affected by the presence of the walls. We investigate the intrinsic variability of the measurements in relation to the randomness of the particle arrangement and how such variability is linked to spatial fluctuations of pressure within the bed. We observe that as D/d decreases, the pressure gradient presents a stronger variability, particularly in relation to where measurements are taken within the pipe. The study also quantifies the difference between measuring the pressure drop at the wall versus averaging over the entire volume, finding a small difference of 2.5% at most. We examine how the mean pressure gradient is affected by the lateral walls, finding that the pressure drop follows a consistent 1/Re scaling regardless of the confinement of the bed. Finally, we observe that the pressure gradient balances the force exerted on the solid spheres with a very weak contribution of the wall friction, showing that the role of confinement corresponds to a global effect on the bed arrangement which in turns affects the mean pressure gradient.

**Keywords** Fixed bed reactors · Porous media · Confinement · Multiphase flows · Computational fluid dynamics

## 1 Introduction

Fixed beds of particles are widely used in the chemical and process industries in a variety of reactors, like heat exchangers, separators, catalytic beds, and many other applications (Barbour et al. 2015; Elouali et al. 2019; Barker 1965). In these cases, the beds are in confined geometries, and the presence of finite reactor walls conducts to additional complexity in modelling such systems, due to supplementary effects caused by wall friction

Jean-Lou Pierson, Lionel Gamet, Mickael Bourgoin and Romain Volk have authors contributed equally to this work.

Extended author information available on the last page of the article
and local porosity variation near the reactor wall (Mueller 2019; Guo et al. 2019). These effects have been widely studied by measuring the global pressure drop in several regimes, spanning from laminar to turbulent ones, both experimentally (Erdim et al. 2015; Foumeny et al. 1993; Clavier et al. 2015; Bağcı et al. 2014) and numerically (Magnico 2003; Reddy and Joshi 2010; Dixon 2021). Nevertheless, there is not yet a clear agreement on a universal model that thoroughly describes all the regimes and possible configurations, especially when wall effects become significant, that takes into account all the parameters involved in the system.

For a fixed bed of spherical particles, the parameters involved are the Reynolds number $Re = Ud/v$ based on the superficial velocity (i.e. the velocity as if the spheres were not present) and the fluid kinematic viscosity $v$; the porosity $\epsilon$ of the bed, defined as the ratio between the volume occupied by the spheres and the total volume of the reactor, $\epsilon = V_{fluid}/V_{total}$; and $D/d$, which is the ratio between the reactor and particle diameter, $D$ and $d$, respectively. When wall effects are negligible, i.e. when $D/d \to \infty$, the pressure drop depends only on Re and $\epsilon$. In this case, the Ergun model (Ergun 1952) is generally used to describe the pressure drop of random arrangements of spherical particles for low fluid volume fraction, and new correlation has been recently derived in Dixon (2023) by considering different datasets. Nevertheless, this model is not suitable for cases where the wall effects are evident, particularly when $D/d < 10$ (Eisfeld and Schnitzlein 2001; Erdim et al. 2015; Hill et al. 2001b; Clavier et al. 2015; Foumeny et al. 1993; Flaischlen et al. 2021). It has also been reported that, in the laminar regime with walls, the Ergun model underestimates the pressure gradient because of the additional reactor wall friction in the low-Re regime, whereas in the higher-Re case the pressure drop is less than that of Ergun since the effects of the local porosity near the walls increase ($\epsilon \approx 1$), therefore acting as a less resistant path for the flow to go through (Eisfeld and Schnitzlein 2001).

In Eisfeld and Schnitzlein (2001), the variability of the coefficients of the Ergun equation was studied by analysing more than 2300 data points with different wall effects. The authors found no clear correlation between the empirical parameters of the Ergun correlation and $D/d$, especially for $D/d < 10$. Models have been proposed where there is a distinction between a bulk and a wall zone, that differentiates the flow going near the reactor walls from the flow in the middle of the bed, where the border effects would be negligible. Such an example is the model proposed by Di Felice and Gibilaro (2004) where the authors define a bulk velocity in terms of $D/d$. In De Klerk (2003) and Zou and Yu (1995) works, the bulk porosity of the bed has been modelled in terms of $D/d$ as well, by fitting various experimental data.

Most of the experiments are done by measuring the pressure drop via pressure sensors installed at the wall of the reactor, and it is worth asking if there is not a bias in the measurement that comes from measuring at some particular locations along the bed, where there can be pressure variations. Not only that, but by measuring the pressure at the wall of the reactor, we are measuring a quantity over a specific part of the experiment, but not a quantity that is averaged over the whole volume of the experiment. Are we indeed measuring the same quantity when we measure at the walls and over the whole volume of the reactor, as it is done in volume-averaging techniques? How do the variations of porosity, pressure and velocity affect the different ways of measuring? There is also an intrinsic variability associated with the randomness of the arrangements of spheres in fixed beds so that the pressure drop will undoubtedly change for each repetition of the same experiment, where the only modified parameter is the bed random packing. Even though there have been many models and studies done (see for example, Table 1 of Erdim et al. (2015)), few have been oriented towards the sources of discrepancies between the different ways of measuring the
pressure drops and if an Ergun-type correlation can be considered still valid. It is not trivial that this law would hold when there are finite border effects, given that there can be a reactor boundary layer that adds up an additional term in the momentum equation.

In order to explore these questions, pressure drops in fixed beds are studied both experimentally and numerically for three different configurations in the inertial regime. These systems are representative of industrial configurations used in thermal energy storage technologies (Sciacovelli et al. 2017; Barbour et al. 2015). The simulations allow us to go further into details and quantify data that can not be easily measured in experiments, such as the difference between averaging the pressure gradient over the whole fluid volume versus measuring it along the wall as it is done in the experiments. The simulations also provide us full access to the pressure field which is used to investigate its local variability and how it affects the estimate of the mean pressure gradient. The article is organized as follows: Sects. 2 and 3 detail the experimental and numerical methodologies used in this work. In Sect. 4, we present an analysis on the variability of the pressure field caused by the random arrangement of spheres in simulations and experiments, and quantify how it impacts the estimate of the mean pressure gradient when measuring pressure at the wall or in the bulk of the flow domain. Section 5 presents the mean pressure drop computed for three different configurations $D/d = [5, 8, 10]$, and we distinguish the contributions on the pressure drop that come from the spheres and the walls. This is done so that we can study whether the difference observed between the Ergun relation for $D/d \to \infty$ and the pressure drops measured for $D/d < 10$ comes from the presence of the walls themselves or from another phenomenon. Finally, Sect. 6 presents the summary and conclusions.

## 2 Experimental Setup

The pressure drop measurements are conducted in the experimental device shown schematically in Fig. 1. It is made of a closed upstream-flow water loop that is mainly composed of a centrifugal pump, a solenoid valve to regulate the flow, a cylindrical test section and a flow meter.

The set-up is made of a clear plexiglass column with a diameter $D = 0.04$m, and different sections of different heights that allow flexibility in the placement of the bed, so as to avoid any nozzle effects. In particular, the test section of the bed is $L_{exp} = 0.4$m high. Two additional 0.2m empty sections are added so as to separate the bed from the nozzle (Sect. 1 and 2 in Fig. 1) and a 0.05 m high honeycomb is placed so as to suppress velocity fluctuations of the incoming flow.

The test section is first filled to its full height with monodispersed stainless steel spherical beads with diameter $d = [7.938, 4.762, 3.969]$ mm and density $8g \ cm^{-3}$ from Marteau and Lemarié (France). In addition, grids are placed at the bottom and top of the bed so as to keep it fixed at any flow conditions.

The setup is then filled with filtered water with temperature control using a thermal bath. Special attention is paid to the removal of trapped bubbles in the loop because the accumulation of bubbles within the pores can lead to wrong pressure measurements. To this end, a degassing tank was added on top of the upstream water loop so that the bubbles can escape the closed water loop when reaching the free surface of the tank. Additionally, the water is degassed by heating it up to 55° C, letting it circulate for nearly 12 hours, and then cooled down to 20° C, the temperature at which the experiments are carried out. This
<figure><img src="image_1.png" /><figcaption>**Fig. 1** Experimental setup sketch. The pressure sensors are mounted along the wall of the test section, and there are two grids on the top and bottom so that the beads remain packed</figcaption></figure>

facilitates suppression of the bubbles, which was corroborated visually before starting each acquisition. In total, each experiment takes about 24 h each time the bed is modified.

A constant speed centrifugal pump is used to drive the water, controlling the flow rate with a solenoid valve. The flow rate $Q$, which can be varied in the range $Q \in [35, 165] \text{ cm}^3/\text{s}$, is measured by a flowmeter *MAG-VIEW MVM-020-QA from Bronkhorst (Netherlands)*. This provides a direct measurement of the mean velocity $U = 4Q/(\pi D^2)$ that is used to define the Reynolds number $Re = Ud/\nu$ of the flow, where $\nu = 10^{-6} \text{ m}^2 \cdot \text{s}^{-1}$ is the kinematic viscosity of water. Given the diameter of the spheres and the height of the bed, the Reynolds number can be varied in the range $Re \in [220, 1100]$ so that the flow remains laminar/inertial for all flow conditions, which was confirmed by the absence of fluctuations in pressure measurements.

Pressure along the bed is measured using an array of ten pressure sensors, flush-mounted along the column, and equidistantly distributed from each other at positions $z_i$, where $z = 0$ corresponds to the bottom of the bed. The pressure probes are high-precision piezo-resistive sensors with a sensitivity of 150mV/bar and a linear deviation of 0.05%/bar. In order to increase the signal-to-noise ratio, the signal is amplified by a homemade voltage amplifier with a gain $G = 10$ and digitized using a high precision data acquisition system (NI Dacq 4472: 8 channels, 21 bits, 100 kHz from *National Instruments*). As a consequence, signals from the flow meter and seven pressure sensors are recorded at a time with a 20 kHz sampling frequency, so that eventual pressure fluctuations due to fluidization or
turbulence can be detected. Before the first experimental campaign with flowing water, the sensitivity of the sensors was checked and calibrated for their offsets by in situ measurements of the static pressure when increasing the water level in the setup. When operating with a flow, the mean pressure $p_i$ measured at altitude $z_i$ contains both contributions from a possible offset, the static pressure, and the flow when $U \neq 0$. As all contributions are additive, the mean pressure due to the flow, further noted $p(z)$ for simplicity, is simply obtained by subtracting the mean pressure at $U = 0$ (which contains the offset and the gravitational contribution) from the current measurement. The mean pressure gradient $\frac{dp}{dz}$, is then calculated by computing the slope of the $p(z)$ curve as shown in Fig. 2 for a particular case of $D/d = 10$ and $Re = 220$ as an example. As can be noted, only sensors placed away from inlet and outlet of the bed are used for the pressure gradient measurement in order to avoid possible biases from the top and bottom boundaries.

## 3 Computational Approach

Similarly to experiments, the fixed bed generation is addressed numerically using a vertical cylindrical domain (for $D/d = 5$ and $D/d = 8, 10$, respectively) and periodic boundary conditions at the top and bottom. The random arrangements of beads in the cylindrical container are built using the Discrete Element Method (DEM) software Grains3D (Wachs et al. 2012). The fixed beds are done in three successive steps: first, the particles are dropped into a non-periodic cylindrical container with a fixed bottom. During this first step, the spheres are driven by gravitational and contact forces until they form a packed bed with a random arrangement. Afterwards, the two periodic ends of the bed are set at the minimum and maximum positions of the settled particles (these will be the inlet and the outlet of the system), and they are given an initial random velocity so that they can move inside the now periodic domain. At the same time, their radius is increased at each time iteration. During expansion, particles experience multiple collisions before they reach the final diameter that satisfies the maximum solid concentration possible, and the DEM simulation is stopped just before the particles are in contact with one another. Lastly, once we have the maximum possible radius so that there are almost no particles touching each other, a third simulation is done with a given initial random velocity such that the particles re-accommodate into their final position, which will be the one that we will use to mesh the bed. The resulting beds are characterized by calculating fluid volume fraction, $\epsilon_f$, as a function of the radial coordinate $r$ for the three cases considered here $D/d = [5.13, 8.03, 10.15]$. This is shown

<figure><img src="image_2.png" /><figcaption>**Fig. 2** Example of the pressure obtained at different heights $z$ at a fixed angular position in the reactor for $D/d = 10$ and $Re = 220$. The pressure gradient is calculated by computing the slope of a linear fit shown with the red dashed line</figcaption></figure>
in Fig. 3, where the results are consistent with those obtained by Mueller (1992); Benenati and Brosilow (1962); Goodling et al. (1983) who considered similar configurations.

Once the bed is characterized, direct numerical simulations (DNS) are performed with the finite-volume solver simpleFoam of the OpenFOAM library (see Weller et al. (1998)). For each bed, the flow domain in between the particles is meshed with the OpenFOAM unstructured mesh utility `snappyHexMesh`, using a standard meshing workflow in three steps. Firstly, the `blockMesh` utility is used to generate a fully hexahedral butterfly O-H topology background grid including the complete geometry, and defining refinement level 0. The `castellatedMesh` step of `snappyHexMesh` is then used to remove the background grid cells outside the fluid region and to refine cells on sphere surfaces (at levels 1 to 2), on intersections of spheres with the periodic boundaries (at level 2) and inside the gap regions between the spheres (with an increment of level 2). The level 1 corresponds to a division of level 0 cells by a factor 2 in all directions, and so on for higher levels. Finally, the `snap` step of `snappyHexMesh` is used to project the remaining refined cell faces on the sphere surfaces. This last step generates polyhedral cells near the walls. No boundary grid layer is used in the present study, as the fluid regime is laminar.

The meshing method presents a limitation for the minimum achievable porosity, given that there is a trade-off between resolving all the small gaps between the spheres, the number of cells necessary for the resolution of the equations and the porosity. If a higher solid fraction is required, it becomes more difficult to mesh the bed, since when the particles are in contact with one another they automatically mesh as a single object. This is an additional difficulty when analysing the data, as it becomes impossible to differentiate the two spheres, for example on quantities like the force on the spheres. With this limitation in mind, the maximum porosity that was achieved was approximately 50% for all cases, as shown in Table 2. Such a bed could be obtained experimentally by a gradual defluidization of a fluidized bed or obtained by sedimentation (Delaney et al. 2011; De Klerk 2003).

Moreover, the system was made periodic along the main flow direction (z-direction), so as to avoid any inlet/outlet boundary conditions effects and to

<figure><img src="image_3.png" /><figcaption>**Fig. 3** Fluid volume fraction (i.e. porosity) as a function of the radius of the bed. Results are compared with those obtained in Goodling et al. (1983), Mueller (1992) and Benenati and Brosilow (1962)</figcaption></figure>
compare it to the experimental measurements. The number of cells ranges between 45 and 70 million, depending on the case, and the axial periodicity is imposed by interpolating the periodic patches using the *cyclicAMI* boundary condition on all variables. In order to optimize the interpolation algorithm, we verified that no sphere was tangent to the top and bottom so as to avoid any strongly skewed cells or non-matching grids which would make the *cyclicAMI* interpolation more difficult. We also refined the edges of the intersection between the spheres and the periodic boundaries so that the borders match as best as possible. With all of this taken into account, the interpolating weight, which is a quantity that is equal to 1 for a perfect match between the periodic walls and 0 for the opposite case, is on average 0.99 for all meshes.

In the end, the whole system is composed of a z-periodic inlet at $z = z_{inlet}$ and a periodic outlet at $z = z_{outlet}$, the individual spheres that make up the bed, and the lateral walls of the reactor. The mesh is non-structured and the system can be described in cylindrical coordinates $(r, \theta, z)$, with $r \in [0, R]$, $\theta \in [0, 2\pi]$ and $z \in [0, H]$, where $R = D/2$ is the reactor's radius and $H$ its height. The geometrical parameters of the simulations are shown in Table 1. In order to verify that the volume meshed is representative and that there are no periodic spurious effects, we computed simulations with double its period $2H$. We verified that there was no significant difference between the results of the simulations with period $H$ and $2H$.

Once the mesh is complete, we solve the steady-state Navier-Stokes equations for the velocity field **u** and pressure field **p**

$$ \nabla \cdot \mathbf{u} = 0, \qquad (1) $$

$$ \rho \nabla \cdot (\mathbf{u}\mathbf{u}) = -\nabla p + \mu \nabla^2 \mathbf{u}, \qquad (2) $$

where $\rho$ is the fluid density and $\mu$ the dynamic viscosity. We make use of the steady-state solver *simpleFoam*, without any turbulent model as the flow remains laminar in the range of parameters considered here. This solver relies on the SIMPLE algorithm (Ferziger et al. 2002), which solves the pressure-velocity coupling in the incompressible Navier-Stokes equations (see Eq. 2) of the fluid by using an iterative method for the pressure $p$ and the velocity field **u**.

We use second-order Gaussian finite volume integration schemes to compute the different terms of Eq. (2). The *gaussLinear* scheme is used for the gradient operator, *gaussLinear corrected* for the laplacian schemes, and *bounded Gauss linearUpwind* is used for the divergence operator. The interpolation between the cell centres and cell faces was done through a linear interpolation scheme, which uses central differences for the interpolation.

<table><caption>**Table 1** Parameters of the simulations</caption><thead><tr><th>D/d</th><th>H/d</th><th>ε</th><th>Number of cells</th><th>ρ</th><th>U</th><th>ν</th><th>Re</th></tr></thead><tbody><tr><td>10.15</td><td>6.35</td><td>0.485</td><td>68.8M</td><td>1</td><td>0.53–7.89</td><td>1e<sup>-4</sup></td><td>20–300</td></tr><tr><td>8.03</td><td>5.02</td><td>0.490</td><td>58.9M</td><td>1</td><td>0.42–6.20</td><td>1e<sup>-4</sup></td><td>20–300</td></tr><tr><td>5.13</td><td>3.85</td><td>0.518</td><td>45.3M</td><td>1</td><td>0.26–3.83</td><td>1e<sup>-4</sup></td><td>20–300</td></tr></tbody></table>

$\epsilon = V_{fluid}/V_{total}$: fluid volume fraction. $Re = Ud/\nu$ Reynolds number based on the superficial velocity $U$, the diameter of the particles $d$ and the kinematic viscosity $\nu$
In order to impose a mean flow with periodic boundary conditions, a forcing term is added to the incompressible momentum equation such as

$$ \rho \nabla \cdot (\mathbf{u}\mathbf{u}) = -\nabla \tilde{p} + \mu \nabla^2 \mathbf{u} - f \hat{\mathbf{z}}, \quad (3) $$

where f is the forcing term imposed in $\hat{\mathbf{z}}$ direction. Under this framework, $\tilde{p}$ is a z-periodic pressure field satisfying $\tilde{p}(z_{inlet}) = \tilde{p}(z_{outlet})$ and is related to the total pressure p by the simple relation

$$ p = \tilde{p} + fz, \quad (4) $$

The equations are then solved in the entire flow domain with periodic boundary conditions in z for $\mathbf{u}$ and $\tilde{p}$ at the top and bottom, with no-slip boundary conditions on rigid boundaries (reactor walls and beads) for the velocity, and Neumann boundary conditions with zero normal gradient ($\nabla \tilde{p} \cdot \mathbf{n} = 0$) for the fluctuating pressure on solid boundaries. Finally, the full pressure gradient is then computed as $\nabla p = \nabla \tilde{p} + f$ which averages to $\langle \nabla p \rangle_V = f$ when the entire flow domain is taken into account.

We also used a residual control, down to $10^{-6}$ for both the velocity and pressure fields. A similar condition was also used in Magnico (2003). Convergence was reached within 5000 iterations, with final residuals ranging between $10^{-8}$ and $10^{-6}$. Given these criteria, it was possible to run the DNS with a range of Reynolds number $Re \in [20-300]$, which is lower than what is achieved in the experiments, and the higher Reynolds values ($Re \sim 200$) of the DNS are the lower ones from the experiments. Parameters are given in Table 1.

The numerical methodology was validated using the Lattice Boltzmann method computations performed by Hill et al. (2001b) as a benchmark case. Ordered arrangements of spheres are computed by placing a sphere in the middle of the three-periodic unit cell (in this case a simple-cubic cell), with a fixed solid volume fraction of 0.408, or equivalently a porosity of $\epsilon = 0.592$. The non-dimensional drag force exerted by the fluid on the sphere is defined as

$$ F^* = \frac{F}{6\pi\mu RU}, $$

where F is the module of the total force felt by the sphere. As shown in Fig. 4, this quantity was calculated and compared to the one obtained by Hill et al. (2001b). The relative error,

<figure><img src="image_4.png" /><figcaption>**Fig. 4** The non-dimensional drag force on a three-periodic simple cubic array of $\epsilon = 0.592$ as a function of $Re$. The results are compared with the work of Hill et. al Hill et al. (2001b). Inset: the relative error shown as a percentage</figcaption></figure>
defined in terms of the force calculated using OpenFOAM $F_{OF}$ and the one calculated by Hill $F_H$, $Error = |F_{OF} - F_{Hill}|/F_{Hill}$ is shown in the inset. We obtain a maximum relative error of 1.8%. In order to test the *cyclicAMI* condition, the sphere in the unit cell was also placed at the border of the domain and separated into two periodic parts, so as to assess that we get the same results. In that case, the error was no larger than 2%.

It is worth noting that even though the SIMPLE algorithm is not usually suitable to handle problems with high Reynolds numbers because it is a stationary solver, the effects of the non-stationary modes in the present configuration are weak. This is verified by comparing the pressure drop results on a fixed bed with $Re = 200$ and $D/d = 10$, obtained with SIMPLE to those given by the non-stationary algorithm PIMPLE. We compared the pressure gradient obtained at each time step with the one obtained from the stationary simulations. The SIMPLE results match down to a relative error of $O(10^{-4})$, thus validating the suitability of the SIMPLE algorithm.

## 4 On the Variability of Pressure Measurements in Packed Beds

When measuring the pressure gradient in experiments, sensors are usually mounted on lateral walls at known positions along the bed, while both surface and volume averages are possible in numerical simulations. We address this question of whether surface and volume averages of the pressure gradient are compatible and how the average pressure gradient varies when changing the particle arrangement. Indeed, as the spheres are randomly packed with a relatively weak scale separation between the particle size $d$ and the domain width $D$, $D/d \le 10$, variability in the particle arrangement may lead to a strong variability in pressure measurements.

Such variability can be observed in Fig. 5 which displays the normalized pressure, $p(z)/\rho U^2$, measured experimentally along the cylindrical reactor wall for 9 independent realizations of the bed, for $D/d = 10$ and a Reynolds number $Re \approx 276$. As can be observed, $p(z)$ decreases almost linearly along the bed but all curves are different for the different realizations of the bed. Not only there are variations of the pressure at a fixed position but the variations are randomly distributed along $z$ with a variability about the ensemble average estimated as $\delta p/\rho U^2 \approx 10$.

<figure><img src="image_5.png" /><figcaption>**Fig. 5** Normalized pressure $p/(\rho U^2)$ as a function of $z/d$ measured experimentally along the tube for nine different random arrangements of spheres with $D/d = 10$ and a Reynolds number $Re = 276$</figcaption></figure>
$$D/d = 10, Re = 200$$

<figure><img src="image_6.png" /><figcaption>**Fig. 6** Pressure field variations $p'$ at two different heights in the bed, for $D/d = 10$, $\epsilon = 0.518$ and $Re = 200$. Results from Direct Numerical Simulations</figcaption></figure>

<table><caption>**Table 2** Range of parameters studied experimentally and numerically</caption><thead><tr><th colspan="3">Experiments</th><th colspan="3">Simulations</th></tr><tr><th>ε</th><th>D/d</th><th>Re</th><th>ε</th><th>D/d</th><th>Re</th></tr></thead><tbody><tr><td>0.386</td><td>10.08</td><td>222–535</td><td>0.485</td><td>10.15</td><td>20–300</td></tr><tr><td>0.396</td><td>8.40</td><td>300–630</td><td>0.490</td><td>8.03</td><td>20–300</td></tr><tr><td>0.425</td><td>5.04</td><td>240–1100</td><td>0.518</td><td>5.13</td><td>20–300</td></tr></tbody></table>

As the investigation of the local pressure field can not be achieved experimentally, we turn to the DNS results, with similar scale separation $D/d = 10$ and Reynolds number $Re = 200$, to get a more detailed picture. As there is a strong pressure gradient along the bed, we define the pressure fluctuation field by subtracting the local average pressure in each plane $z = \text{cte}$. It reads:

$$p'(r, \theta, z) = p(r, \theta, z) - \frac{1}{\pi R^2} \int_0^{2\pi} \int_0^R p(r, \theta, z)r dr d\theta = p(r, \theta, z) - \langle p \rangle_{r,\theta}(z), \quad (5)$$

where $\langle \rangle_{r,\theta}$ denotes a spatial average over cylindrical coordinates $r$ and $\theta$ for $z = \text{cte}$. Moreover, $\langle \rangle_{\theta}$ denotes an average over $\theta$ for $r, z = \text{cte}$, and $\langle \rangle_{V}$ an average over the whole volume. Figure 6 shows the pressure fluctuations field $p'/\rho U^2$, at two planes at two different heights in the bed. As observed in the experiment, pressure fluctuates with typical amplitude a fraction of $\rho U^2$, although the fluid volume fraction in the DNS is about 25% larger as compared to the experiment (Table 2). We find that pressure not only fluctuates strongly in the bulk but also along the lateral boundary due to local porosity effects and solid-fluid interactions, which depend on the particle distribution in the bed. The fluctuation field strongly depends on the vertical position $z$, which may impact measurements of
the pressure gradient due to the geometry of the bed. Moreover this pressure fluctuations field will change each time the packing is changed, adding variability to the results.

To illustrate how variations in the bed impact the measurements of the pressure gradient, we first turn to the set of independent experiments done with $D/d = 10$ for which we displayed the results in Fig. 5 for $Re = 276$. For each realization of the bed, we perform a series of measurements of $p(z)$ at different Reynolds number (Table 2) from which we estimate the mean pressure gradient, $dp/dz$, at a given Reynolds number by a linear fit of $p(z)$ along the bed. As it is usually done, we define a non dimensional pressure gradient

$$G_p = -\frac{d}{\rho U^2} \frac{dp}{dz}$$

that we can plot against $1/Re$.

The results are displayed in Fig. 7 in which we note that all the measurements follow the same trend. The normalized gradient $G_p$ is nearly linear as a function of $1/Re$, each curve converging towards a constant at high $Re$ which corresponds to the inertial regime beyond the purely viscous (Darcy) regime.

As shown in Fig. 7, the precise values of the Reynolds number are not exactly the same for each realization so that estimating the variability is not straightforward. We therefore compute a linear fit of $G_p = \beta + \alpha/Re$ for each realization from which we get the ensemble average evolution $\overline{G_p} = \bar{\beta} + \bar{\alpha}/Re$, where $\bar{\alpha}$, $\bar{\beta}$ are averages of different values of $\alpha$ and $\beta$, drawn as a solid line in Fig. 7. We then obtain the standard deviation from the ensemble average,

<figure><img src="image_7.png" /><figcaption>Fig. 7 Pressure drop measured experimentally for nine different random arrangements of spheres with $D/d = 10$ (shown in different colors and symbols). The red solid line is the average linear fit $\overline{G_p} = \bar{\beta} + \bar{\alpha}/Re$ obtained by averaging the linear fits $G_p = \beta + \alpha/Re$ of the different realizations</figcaption></figure>
$$ \sigma_{G_p} = \left( \overline{\left( G_p \left( \frac{1}{Re} \right) - \overline{G_p} \left( \frac{1}{Re} \right) \right)^2} \right)^{\frac{1}{2}} $$

and find that the standard deviation averaged over all the values of the Reynolds number is $\sigma_{G_p} = 0.43$ with a quite small mean relative error $\sigma_{G_p}/\overline{G_p} = 2.58\%$.

As mentioned earlier, the pressure field not only changes for each different arrangement, but pressure exhibits spatial fluctuations within the bed (Fig. 6). We now explore how the estimated pressure gradient will depend on location by using the results from Direct Numerical Simulations. As measurements are usually performed on a vertical line along the bed, we first study how pressure varies along z, at the wall of the reactor, for different values of the angle $\theta$. The pressure lines at four different angles are shown in Fig. 8a for a particular case of $D/d = 10$ at $Re = 200$. The pressure lines show that the pressure as a function of the height of the bed varies depending on the angle. The pressure gradient for different angles $\langle G_p(r = R, \theta) \rangle_z$ is shown in Fig. 8b for 4 different values of the Reynolds numbers. This was calculated for different angles $\theta$ by a linear fit of the curves as is done experimentally (see Fig. 2). This quantity is normalized by the pressure gradient averaged over the whole volume, $\langle G_p \rangle_V$, which is equal to the forcing term in the DNS, $\langle G_p \rangle_V = f$. In the present case the different curves nearly collapse onto the same curve for all Reynolds numbers, which is due to the fact that the 4 simulations were performed using the same bed, so that the spatial organization of the pressure field does not change much when changing $Re$ as it is strongly correlated to the structure of the bed.

The relative variations of the pressure gradient at the walls for $D/d = 5, 8$ and $10$ are shown in Fig. 9, noted as $G'_p(r = R)$. This was calculated by computing the standard deviation of the superficial pressure drop described above,

$$ G'_{p}(r=R) = \left( \langle G_{p}(r=R) \rangle_{z}^{2} - \langle G_{p}(r=R) \rangle_{z,\theta} \right)^{1/2} $$

<figure><img src="image_8.png" /><figcaption>Fig. 8 a Pressure versus $z/d$ at the wall of the reactor ($r = R = D/2$) at different angles $\theta$, for a fixed bed with $D/d = 10, Re = 200$. b Pressure gradient at the wall of the reactor versus $\theta$, normalized by the perimeter of the bed, for the same random arrangement as in a. Results from DNS</figcaption></figure>
<figure><img src="image_9.png" /><figcaption>**Fig. 9** Standard deviations of the mean pressure gradient at different angles of the bed, $G'_p(r=R)$ compared to its mean over all angles, $\langle G_p(r=R) \rangle_{z,\theta}$</figcaption></figure>

The value of the standard deviation for $D/d = 10$ is 0.38, which is comparable to the experimental repeatability result which we recall to be 0.42. It is worth noting that the experimental repeatability result was averaged over a higher $Re$ range. This means that a direct analogy between measuring at different angles of the bed and changing the bed several times can be made, with an overestimation of the variability by using the former. Furthermore, the relative deviations increase for the smaller scale separations: it is on average 4%, 7% and 10% for $D/d = 10, 8$ and $5$, respectively. Assuming that the tendency is right, this means that the relative variability can be of order 10 for $D/d < 5$.

Finally, let us now compare the mean pressure gradient measured at the walls $\langle G_p(r=R, \theta, z) \rangle_{\theta}$ with the volumetric average $\langle G_p \rangle_V$ as a function of $Re$. This is of interest because the pressure drop is typically measured at the walls in experiments, whereas macroscopic laws are usually derived by volume averages over the whole domain. Figure 10 displays how the ratio $\langle G_p(r=R, \theta, z) \rangle_{\theta} / \langle G_p \rangle_V$ changes with $Re$ for the different scale separations investigated. It can be observed that the value does not change significantly with $Re$. It does not stray that much from unity, although there is a clear effect that depends on $D/d$, with an overestimation of the pressure gradient at the walls when the scale separations are smaller $D/d = 5$, and an underestimation for larger $D/d$. Nevertheless, this variation is found of the order of 1%, which is less than the error that

<figure><img src="image_10.png" /><figcaption>**Fig. 10** Pressure gradient calculated at the walls averaged over all $\theta$ and $z$ compared to the pressure gradient averaged over the whole volume $\langle G_p \rangle_V$</figcaption></figure>
one commits when measuring the pressure drop experimentally for different beds, or at a particular angle, as quantified in Fig. 9.

We shall therefore conclude that the error made by computing the pressure gradient at the wall instead of using a volumetric average is by far smaller than the typical variability one gets in the pressure gradient measurement when using only one realization of the bed which can reach up to 10% when scale separation is too small.

# 5 Mean Pressure Gradient

We now turn to the interpretation of the results obtained experimentally and numerically for the pressure gradient as a function of $Re$ for the three different beds with $D/d = 5, 8$ and $10$ (Table 2).

## 5.1 Dimensional Analysis

In order to show a robust relation between the pressure gradient in the bed, which is proportional to the total pressure loss $\Delta p$ along the bed, we turn to dimensional analysis. The pressure loss is a function of 7 physical parameters so that 8 dimensional quantities are involved:

$$ \Delta p, \rho, U, \mu, d, D, L, V_s, $$

where $L$ is the relevant length in the tube ($L = H$ in the simulation while $L$ is the maximal distance between the sensors in the experiment so that the pressure gradient is $dp/dz = \Delta p/L$ in all cases), and $V_s$ is the volume occupied by all the spheres. As there are 3 physical dimensions involved (mass, time and length) and 8 dimensional quantities, the functional relation between the pressure loss and the other quantities can be reduced to a relation between 5 independent non-dimensional numbers that we choose to be

$$ G_p \equiv \frac{d\Delta p/L}{\rho U^2}, Re, \varepsilon, D/d, L/D, $$

that can be written as:

$$ G_p = f(Re, \varepsilon, D/d, L/D). \qquad (6) $$

One may wonder if $L/D$ is really involved in this relation. The answer comes from an analysis of the system. If $L$ is long enough so that pressure fluctuations decorrelate along the bed, then inlet-outlet condition do not matter and the bed can be considered as infinite. In such condition, $L/D$ can be dropped in the above expression to get:

$$ G_p = f(Re, \varepsilon, D/d). \qquad (7) $$

Note that the condition $H/D \gg 1$ is easily satisfied in the experiment as we have $L/D \sim 10$ with a linear pressure loss along the bed (Fig. 5). We verified that $L = H$ was not involved as well in the direct numerical simulations by checking that $G_p$ remains the same when doubling the height of the bed.

The functional form of Eq. 7 may be complex but it takes a simple form in limit cases. Let us for instance consider a pipe flow, for which $\varepsilon$ and $D/d$ are not involved. In such
case, the pressure loss is 1/Re in the Poiseuille regime (Idelchik 1987). Back to the porous medium in this regime of vanishingly small Re, the loads are still linearly related to the mean velocity so that $G_p$ must be proportional to 1/Re so that we write

$$G_p = \alpha(D/d, \varepsilon) \frac{1}{Re} \quad (Re \ll 1),$$

where $\alpha(D/d, \varepsilon)$ is an unknown function.

On the other hand, in the high Reynolds number regime, the loads are proportional to $U^2$ so that $G_p$ will remain constant (Idelchik 1987), which allows us to write

$$G_p = \beta(D/d, \varepsilon) \quad (Re \gg 1).$$

Guessing about the evolution $G_p$ as a function of Re in the intermediate regime, which does not correspond to the viscous nor the inertial regime, is more difficult. However, one may guess that some inertial contribution may first appear as a correction of the viscous term when increasing Re so we choose to explore how $G_p$ varies as a function of 1/Re and write $G_p = \alpha/Re + \beta$. This is consistent with the work carried out in Hill et al. (2001a), where they observed that the inertial contribution is a constant with our current normalization.

There is no a priori reason for which $\alpha$ and $\beta$ should not depend on the Reynolds number in the inertial regime (Reynolds number of order 100) as there could be boundary layer or transitional effects. Also, $\alpha$ and $\beta$ may depend of $D/d$ for wall bounded porous medium. In the particular case where $D/d \to \infty$, the Ergun equation is used to describe the pressure drop on non-confined fixed beds of spherical particles (Ergun 1952), which reads:

$$G_p = \frac{(1 - \varepsilon)}{\varepsilon^3} \left( A \frac{1 - \varepsilon}{Re} + B \right). \quad (8)$$

The equation depends on $\varepsilon$ and Re, and A and B are empirical constants that do not depend on Re and are reported to be 150 and 1.75, respectively (Ergun 1952). This equation will now serve as a basis to interpret the present results which have been obtained with weaker scale separations $D/d < 10$.

## 5.2 Results

Figure 11 displays the evolution of the normalized pressure gradient $G_p = (\Delta p/L)/(\rho U^2/d)$ as a function of 1/Re both from experiments and numerical simulations (see Table 2 for the range of Reynolds number for each run). The pressure gradient measured experimentally corresponds to the one measured on the wall of the reactor, at a fixed angle $\theta_0$: $G_p(R, \theta_0, z)$, and the numerical one is presented as the gradient averaged over the wall surface, that is, the average of the gradient shown in Fig. 8b, $\langle G_p(r = R, \theta, z) \rangle_{\theta,z}$. The error bars for both results correspond to the variability at the walls measured numerically, which increase with the confinement effects.

All the experimental results seem to follow a linear trend of $G_p$ with 1/Re. The fact that $G_p = \alpha/Re + \beta$ is not a trivial result as there is no reason why the result obtained at moderate Reynolds number should be a superposition of the limit cases Re ~ 0 and Re ~ $\infty$. It is worth noting that the DNS results do not follow the same linear trend, as the Reynolds numbers are lower than the ones achieved in the experiments and we do not expect the same inertial contribution as the one observed in the experiments. Nevertheless, it is visible in this figure that, both in the experiment and the DNS, the non dimensional pressure
<figure><img src="image_11.png" /><figcaption>**Fig. 11** Non-dimensional pressure gradients $G_p = (\Delta p/L)/(\rho U^2/d)$ as a function of $1/Re$. Left: Experimental results, $\epsilon \approx 0.4$. Right: Numerical results, with $\epsilon \approx 0.5$. The error bars for both results correspond to the variability at the walls measured numerically, which increase with the confinement effects</figcaption></figure>

gradient is a decreasing function of $D/d$ in the range of parameters investigated which is a clear signature of the presence of the wall. This is probably due to the fact that as $D/d$ lowers the porosity of the bed increases (De Klerk 2003), so then the bed offers less resistance to the flow.

The slope and intercept of the linear fits are shown in Table 3, multiplied by a factor depending on the porosity so as to have a direct comparison with A and B from Eq. (8), noted as $A_{ergun}$ and $B_{ergun}$ in the table.

Both the numerical and experimental values are within the values reported in Eisfeld and Schnitzlein (2001). The overall trend is that the value of A increases with decreasing $D/d$, and they report that B shows a "slight opposite trend", although there is not a clear correlation in the results. It is worth noting that in all the experiments studied, the contributions of $\epsilon$ and $D/d$ were not separated, and the porosities of the experiments are not specified. In fact, as it has been modelled in the previous works (De Klerk 2003; Zou and Yu 1995) both variables depend on the other, so their contributions might not be separable. In other words, one might not be able to write, e.g. $\alpha(\epsilon, D/d) = \alpha_{\epsilon}(\epsilon) \times \alpha_{D/d}(D/d)$. Moreover, there are no error bars reported in the results, which we have proved to be important due to the variability, especially for lower $D/d$.

We now compare in Fig. 12 our experimental results for $D/d = 5$, which correspond to the smallest scale separation investigated, to those from other works (Foumeny et al. 1993; Erdim et al. 2015), with the Ergun model (Eq. (8) and the recent correlation obtained by Dixon (2023), two different models for the $D/d \to \infty$ limit. As can be

<table><caption>**Table 3** Experimental and numerical values of the slope and intersection of the $G_p(1/Re)$ linear fits, $\alpha$ and $\beta$ respectively</caption><thead><tr><th>D/d</th><th>A<sub>Ergun</sub></th><th>ε<sup>3</sup><sub>(1-ε)</sub> α<sub>exp</sub></th><th>ε<sup>3</sup><sub>(1-ε)</sub> α<sub>dns</sub></th><th>B<sub>Ergun</sub></th><th>ε<sup>3</sup><sub>(1-ε)</sub> β<sub>exp</sub></th><th>ε<sup>3</sup><sub>(1-ε)</sub> β<sub>dns</sub></th></tr></thead><tbody><tr><td>5</td><td>150</td><td>141</td><td>139</td><td>1.75</td><td>1.05</td><td>1.21</td></tr><tr><td>8</td><td>150</td><td>178</td><td>125</td><td>1.75</td><td>1.11</td><td>1.15</td></tr><tr><td>10</td><td>150</td><td>219</td><td>130</td><td>1.75</td><td>1.00</td><td>1.11</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**Fig. 12** Pressure drop for cases with $D/d = [5, 5.04]$, compared with the Ergun and Dixon models for ($D/d \to \infty$) and results obtained by Erdim et al. (2015). The error bars correspond to the estimated variability of the experimental measurements</figcaption></figure>

observed, the experimental curves follow the same type of $G_p$ linear relation (within the variability errorbars) as a function of $1/Re$, with good agreement between the present data and those of Erdim et al. (2015). Besides, all experimental sets show a clear deviation from the Ergun relation (red dashed line), which overestimates by a factor 3 the pressure loss along the bed. By contrast, the Dixon relation is much closer to the experimental data, especially at lower $1/Re$ (higher Reynolds), and it takes account the slight curvature of the results within the variability error.

In both cases the deviations can be attributed to the presence of the walls (Foumeny et al. 1993; Erdim et al. 2015; Mehta and Hawley 1969; Eisfeld and Schnitzlein 2001; Di Felice and Gibilaro 2004). Whether this is due to some friction at the outer boundary, or to a more global effect on the particle distribution due to the confinement, is an open question. In order to answer this question, we will turn to numerical results to get a more local analysis, which is that of the forces involved in the spheres. This will allow us to separate the contribution from the walls and from the spheres. This is presented in the following section.

## 5.3 Force Contributions

The DNS data allow for the direct computation of the forces exerted on the spheres by numerical integration of the loads $\mathbf{F} = \int_S \boldsymbol{\sigma} \cdot \mathbf{n} \, dS$ where $\boldsymbol{\sigma}$ is the stress tensor and $S$ is the surface of the particle.
We calculate the forces for all the spheres as the sum of all individual forces, and the forces over the walls of the reactor. Because of conservation of momentum, it is straightforward that

$$ \left\langle \frac{dp}{dz} \right\rangle = \frac{1}{V} \left( \sum F_{\text{Spheres}} + F_{\text{Walls}} \right), \qquad (9) $$

where $F$ represents the z-component of the force. This means that the pressure gradient is a result of the different contributions of the forces of the system. We can separate the two contributions in order to see which one is more influential. Equation (9) was verified, down to a 1% difference for all simulations.

Figure 13a shows the total force per unit of volume exerted on the bed and the wall $F_{\text{Total}}/V$ (labelled as Total), the sum of the forces over all the spheres only, labelled as Spheres, and the forces at the walls, labelled as Walls, all normalized by $\rho U^2/d$, which is the same normalization factor used in $G_p$.

As was mentioned, the total force per unit of volume is the pressure drop, and by separating the contributions we can see which is more important. It is evident that most of the contribution comes from the spheres and not from the walls. In particular, for $D/d = 5$ and $Re = 200$, $F_{\text{Spheres}} = 22 \times F_{\text{Walls}}$, and $F_{\text{Total}} = 1.05 \times F_{\text{Spheres}}$. Not only that, but both $F_{\text{Walls}}$ and $F_{\text{Spheres}}$ follow the same $1/Re$ scaling, which validates the robustness of the scaling (the walls contribution was re-plotted in the inset for clarity), and that there is not an additional term due to the walls in the ranges of parameters studied, thereby revalidating that $\alpha$ and $\beta$ only depend on $\epsilon$ and $D/d$.

Figure 13b shows $F_{\text{Total}}$ and $F_{\text{Spheres}}$ for $D/d = 10, 8$ and $5$. All of them follow the same scaling and once again most of the contribution for the pressure gradient comes from the spheres, not the walls. It is worth noting though that the wall contribution increases with $D/d$. That is, there is a dependency with the confinement, however small (as it was said, for $D/d = 5$ the spheres account for 95% of the contribution).

<figure><img src="image_13.png" /><figcaption>Fig. 13 a Force contributions on the pressure drop for $D/d = 5$, where most of it comes from the spheres. Inset: The force at the walls is enlarged for clarification. b Total force of the system and total force over the spheres for $D/d = 5, 8$ and $10$. We observe that for all the cases considered most of the contribution comes from the spheres arrangements and not from the walls themselves</figcaption></figure>
<table><caption>**Table 4** Forces per surface area in units of N/m² for D/d = 10 and 5</caption><thead><tr><th colspan="5">Forces per surface area [N/m²]</th></tr></thead><tbody><tr><td colspan="5">D/d = 10</td></tr><tr><td>Re</td><td>200</td><td>150</td><td>100</td><td>60</td></tr><tr><td>Wall</td><td>6.45</td><td>4.32</td><td>2.51</td><td>1.31</td></tr><tr><td>Spheres</td><td>69.48</td><td>44.21</td><td>23.74</td><td>11.17</td></tr><tr><td colspan="5">D/d = 10</td></tr><tr><td>Re</td><td>200</td><td>150</td><td>100</td><td>60</td></tr><tr><td>Wall</td><td>1.41</td><td>1</td><td>0.59</td><td>0.31</td></tr><tr><td>Spheres</td><td>14.52</td><td>9.29</td><td>5.05</td><td>2.43</td></tr></tbody></table>

Finally, we calculated the forces per surface area for D/d = 10 and 5 (shown in Table 4), and found that the force exerted by the walls is approximately 10% of that exerted by the spheres. Its contribution increases with increasing Re as well, which is consistent with what is shown in Fig. 13a and b. This validates the fact that the wall contribution is not as significant as the one provided by the spheres when taking into account the force per unit area, as is the case of the force per unit volume (which is representative of the pressure drop). Nevertheless, it is worth noting that the ratio $F_{spheres}/F_{walls}$ does not change as significantly with D/d as the force per unit volume.

All of this leads to the conclusion that the differences observed between the pressure drops measured in confined beds with D/d < 10 and non-confined beds (i.e. those that follow the Ergun equation) are not a result from the walls themselves, but from how the spheres are arranged because of the presence of the walls. In other words, the change in geometry imposed into the system by a reactor with D/d < 10 is what causes a difference in the forces felt by the spheres which is by itself generated by the differences in the mean porosity of the bed generated by the walls. For D/d > 10, we would expect the wall effects to become fully negligible with its force contribution being nearly zero, and to the porosity of the bed to tend towards the porosity of random beds of spheres without any border effects, which is usually $\epsilon \approx 0.4$. This would be the limit of non-confined beds, where we cite a recent review by Dixon (2023).

## 6 Conclusions

Fixed beds of randomly arranged particles present an intrinsic variability that is linked to the arrangements themselves. This variability is reflected in the pressure field inside the bed and thus the pressure gradient might vary as well. We have studied this with the aid of numerical simulations and experiments, where we have observed that, as D/d decreases, the pressure gradient presents a stronger variability depending on where we measure it at the wall of the reactor. Moreover, we were able to quantify the difference between measuring the pressure drop at the wall of the reactor and the one averaged over the whole volume. We observed that there is a small difference between measuring at the wall of the reactor and averaging over the whole volume, of the order 2.5% at the most. This is of particular interest because the theoretical developments are usually derived for bulk-averaged quantities, while pressure measurements are usually performed experimentally at
the lateral boundaries of the domain. Our study showed that the variability corresponds to an intrinsic error that can reach 10% when *D/d* is small.

We also studied how the mean pressure gradient is affected when the walls are present. As a first observation, all pressure drops follow the same 1/Re scaling, which is consistent with the Ergun correlation. Furthermore, all pressure drops were found to increase with scale separation due to a decrease in the average bed porosity, a decreasing function of *D/d*. We quantified the contribution of the walls and the solid spheres by doing a simple force balance of the system. The balance of forces satisfies that the pressure drop of the system is equal to the sum of the fluid-wall and fluid-spheres forces. We find that most of the contribution (between 95 and 98%) comes from the spheres, whereas the walls themselves have a small effect on the pressure gradient. This means that the difference documented between the pressure drop in confined beds and the infinite case described by the Ergun equation comes from how the spheres are rearranged and distributed at the presence of the walls rather than the presence of the walls themselves. The wall forces contribution to the pressure gradient is within the variability error of the measurement.

**Funding** Open Access funding enabled and organized by Projekt DEAL.

**Declarations**

**Conflicts of interest** The authors have no relevant financial or non-financial interests to disclose.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.

**References**

Bağcı, Ö., Dukhan, N., Özdemir, M.: Flow regimes in packed beds of spheres from pre-darcy to turbulent. Trans. Porous Med. **104**(3), 501–520 (2014). https://doi.org/10.1007/s11242-014-0345-0

Barbour, E., Mignard, D., Ding, Y., et al.: Adiabatic compressed air energy storage with packed bed thermal energy storage. Appl. Energy **155**, 804–815 (2015). https://doi.org/10.1016/j.apenergy.2015.06.019

Barker, J.J.: Heat transfer in packed beds. Indust. Eng. Chem. **57**(4), 43–51 (1965). https://doi.org/10.1021/ie50664a008

Benenati, R., Brosilow, C.: Void fraction distribution in beds of spheres. AIChE J. **8**(3), 359–361 (1962)

Clavier, R., Chikhi, N., Fichot, F., et al.: Experimental investigation on single-phase pressure losses in nuclear debris beds: identification of flow regimes and effective diameter. Nuclear Eng. Des. **292**(July), 222–236 (2015). https://doi.org/10.1016/j.nucengdes.2015.07.003

De Klerk, A.: Voidage variation in packed beds at small column to particle diameter ratio. AIChE J. **49**(8), 2022–2029 (2003). https://doi.org/10.1002/aic.690490812

Delaney, G.W., Hilton, J.E., Cleary, P.W.: Defining random loose packing for nonspherical grains. Phys. Rev. E Statist. Nonli. Soft Matter Phys. **83**(5), 1–4 (2011). https://doi.org/10.1103/PhysRevE.83.051305

Di Felice, R., Gibilaro, L.: Wall effects for the pressure drop in fixed beds. Chem. Eng. Sci. **59**(14), 3037–3040 (2004) https://doi.org/10.1016/j.ces.2004.03.030, www.sciencedirect.com/science/article/pii/S0009250904002714
Dixon, A.G.: General correlation for pressure drop through randomly-packed beds of spheres with negligible wall effects. AIChE J. n/a(n/a):e18035. (2023) https://doi.org/10.1002/aic.18035

Dixon, A.G.: Particle-resolved cfd simulation of fixed bed pressure drop at moderate to high reynolds number. Powder Technol. 385, 69–82 (2021) https://doi.org/10.1016/j.powtec.2021.02.052, www.sciencedirect.com/science/article/pii/S0032591021001625

Eisfeld, B., Schnitzlein, K.: The influence of confining walls on the pressure drop in packed beds. Chem. Eng. Sci. 56(14), 4321–4329 (2001). https://doi.org/10.1016/S0009-2509(00)00533-9

Elouali, A., Kousksou, T., El Rhafiki, T., et al.: Physical models for packed bed: Sensible heat storage systems. Journal of Energy Storage 23, 69–78 (2019) https://doi.org/10.1016/j.est.2019.03.004, www.sciencedirect.com/science/article/pii/S2352152X19300179

Erdim, E., Akgiray, Ö., Demir, I.: A revisit of pressure drop-flow rate correlations for packed beds of spheres. Powder Technol. 283, 488–504 (2015). https://doi.org/10.1016/j.powtec.2015.06.017

Ergun, S.: Fluid Flow Through Packed Columns. Contribution (Carnegie Institute of Technology, Coal Research Laboratory), (1952) https://books.google.de/books?id=37sgywEACAAJ

Ferziger, J.H., Perić, M., Street, R.L.: Computational Methods for Fluid Dynamics, vol. 3. Springer (2002)

Flaischlen, S., Kutscherauer, M., Wehinger, G.D.: Local structure effects on pressure drop in slender fixed beds of spheres. Chem. Ingen. Tech. 93(1–2), 273–281 (2021). https://doi.org/10.1002/cite.202000171

Foumeny, E., Benyahia, F., Castro, J., et al.: Correlations of pressure drop in packed beds taking into account the effect of confining wall. Int. J. Heat Mass Trans. 36(2), 536–540 (1993). https://doi.org/10.1016/0017-9310(93)80028-S

Goodling, J.S., Vachon, R.I., Stelpflug, W.S., et al.: Radial porosity distribution in cylindrical beds packed with spheres. Powder Technol. 35(1), 23–29 (1983). https://doi.org/10.1016/0032-5910(83)85022-0

Guo, Z., Sun, Z., Zhang, N., et al.: Mean porosity variations in packed bed of monosized spheres with small tube-to-particle diameter ratios. Powder Technol. 354, 842–853 (2019). https://doi.org/10.1016/j.powtec.2019.07.001

Hill, R.J., Koch, D.L., Ladd, A.J.: The first effects of fluid inertia on flows in ordered and random arrays of spheres. J. Fluid Mechan. 448, 213–241 (2001). https://doi.org/10.1017/s0022112001005948

Hill, R.J., Koch, D.L., Ladd, A.J.C.: Moderate-Reynolds-number flows in ordered and random arrays of spheres. J. Fluid Mechan 448, 243–278 (2001). https://doi.org/10.1017/s0022112001005936

Idelchik, I.: Handbook of hydraulic resistance, 2nd edition. Journal of Pressure Vessel Technology-transactions of The Asme - J PRESSURE VESSEL TECHNOL 109. (1987) https://doi.org/10.1115/1.3264907

Magnico, P.: Hydrodynamic and transport properties of packed beds in small tube-to-sphere diameter ratio: pore scale simulation using an eulerian and a lagrangian approach. Chem. Eng. Sci. 58(22), 5005–5024 (2003). https://doi.org/10.1016/S0009-2509(03)00282-3

Mehta, D., Hawley, M.C.: Wall effect in packed columns. Indust. Eng. Chem. Process Des. Develop. 8(2), 280–282 (1969). https://doi.org/10.1021/i260030a021

Mueller, G.E.: Radial void fraction distributions in randomly packed fixed beds of uniformly sized spheres in cylindrical containers. Powder Technol. 72(3), 269–275 (1992). https://doi.org/10.1016/0032-5910(92)80045-X

Mueller, G.E.: A modified packed bed radial porosity correlation. Powder Technol. 342, 607–612 (2019). https://doi.org/10.1016/j.powtec.2018.10.030

Reddy, R.K., Joshi, J.B.: Cfd modeling of pressure drop and drag coefficient in fixed beds: Wall effects. Particuology 8(1):37–43. (2010) https://doi.org/10.1016/j.partic.2009.04.010, https://www.sciencedirect.com/science/article/pii/S1674200109001254, particulate Flows and Reaction Engineering

Sciacovelli, A., Vecchi, A., Ding, Y.: Liquid air energy storage (laes) with packed bed cold thermal storage - from component to system level performance through dynamic modelling. Applied Energy 190, 84–98 (2017) https://doi.org/10.1016/j.apenergy.2016.12.118, www.sciencedirect.com/science/article/pii/S0306261916319018

Wachs, A., Girolami, L., Vinay, G., et al.: Grains3d, a flexible dem approach for particles of arbitrary convex shape – part i: Numerical model and validations. Powder Technology 224, 374–389 (2012) https://doi.org/10.1016/j.powtec.2012.03.023, www.sciencedirect.com/science/article/pii/S003259101200191X

Weller, H., Tabor, G., Jasak, H., et al.: A tensorial approach to computational continuum mechanics using object orientated techniques. Comput. Phys. 12, 620–631 (1998). https://doi.org/10.1063/1.168744

Zou, R.P., Yu, A.B.: The packing of spheres in a cylindrical container: the thickness effect. Chem. Eng. Sci. 50(9), 1504–1507 (1995). https://doi.org/10.1016/0009-2509(94)00483-8

**Publisher's Note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Authors and Affiliations

Florencia Falkinhoff<fnref n="1" />,<fnref n="2" />,<fnref n="3" /> Ⓒ · Jean-Lou Pierson<fnref n="1" /> · Lionel Gamet<fnref n="1" /> · Mickäel Bourgoin<fnref n="2" /> · Romain Volk<fnref n="2" />

<figure><img src="image_14.png" /></figure>

Florencia Falkinhoff
florencia.falkinhoff@ds.mpg.de

<figure><img src="image_15.png" /></figure>

Jean-Lou Pierson
jean-lou.pierson@ifpen.fr

Lionel Gamet
lionel.gamet@ifpen.fr

Mickäel Bourgoin
mickael.bourgoin@ens-lyon.fr

Romain Volk
romain.volk@ens-lyon.fr

<footnote n="1">¹ IFP Energies Nouvelles, 69360 Solaize, France</footnote>

<footnote n="2">² ENS de Lyon, CNRS, Laboratoire de Physique, UMR 5672, F-69342 Lyon, France</footnote>

<footnote n="3">³ Max Planck Institute for Dynamics and Self-Organization, Göttingen, Germany</footnote>