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

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%\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$2.5\%$$\end{document} 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.


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, Mick¨ael 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∕ based on the superficial velocity (i.e. the velocity as if the spheres were not present) and the fluid kinematic viscosity ; the porosity of the bed, defined as the ratio between the volume occupied by the spheres and the total volume of the reactor, = 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 → ∞ , the pressure drop depends only on Re and .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 ( ≈ 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 dis- tinction 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 → ∞ 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.

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 differ- ent 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.2 m 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 fluctua- tions 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 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 ∈ [35, 165] cm 3 ∕s , is measured by a flowmeter MAG-VIEW .This provides a direct measurement of the mean velocity U = 4Q∕( D 2 ) that is used to define the Reynolds number Re = Ud∕ of the flow, where = 10 −6 m 2 .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 ∈ [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, flushmounted 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 sensibility 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 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 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 ≠ 0 .As all contributions are addi- tive, 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 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.

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, 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 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 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 Open-FOAM 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 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) 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, , z) , with r ∈ [0, R] , ∈ [0, 2 ] and z ∈ [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 where is the fluid density and 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.
(1) In order to impose a mean flow with periodic boundary conditions, a forcing term is added to the incompressible momentum equation such as where f is the forcing term imposed in ẑ direction.Under this framework, p is a z-periodic pressure field satisfying p(z inlet ) = p(z outlet ) and is related to the total pressure p by the sim- ple relation The equations are then solved in the entire flow domain with periodic boundary conditions in z for u and p at the top and bottom, with no-slip boundary conditions on rigid bounda- ries (reactor walls and beads) for the velocity, and Neumann boundary conditions with zero normal gradient ( ∇p ⋅ n = 0 ) for the fluctuating pressure on solid boundaries.Finally, the full pressure gradient is then computed as ∇p = ∇p + f which averages to ⟨∇p⟩ 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 ∈ [20 − 300] , which is lower than what is achieved in the experiments, and the higher Reynolds values ( Re ∼ 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 = 0.592 .The non-dimensional drag force exerted by the fluid on the sphere is defined as 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, 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.

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 ≤ 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)∕ U 2 , measured experimentally along the cylindrical reactor wall for 9 independ- ent realizations of the bed, for D∕d = 10 and a Reynolds number Re ≈ 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 p∕ U 2 ≈ 10.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 num- ber 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 = cte .It reads: where ⟨ ⟩ r, denotes a spatial average over cylindrical coordinates r and for z = cte .Moreover, ⟨ ⟩ denotes an average over for r, z = cte , and ⟨ ⟩ V an average over the whole volume.Figure 6 shows the pressure fluctuations field p � ∕ U 2 , at two planes at two dif- ferent heights in the bed.As observed in the experiment, pressure fluctuates with typical amplitude a fraction of 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 (5)  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 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 = + ∕Re for each realization from which we get the ensemble average evolution G p = + ∕Re , where , are averages of different val- ues of and , drawn as a solid line in Fig. 7.We then obtain the standard deviation from the ensemble average, dp dz G p = 0.43 with a quite small mean relative error G p ∕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 .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 ⟨G p (r = R, )⟩ z is shown in Fig. 8b for 4 different values of the Reynolds numbers.This was calculated for different angles 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, ⟨G p ⟩ V , which is equal to the forcing term in the DNS, ⟨G p ⟩ V = f .In the present case the dif- ferent 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, 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 ⟨G p (r = R, , z)⟩ with the volumetric average ⟨G p ⟩ V as a function of Re.This is of inter- est 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 ⟨G p (r = R, , z)⟩ ∕⟨G p ⟩ V changes with Re for the differ- ent 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 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, ⟨G p (r = R)⟩ z, Fig. 10 Pressure gradient calculated at the walls averaged over all and z compared to the pressure gradient averaged over the whole volume ⟨G p ⟩ V 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.

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).

Dimensional Analysis
In order to show a robust relation between the pressure gradient in the bed, which is proportional to the total pressure loss Δ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: where L is the relevant length in the tube ( L = H in the simulation while L is the maxi- mal distance between the sensors in the experiment so that the pressure gradient is dp∕dz = Δ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 that can be written as: 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: Note that the condition H∕D ≫ 1 is easily satisfied in the experiment as we have L∕D ∼ 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 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 where (D∕d, ) 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 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 = ∕Re + .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 and 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, and may depend of D/d for wall bounded porous medium.In the particular case where D∕d → ∞ , the Ergun equation is used to describe the pressure drop on non-confined fixed beds of spherical particles (Ergun 1952), which reads: The equation depends on 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.

Results
Figure 11 displays the evolution of the normalized pressure gradient G p = (Δp∕L)∕( 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 0 : G p (R, 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, ⟨G p (r = R, , z)⟩ ,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 = ∕Re + is not a trivial result as there is no reason why the result obtained at moder- ate Reynolds number should be a superposition of the limit cases Re ∼ 0 and Re → ∞ .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 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 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. ( , D∕d) = ( ) × D∕d (D∕d) .Moreo- ver, 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 corre- spond 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 → ∞ limit.As can be  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.

Force Contributions
The DNS data allow for the direct computation of the forces exerted on the spheres by numerical integration of the loads F = ∫ S ⋅ n dS where 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 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 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 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 Spheres = 22 × F Walls , and F Total = 1.05 × F Spheres .Not only that, but both F Walls and F 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 and only depend on and D/d.
Figure 13b shows F Total and F 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).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 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 sig- nificantly 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 ≈ 0.4 .This would be the limit of non-confined beds, where we cite a recent review by Dixon (2023).

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 bulkaveraged 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:// creat iveco mmons.org/ licen ses/ by/4.0/.

Fig. 4
Fig.4The non-dimensional drag force on a three-periodic simple cubic array of = 0.592 as a function of Re.The results are compared with the work of Hill et.alHill et al. (2001b).Inset: the relative error shown as a percentage

Fig. 5
Fig. 5 Normalized pressure p∕( 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

Fig. 6
Fig. 6 Pressure field variations p ′ at two different heights in the bed, for D∕d = 10 , = 0.518 and Re = 200 .Results from Direct Numerical Simulations

Fig. 7
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 G p = + ∕Re obtained by averaging the linear fits G p = + ∕Re of the different realizations

Fig. 8 a
Fig. 8 a Pressure versus z/d at the wall of the reactor ( r = R = D∕2 ) at different angles i , for a fixed bed with D∕d = 10, Re = 200 .b Pressure gradient at the wall of the reactor versus , normalized by the perimeter of the bed, for the same random arrangement as in a. Results from DNS

Fig. 11
Fig. 11 Non-dimensional pressure gradients G p = (Δp∕L)∕( U 2 ∕d) as a function of 1/Re.Left: Experimental results, ≈ 0.4 .Right: Numerical results, with ≈ 0.5 .The error bars for both results correspond to the variability at the walls measured numerically, which increase with the confinement effects

Fig. 12
Fig. 12Pressure drop for cases with D∕d = [5, 5.04] , compared with the Ergun and Dixon models for ( D∕d → ∞ ) and results obtained byErdim et al. (2015).The error bars correspond to the estimated variability of the experimental measurements

Fig. 13 a
Fig.13a 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

Table 1
Parameters of the simulations = V fluid ∕V total : fluid volume fraction.Re = Ud∕ Reynolds number based on the superficial velocity U, the diameter of the particles d and the kinematic viscosity

Table 2
Range of parameters studied experimentally and numerically

Table 4
Forces per surface area in units of N∕m 2 for D∕d = 10 and 5 Forces per surface area [N/m 2 ]