This thesis presents an experimental and numerical study focused on the hydrodynamics in fixed beds of spherical particles that are randomly distributed inside a cylindrical container (e.g. a reactor), so that border effects are not negligible. Fixed beds are a particular example of a porous medium, and we are particularly motivated by the use of fixed beds in the context of AA-CAES technologies as a thermal energy storage unit.

A doctoral thesis is no easy task and I couldn't have done it without the support of a lot of people.

I would like to start thanking Romain, Lionel, Jean-Lou and Mica for having given me this opportunity and for trusting me during this PhD. You are all amazing people, and thank you for all your help in the last couple of years. I knew that I could always have your back and you were all there for me when I needed it. Each of you have helped me in your own way and I deeply care for all of you. Not only are you all excellent researchers, but also excellent people.

Jean-Lou, I think you are a source of inspiration for young scienctists; your knowledge and pedagogy are admirable. Thank you for always taking the time to help me and for all the scientific discussions that we had. Thanks for always respecting and valuing my opinion and for pushing me when needed.

Lionel, the OpenFOAM master, your good mood was and is always welcome, along with your homemade cakes. You were the first one to receive me and I'll never forget how welcoming you were and how you helped me during that first week in IFPEN. I appreciate you always checking up on me to see if everything was okay, especially during the end of the thesis.

Romain, I have learnt a great deal from you, from knowing how things work in the academic system to some life lessons (whether you meant to or not haha). Thank you for feeding my scientific curiosity when I asked random questions and for your dedication. And finally, thank you for always having my back.

Mica, there are not enough words to tell you how thankful I am to you. You were the one who first trusted me from B.A., you opened the doors for this opportunity and once again you helped me with my next steps as a scientist. Thank you for absolutely everything, from being my garant to all the scientific discussions. I knew that if anything happened to me, you would be there for me, and that meant a lot. I would like to thank IFP Energies Nouvelles for having financed this project and for having given me this opportunity. I'd like to particularly thank Olivier Vincke for trusting me, being accessible and for being present when needed.

I would like to thank next to Alexandre Ponomarenko, who helped me with the PTV codes and we worked side by side during one year, and to Nicolas Fintzi for being the best office mate and his help in OpenFOAM. I would also like to thank the Electronic and Mechanic shops from the Laboratoire de Physique, in particular Olivier Razebassia who constructed both experimental setups and was always available when there was any problem, and to Franck Ropars, who helped with the electronics of the pressure sensors.

I would like to thank all the people from the lab at ENS. Every face became familiar and a small part of my day, so to everyone I said "hi" to in the corridors, thank you. It is a very nice V place to work in and the atmosphere is always welcoming and warm. Thank you Jean-Christophe for making the lab such a nice environment to work in.

Living abroad and away from family means that friends become our main support and the people with whom we share everything. I had the fortune of meeting amazing people and becoming friends with them. Marcelo, we've been through hell together and we made it out alive! Your friendship is invaluable and you are like a brother to me. Thank you for absolutely everything: for the laughs, the book talks, doing sports with me and always being there for me. I will always have your llama nearby in order to remind me that I need to have balance in my life. Our talks have literally changed my life! Sebastian, it wasn't so difficult for us to become friends because we both have so similar personalities. Thank you for all the game nights, the shared laughs, the coffee breaks, the beers, the long talks, etc, etc, etc. You have a bright future ahead, and I am sure we will continue to see each other. You are a truly amazing friend.

Facu! First of all, thank you for all your help in my first few weeks in Lyon. After that we quickly became friends and you turned out to be one of the most interesting and loyal people I have ever met. I am sure you will get far in life and that we will cross paths once again. I wish you the best in life, you deserve it. I'd also like to thank Simon for always being there, your company, and introducing the French culture to me. Kris, my partner in crime! Lyon wouldn't have been the same without you. You take me out of my comfort zone and we have the best memories here. Poli, thank you for the hikes in the Alps. One of the best memories that I have from my stay here is when we randomly decided to go to Grenoble after a hike to eat a pizza at a particular restaurant (twice!). Santi, thank you for your kind spirit, you are an amazing friend. Mehdi, you always know the right words to say and you are one of the best people I know. Ele and Bea, you've been here since the beginning, and despite the distance we remain good friends. Cheers to our nights out in the middle of the week! Thomas, thank you for always having a smile on your face and for always teaching me more French. Vamos a la playa! I will never forget our adventures in Greece. Alberto, the other third of the IFPEN office, thank you for your good stories and your good mood. You and Nico made the workdays way more fun! I thank my friends from back home as well, who have accompanied me during this period.

Fixed beds, as porous media in general, present a multi-scale problem: the hydrodynamics at the smallest scales can be very different from those at the largest scales. In order to bridge this scale hierarchy upscale techniques are usually used such as volume averaging, that link the small-scale fluctuations with the large-scale dynamics. These techniques are typically applied at a mesoscale. However, the equations involved at this scale do not form a closed system, so that a closure model is needed.

Moreover, the flows that saturate fixed beds can get incredibly complex, as the spheres distribution and the walls produce tortuous paths and local changes in porosity, and there are solid-fluid interactions that also drive and affect the hydrodynamics. In order to explore how these effects can affect upscaling techniques, we studied the behavior of an inertial confined flow at all three scales: the macro, the micro and the meso. We did this by computing numerical simulations and doing two different experimental campaigns. Both the experimental and numerical methods allowed us to study hydrodynamical effects of the three non-dimensional parameters involved in the system: the Reynolds number Re, the porosity ε and the sphere-to-reactor diameter ratio D/d. The first experimental and numerical campaigns were dedicated to the study of the flow past fixed beds from the global (macro) point of view. We studied the variability of the pressure field and pressure gradient and how they are affected by confinement effects.

The second experimental campaign consisted on studying the hydrodynamics at the pore (micro) scale by using refractive index-matching Particle Tracking Velocimetry. We observed that even though the flow is globally laminar, its behavior at the microscale is comparable to that of a fully turbulent flow. We then complemented the experiments with numerical data, which allowed us to study the local homogeneity and isotropy of the flow.

We finally linked the results obtained at the global and pore scales to the mesoscale using numerical results, which reflected the multi-scale nature of the system. We explored the different terms involved in the non-closed volume-averaged equations in order to help find a closure model that takes into account the wall effects, which include the local solid-fluid interactions and the Reynolds stress tensor.

I II

Résumé

Cette thèse présente une étude expérimentale et numérique portant sur l'hydrodynamique dans des lits fixes de particules sphériques qui sont distribuées aléatoirement à l'intérieur d'un conteneur cylindrique (par exemple, un réacteur), de sorte que les effets de bord ne sont pas négligeables. Les lits fixes sont un exemple particulier de milieu poreux, et nous sommes particulièrement motivés par l'utilisation de lits fixes dans le cadre des technologies AA-CAES en tant qu'unités de stockage d'énergie termique.

Les lits fixes, comme les milieux poreux en général, présentent un problème multi-échelle : l'hydrodynamique à la plus petite échelle peut être très différente de celle à la plus grande échelle. Afin de surmonter cette hiérarchie d'échelle, des méthodes de changement d'échelle sont généralement utilisées, telles que le calcul de la moyenne volumique, qui relie les fluctuations à petite échelle avec la dynamique à grande échelle. Ces techniques sont généralement appliquées à une mésoéchelle. Cependant, les équations impliquées à cette échelle ne forment pas un système fermé, de sorte qu'un modèle de fermeture est nécessaire.

De plus, les écoulements qui saturent les lits fixes peuvent devenir incroyablement complexes, car la distribution des sphères et les parois produisent des chemins tortueux et des changements locaux dans la porosité, et il existe des interactions solide-fluide qui dirigent et affectent également l'hydrodynamique. Afin d'explorer comment ces effets peuvent affecter les techniques de changement d'échelle, nous avons étudié le comportement d'un écoulement inertiel confiné aux trois échelles : la macro, la micro et la méso. Pour ce faire, nous avons effectué des simulations numériques et réalisé deux campagnes expérimentales. Les méthodes expérimentales et numériques nous ont permis d'étudier les effets hydrodynamiques des trois paramètres non dimensionnels impliqués dans le système : le nombre de Reynolds Re, la porosité ε et le rapport diamètre de la sphère/réacteur D/d.

Les premières campagnes expérimentales et numériques ont été consacrées à à l'étude de l'écoulement à travers des lits fixes d'un point de vue global (macro). Nous avons étudié la variabilité du champ de pression et du gradient de pression et comment ils sont affectés par les effets de confinement.

La deuxième campagne expérimentale a consisté à étudier l'hydrodynamique à l'échelle des pores (micro) par la vélocimétrie par suivi de particules en adaptation d'indice optique. Nous avons observé que même si l'écoulement est globalement laminaire, son comportement à la micro-échelle est comparable à celui d'un écoulement turbulent. Nous avons ensuite complété les expériences par des données numériques, qui nous ont permis d'étudier l'homogénéité et l'isotropie locales de l'écoulement.

Nous avons finalement relié les résultats obtenus à l'échelle globale et à l'échelle des pores à la méso-échelle en utilisant des résultats numériques, qui reflètent la nature multi-échelle du système. Nous avons exploré les différents termes impliqués dans les équations non fermées moyennées en volume afin d'aider à trouver un modèle de fermeture qui prend en compte les effets de paroi, qui incluent les interactions locales solide-fluide et le tenseur des contraintes de Reynolds.

III

IV Chapter 1

Motivation

Due to the energy crisis and environment issues such as global warming enhanced by the combustion of fossil fuels, the need to change to renewable energy resources becomes more important and urgent each year. Electrical energy is a commodity used everywhere. Many developed countries are switching to new and renewable sources of energy in order to shift towards a more sustainable energy scenario, as CO2 emissions must decrease by 90% by 2050 to limit global warming below 2℃. Wind, solar and other alternative energy sources are being developed but it is well known that natural sources of energy are intermittent. For example, both solar and wind power depend on meteorological conditions, with solar peaking in the middle of the day, and wind energy peaking in the middle of the night (see figure 1.1). This leads to another major inconvenience: they fluctuate independently from demand. That is, we can have a high electricity generation during an off-peak hour (low demand/high supply) and that energy can remain unused, whereas the opposite can happen as well, where not enough electricity is generated at peak periods. This intermittency can be addressed by developing energy storage units and techniques. This can allow us to save over-generation of electricity to be used at times of high demand [START_REF] Schoenung | Energy storage for a competitive power market[END_REF][START_REF] Ibrahim | Energy storage systems-Characteristics and comparisons[END_REF][START_REF] John | Energy Storage and Its Use With Intermittent Renewable Energy[END_REF]. In particular, we shall turn to Compressed Air Energy System (CAES), which is a kind of storage technology. Energy storage by compressed air has existed since the 1950s, but its energetic efficiency is of about 50%. These are systems where energy is stored in the form of high pressure compressed air and consumed in a different form of energy that has been converted from the compressed air. CAES technology is based on the principle of traditional gas turbine plants, consisting of a turbine, a compressor and a combustor. Gas with high temperature and high pressure, which is formed by the mixing compressed air and fuel in the combustion chamber, drives the turbine which in turn drives a generator to generate electricity. For a CAES plant, there are two different stages of operation: compression and expansion. Since the two stages do not run simultaneously, there is higher system efficiency (48% -54%) than in traditional gas turbine systems (Wang et al., 2017).

The main problem with these systems is that a large amount of the heat generated during the compression process is lost, and the fact that a combustion chamber is needed to heat the air stored in the compression phase, which conducts to a non "green" process. An alternative solution can be achieved using AA-CAES (Advanced Adiabatic Compressed Air Energy Storage), where the heat released from the compression stages can be stored in adiabatic containers and reused during the expansion stages (see figure 1.2 for a schematic representation of the process). In this case, the efficiency can be driven up to 50 -75% (Barbour et al., 2015). These systems can be coupled to work, for example, with wind power generation, as reported in [START_REF] Capecelatro | Controllable and affordable utility-scale electricity from intermittent[END_REF]. It is in this context that IFPEN proposes a system consisting of existing components like compressors and turbines, and new components such as thermal energy storage (TES) systems, where in this case, the TES unit will consist of packed, fixed beds of spherical particles. These fixed beds will be the main focus of this study. Such a system has been explored in Barbour et al. (2015); [START_REF] Sciacovelli | Liquid air energy storage (LAES) with packed bed cold thermal storage -From component to system level performance through dynamic modelling[END_REF], although the use of packed bed as TES unit is still in its infancy. This is why we want to have a thorough understanding of the hydrodynamics in fixed beds so that they can later be applied to modelling problems that can be of use for the development of new AA-CAES technologies.

Chapter 2

Porous media: from the macro-to the microscale Fixed beds of spherical solid particles are a particular case of a porous medium where the fluid saturates the pores between the spheres. This makes them part of a highly multi-disciplinary and wide object of study, as they are present in a wide range of phenomena, from typical industrial applications such as packed bed reactors [START_REF] Wehinger | On the use of refractive-indexmatched hydrogel for fluid velocity measurement within and around geometrically complex solid obstructions[END_REF][START_REF] Petrazzuoli | Numerical prediction of peclet number in small-sized fixed bed reactors of spheres[END_REF], nuclear debris beds (Clavier et al., 2015) and solar thermal power generation (Kumar & Kim, 2017;[START_REF] Gautam | A review on sensible heat based packed bed solar thermal energy storage system for low temperature applications[END_REF], to cryogenics engineering applications [START_REF] Soulaine | Pore-scale modelling of multiphase reactive flow: application to mineral dissolution with production of CO 2[END_REF], medicine [START_REF] Peyrounette | Multiscale modelling of blood flow in cerebral microcirculation: Details at capillary scale control accuracy at the level of the cortex[END_REF] and biology [START_REF] Jensen | Sap flow and sugar transport in plants[END_REF], just to name a few. So apart from the current motivations for AA-CAES systems, their study can contribute to a whole range of phenomena.

Porous media are two-phase systems where there is a solid and a fluid phase. This generates large-scale heterogeneities so that the physics can get incredibly complex. Fluidsolid interactions need to be accounted for and the heterogeneity can give rise to velocity fluctuations over different scales. Because of this, it becomes a fundamental problem that is highly rich: porous media are intrinsically multi-scale, that goes from the smallest scales (the pore-or micro-scale) where local hydrodynamic effects can be observed and fluctuations can arise, to the biggest ones where the medium is driven by the pressure drop and is considered as a continuum. Because of this, modelling porous media requires immense computational power and resources if one wants to resolve all of the scales. Even if a complete solution could be generated, it could provide a far more detailed solution than the one needed at the industrial level.

In order to close this gap between the small and large scales upscaling techniques can be developed, where the system is described at a mesoscale where both micro and macro effects appear in the equations of the fluid and solid phases, and the information from this hierarchy of scales can be propagated [START_REF] Cushman | A primer on upscaling tools for porous media[END_REF].

Object of study

As it has been mentioned, our main object of study will be fixed beds. An example of a simulated bed and a real-life bed are shown in figure 2.2. They are composed of randomly distributed spherical particles that constitute the solid phase and will be saturated by the [START_REF] Wood | Modeling Turbulent Flows in Porous Media[END_REF] to show the multi-scale hierarchy present in porous media. We can start off from the global scale, where averaged quantities dominate the flow behavior (e.g., the superficial velocity U and the pressure drop generated by the pressure), to a meso-scale where we can look at the dynamics in a control volume V, to the smallest scales where only the fluid phase is present. fluid phase. Moreover, they are confined in a cylindrical container (for example, let it be a reactor) so that wall effects might not be negligible, and the height of the bed L is much larger than the sphere diameter so that L/d 1. Taking all of this into account, there are three non-dimensional parameters of interest: the first one is the porosity of the bed ε, defined as the ratio between the fluid volume fraction V f luid and the total volume V total ,

ε = V f luid V total .
This quantity defines the amount of fluid fraction that makes up the bed, and is typically around 0.4 for randomly packed beds of spheres. The second parameter of interest is the ratio between the reactor diameter D and the sphere diameter d, D/d, which accounts for the scale separation of the bed and are going to let us quantify when the border effects are present. When D/d is small, the bed is tightly confined and we would expect to observe wall effects, whereas when D/d 1 the bed is not confined and it could be a priori treated as a tri-periodic bed of particles.

Finally, the last parameter of interest is the Reynolds number Re, which we will define in terms of the fluid viscosity ν, the sphere diameter and the superficial velocity U , which is the velocity we would have if the spheres were not present: Re = dU ν . The different flow regimes present in porous media can be divided into different Reynolds ranges, and citing [START_REF] Wood | Modeling Turbulent Flows in Porous Media[END_REF] we can identify them as1

• Darcy regime (Re < 10): Flow field well approximated by Stokes flow. The viscous forces dominate over the inertia forces and only the pore-level geometry influences the flow.

• Inertial regime (10 < Re < 150): Inertial effects begin to manifest, with pore spaces dominated by inertial cores; the velocity distribution in the pores differs markedly from that of the Darcy regime. Steady vortical structures have been observed.

• Unsteady laminar (transitional) regime (150 < Re < 300): Transitional flows which may be generated by laminar wake instabilities.

• Turbulent regime (300 < Re < 1360): The structure of the flows within the pore space begins to resemble more conventional turbulent flows.

• Asymptotic regime (Re > 1360): The value of Re p is high enough such that the turbulence can be approximated as being locally isotropic over most of the pore space.

We are particularly interested in flows in the inertial and transitional regimes. It is worth noting for example that in the unsteady laminar regime we already have a clue of the scale hierarchy present: laminar wakes can generate instabilities, which means that globally, we can have a laminar regime, whereas the behaviour at the pore scale might be dominated by small-scale fluctuations. Citing [START_REF] Heyman | Scalar Signatures of Chaotic Mixing in Porous Media[END_REF], steady laminar flows through porous media spontaneously generate Lagrangian chaos at the pore-scale.

Upscaling in porous media

There are several methods that can be used for upscaling in porous media, which include time-averaging (Ishii & Hibiki, 2010;[START_REF] Antohe | A general two-equation macroscopic turbulence model for incompressible flow in porous media[END_REF], volume averaging [START_REF] Anderson | Fluid mechanical description of fluidized beds[END_REF][START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF][START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] and ensemble averaging [START_REF] Drew | Mathematical modeling of two-phase flow[END_REF][START_REF] Drew | Averaged equations for two-phase flows[END_REF][START_REF] Zhang | Momentum and energy equations for disperse twophase flows and their closure for dilute suspensions[END_REF][START_REF] Pai | A comprehensive probability density function formalism for multiphase flows[END_REF], but despite the amount of theories and methods that have been developed for multi-phase flows there are still many open questions that have not been yet answered or explored [START_REF] Sreenivasan | Multiphase flows: Rich physics, challenging theory, and big simulations[END_REF].

In particular, one the most widely used methods for upscaling in porous media is the volumeaveraging approach. This formalism is similar to homogeneization techniques [START_REF] Davit | Homogenization via formal multiscale asymptotics and volume averaging: How do the two techniques compare[END_REF] and LES or RANS approaches where the Navier-Stokes equations are expressed in terms of average quantities [START_REF] Pope | Turbulent Flows[END_REF].

This technique is based in capturing small-scale quantities and to relate them to averaged and global variables, such as the mean pressure drop in a bed. In order to do this, it is necessary to have a description of both the fluid and solid phases of the medium and to account for the fluid-solid interactions. The well-known problem with these kinds of methods is that the equations are not closed: there are more unknowns than equations in the system. However, it is possible to manipulate the averaged equations so that the closure problem is reduced to a small number of well-defined and known quantities that can be explicitly related to variables describing the detailed motion of both fluid and particles, so that some microscopic quantities are needed to solve the closure problem [START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF].

In [START_REF] Anderson | Fluid mechanical description of fluidized beds[END_REF]; [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF][START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF] they define an averaging operator that includes a filtering function g and is applied to the whole system. The average of a variable q in a system of volume V is defined as q (x, t) = V q(y, t)g(|x -y|)dV y , which is reminiscent of LES spatial filtering. They average the fluid momentum equation (the Navier-Stokes equation) and the particle momentum equation. The fluid-particle interactions are reflected by the solid-fluid force and the first moments of the force, and they are computed from the particles center of mass. In this case, the averaging is valid if there is a large scale separations, that is, d r L, where d is the sphere diameter, L is the global scale of the medium and r and intermediate scale.

Another widely used volume-averaging method is the one proposed in [START_REF] Quintard | Two-phase flow in heterogeneous porous media: The method of large-scale averaging[END_REF]; [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF][START_REF] Whitaker | The closure problem for two-phase flow in homogeneous porous media[END_REF][START_REF] Whitaker | The Method of Volume Averaging[END_REF]. In this case the averaging is done over a representative elementary volume (REV) and the averaging operator is defined as

q (x, t) = 1 V V q(y, t)dV y ,
where V is the volume of the REV. Here the porous medium is assumed to be made up of periodic REVs and there are also length-scale constraints: the scales associated to the fluid phase have to be much smaller than the one associated to the REV scale L: r L. Moreover, in this formalism local homogeneity is assumed.

Both of these techniques have a strong legacy in the community, as they have been revisited and improved. For example, in [START_REF] Capecelatro | Controllable and affordable utility-scale electricity from intermittent[END_REF] they extended the spatial filtering developed in [START_REF] Anderson | Fluid mechanical description of fluidized beds[END_REF] to formulate a volume-filtered Euler-Lagrange method, and in [START_REF] Jamshidi | On the closure problem of the effective stress in the eulerian-eulerian and mixture modeling approaches for the simulation of liquid-particle suspensions[END_REF] for the modelling of liquid-particle suspensions in an Euler-Euler approach. In Soulaine & Quintard (2014) they use the formalism developed in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] to extend the model in turbulent flows in porous media, and in [START_REF] Wang | Volume-averaged macroscopic equation for fluid flow in moving porous media[END_REF] they apply the same formalism for moving porous media.

These formalisms allow for the equations to be closed, such as in [START_REF] Soulaine | Michel 2014 On the use of a darcy-forchheimer like model for a macro-scale description of turbulence in porous media and its application to structured packings[END_REF] where they close the equations in bi-structure porous media, in [START_REF] Lasseux | Two-phase inertial flow in homogeneous porous media: A theoretical derivation of a macroscopic model[END_REF] they close the system for two saturating fluids, and in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] the Darcy-Forchheimer relation is recovered. However, these kinds of models present a limitation for the present object of study: we do not know, a priori if the separation of scales is valid in our system. Additionally, we do not know if, for example, were we to use the volume averaging technique developed in [START_REF] Whitaker | The Method of Volume Averaging[END_REF], we could represent the fixed bed in terms of periodic REVs, as is typically done in tri-periodic beds of spheres [START_REF] Mehrabadi | Pseudo-turbulent gas-phase velocity fluctuations in homogeneous gas-solid flow: fixed particle assemblies and freely evolving suspensions[END_REF][START_REF] Hardy | Simulation of the flow past random arrays of spherical particles: Microstructure-based tensor quantities as a tool to predict fluid-particle forces[END_REF][START_REF] Tenneti | Particle-Resolved Direct Numerical Simulation for Gas-Solid Flow Model Development[END_REF].

These problems arise because of the geometry of our system. Because we are dealing with confined fixed beds, we expect to have an axial symmetry but radial inhomogeneities are most likely to arise due to finite wall effects. This breaks down the local homogeneity hypothesis proposed in [START_REF] Whitaker | The Method of Volume Averaging[END_REF]. Not only that, but because we have finite D/d we might not be able to satisfy the hypothesis d r L in the radial direction.

In order to explore the finite wall effects (or confinement effects) we are going to study the system at both the global scale and porescale, so as to have a thorough understanding and build an intuition as to how the system will behave when D/d decreases, and if this behavior changes wiht the Reynolds number.

Global border effects

Pressure losses constitute a key parameter governing fluid penetration in packed beds, and understanding various flow regimes in porous media and the associated pressure drop, is critical for its applications. For example, flow details directly influence convection heat transfer, chemical reaction rates, and filtration effectiveness. The pressure drop due to porous-media presence in the path offlow determines the required pumping power in applications utilizing porous media. Fluid flow in porous media is complex mainly due to the presence of the often-random solid structure of the media in the path of the fluid. This structure drastically changes the flow field: it destroys boundary layers and compels the fluid to travel only through winding and tortuous open flow passages. Therefore, the pressure drop in porous media is much higher compared to open-area flows [START_REF] Bağcı | Flow Regimes in Packed Beds of Spheres from Pre-Darcy to Turbulent[END_REF].

In 1952 S. Ergun presented an empirical relation that describes the pressure drop for fixed beds of spherical particles that depends two empirical parameters A and B that take into account the porosity and the particle diameter, and has a leading term that goes with U 2 and a second one that goes with U [START_REF] Ergun | Fluid Flow Through Packed Columns[END_REF]. Even though it has been corrected and readapted to different particle shapes [START_REF] Macdonald | Flow through porous media-the ergun equation revisited[END_REF][START_REF] Allen | Packed bed pressure drop dependence on particle shape, size distribution, packing arrangement and roughness[END_REF], it captures the right correlation between the pressure drop and the Reynolds number and it is usually used as a proxy when calculating the pressure drop. Nevertheless, the Ergun equation applies only for packed beds with negligible wall effects, and deviations occur for finite D/d (Di Felice & Gibilaro, 2004;[START_REF] Cheng | Wall effect on pressure drop in packed beds[END_REF]Eisfeld & Schnitzlein, 2001;[START_REF] Ribeiro | Mean Prosity and Pressure Drop Measurements in Packed Beds of Monosized Spheres: Side Wall Effects[END_REF]Mehta & Hawley, 1969). And the difference is Reynolds-dependent: in the laminar regime, 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 & Schnitzlein, 2001).

No clear correlation has been yet derived for when there are finite wall effects, as we would have to add the wall dependence (for example, in the form of a D/d explicit dependence) into the model [START_REF] Erdim | A revisit of pressure drop-flow rate correlations for packed beds of spheres[END_REF][START_REF] Cerantola | The existence of universal pressure loss and heat transfer correlations for packed beds[END_REF]. In [START_REF] Erdim | A revisit of pressure drop-flow rate correlations for packed beds of spheres[END_REF] they provide a non-exhaustive list of thirty five different models in the literature that propose different pressure drops for different D/d and Re ranges. To derive such a model is not trivial, as there are a priori three non-dimensional parameters, D/d, Re and ε, and there is the additional complication that each arrangement is random in nature, so that there is an intrinsic variability associated with the geometry presented, which can have effects on for example the tortuosity of the flow and might affect the pressure drop measurement, and this does not a priori guarantee the repeatability of the experiments.

In Eisfeld & Schnitzlein (2001), the variability of the coefficients of the Ergun equation was studied by analysing more than 2300 data points with different border effects. The authors found no clear correlation between the empirical parameters A and B 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 & Gibilaro (2004) where the authors define a bulk velocity in terms of D/d. In De Klerk (2003) and Zou & Yu (1995) works, the bulk porosity of the bed has been modelled in terms of D/d as well, by fitting various Mueller (1992). The works referred to in the legend correspond respectively to [START_REF] Roblee | Numerical determination of the volumetric heat transfer coefficient in fixed beds of wood chips[END_REF]; Benenati & Brosilow (1962a); [START_REF] Ridgway | Voidage fluctuations in randomly-packed beds of spheres adjacent to a containing wall[END_REF]; Goodling et al. (1983); Mueller (1992); [START_REF] Giese | Measured and modeled superficial flow profiles in packed beds with liquid flow[END_REF] experimental data.

The wall effects are not only observed in the pressure drop, but also in the global and local porosity of the bed. In fact, the change in the pressure drop can be considered as a direct consequence of the change in porosity. The container used to hold a packing of particles will induce a local area that will make both the micro-and macro-structural properties of particles near the wall different from those far away from the wall (Zou & Yu, 1995). That is, when decreasing the tube-to-particle diameter ratio, the wall exerts a significant influence on the structure in a packed bed: it not only changes the local porosity but also its mean porosity. Figure 2.3 (taken from Guo et al. (2019)) shows different predictions of mean porosity values for packed beds with D/d < 30, where it can be seen that predictions for small tube-to-particle diameter ratios (typically D/d < 10) could vary significantly when using different empirical corrections, and that the porosity of the bed is higher for lower scale separations. In Guo et al. (2019) Not only does the bulk porosity change with D/d, but the wall effects are also evident at a more local scale in the radial porosity profile. In [START_REF] Mueller | A modified packed bed radial porosity correlation[END_REF]Mueller ( , 1992)); De Klerk (2003) they study the variation of ε as a function of the bed radius r, and as shown in figure 2.4 (taken from Mueller (1992)) for different D/d ratios, the porosity near the walls tends to unity, meaning that there is more fluid passing through near the walls. Moreover, an oscillatory behaviour is observed, which is due to the geometry of the bed imposed by the arrangement of the spherical particles. This local variation can have an effect on local quantities: for example, the flow near the walls can be smoother than the one in low-porosity areas, and velocity fluctuations can be enhanced in those areas. This difference also affect the flow distribution: because there is less resistance near the walls (as ε = 1), the flow will have a tendency to flow near them, as opposed to the center of the bed. This leads us to the local effects in fixed beds, which are presented in the next section.

Local hydrodynamics

As it has been mentioned, a flow passing through a porous medium can have a very different behavior at the pore-scale than the one observed at the global scale. In particular, the interaction of the fluid velocity fluctuations with solid particles, as well as particle-particle interactions play an important role in the formation of complex flow structures, such as the ones observed in [START_REF] Patil | Turbulent flow characteristics in a randomly packed porous bed based on particle image velocimetry measurements[END_REF]. By "local" we mean the scales at which the velocity fluctuations are significant and play a part in the dynamics of the flow, because as it has been mentioned, they can be smoothed out using volume-averaging or homogeneization techniques and they might not play a significant role in the global sense (for example, the Darcy-Forchheimer relation works but it does not include the velocity fluctuations).

Nevertheless, the fluctuations and local velocity fields should not need to be neglected at the smaller scales, as they are responsible for various phenomena. For example, mixing in porous media is a clear example of how the velocity fluctuations at the smallest scales can be propagated to the larger ones [START_REF] Villermaux | Mixing by porous media[END_REF][START_REF] Dentz | Mixing, spreading and reaction in heterogeneous media: A brief review[END_REF], and there have been numerous experimental studies on the subject [START_REF] Heyman | Scalar Signatures of Chaotic Mixing in Porous Media[END_REF][START_REF] Jensen | Sap flow and sugar transport in plants[END_REF][START_REF] Le Borgne | The lamellar description of mixing in porous media[END_REF][START_REF] Souzy | Stretching and mixing in sheared particulate suspensions[END_REF], as it is of interest for transport through porous media.

Additionally, experimental studies have been done to study the local velocity field, many at low Reynolds number Re < 1, where we include the experimental works done in [START_REF] Morales | Stochastic dynamics of intermittent pore-scale particle motion in three-dimensional porous media: Experiments and theory[END_REF] and [START_REF] Holzner | Intermittent lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF] where they observe via the probability density functions intermittent dynamics in the local velocity increments and acceleration fields respectively, and in [START_REF] Holzner | Intermittent lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF] they observe that preferential paths develop where the velocity is high next to regions where velocities are much lower; this is probably related to local porosity effects. In [START_REF] Datta | Spatial Fluctuations of Fluid Velocities in Flow through a Three-Dimensional Porous Medium[END_REF] they study the local correlations of the fluctuations of the velocity in the flow for Re < 10 -3 , and they show that the flow is correlated at the pore scale and that the structure of the correlations is determined by the geometry of the pore space, which also affect the tortuosity of the medium, which might affect the preferential paths of the flow.

In [START_REF] Finn | Characteristics of Porescale Vortical Structures in Random and[END_REF] they perform direct numerical simulations for confined packed beds for moderate Reynolds number flows (10 < Re < 200) and they show how different pore scale vertical structures are formed, and vorteces can already be observed at these Reynolds numbers, with different dynamics near the wall of the tube, where the local Reynolds number because of the local increase in porosity. In [START_REF] Patil | Turbulent flow characteristics in a randomly packed porous bed based on particle image velocimetry measurements[END_REF] they perform an experimental Particle Image Velocimetry study for packed beds with a big Reynolds number range, with 228 < Re < 2164, where they observe different structures and fluid patterns at the pore scale, even for Reynolds numbers close to 200, such as "channel-like", "impinging", "recirculation" and "jet-like" regions, each describing a different pattern flow. They also calculate the local velocity correlations, which appear to be once again correlated at the pore-scale, for distances lower than 0.4d.

All these previous studies show how the micro-scale dynamics can differ from the global averaged ones, given that for instance we can have an intermittent behavior even for Re 1, reflected by strong acceleration event and high velocities localized in pore throats, whiel acceleration events are weak in pores where there are almost stagnant velocities; the small scales also affect the global transport of mixtures in porous media, and complex flow structures appear for non-turbulent flow regimes.

Thesis overview

This thesis consists of both experimental and numerical work. Both of these techniques enable us to have different outlooks on the problem, as experiments allow us to for instance validate the simulations, to develop an intuition in the hydrodynamics of the flow and to explore higher Reynolds numbers than the ones achieved in numerically. Simulations on the other hand, will allow us to have access to well-resolved quantities and further information that is not easily measured in the experiments, such as the drag force over the spheres in different packed beds. The parameter space that we will study is the one presented above: the Reynolds number Re, the porosity ε and the particle-to-tube ratio D/d.

The thesis is organized as follows: in chapter 3 we will detail the numerical methodology, including the meshing technique, limitations encountered the validation cases. In chapter 4 the experimental methods are presented. We did two different experimental campaigns, including the building of the set ups from scratch. The characterization of the instruments is presented, as well as the error analysis and the post-processing techniques.

Chapter 5 is be presented in the form of a paper that will be submitted in January 2023. This includes results concerning the global aspect of fixed beds, namely the pressure drop and how it is affected by the confining walls. This concerns the first experimental campaign and a first set of simulations. We will use the experiments to study the pressure drop in a large Reynolds number range (200-1000) and the repeatability of the results. The simulations allow us to have access to the pressure field, and we use it to study the variability of the pressure and how it translates to the global pressure drop.

Chapter 6 concerns local scale results, and includes the second experimental campaign, and a second set of simulations. We will do a thorough experimental study using refractive index-matched particle tracking velocimetry (PTV) measurements to have an access to the local velocity field. This allows us to compute the space correlations of the velocity fluctuations and to perform one-and two-point statistics. The aim of this campaign is to study the local velocity and acceleration fields in a transitional regime (Re ∈ [160,211]).

Chapter 7 is an exploratory chapter where we find some preliminary results in local volume-averaging techniques in confined packed beds, and we use results from the previous chapters to drive the discussion and to analyze these findings. This consists of a third set of simulations where we explore further effects in porosity. The different closure terms are explored in this chapter when border effects are evident. The three chapters containing new results are summarized in figure 2.5, along with the parameters that we varied. Lastly, a final summary, the conclusions and the perspectives are shown in chapter 8

Chapter 3

Numerical Methodology

Governing equations

The system is made up of two phases: a solid and a fluid one. We consider the flow, u, of a newtonian fluid through a packed bed (figure 2.1) so that

∇ • u = 0. (1a)
The momentum equation is described by the incompressible Navier Stokes (NS) equations for a fluid with velocity u, pressure p, dynamic viscosity µ and density ρ, which read:

ρ ∂u ∂t + ρu • ∇u = -∇p + µ∇ 2 u, (1b)
The fluid flows around the solid medium, which is made up of random arrangements of spheres. The spheres are fixed with respect to the fluid and the fluid exerts a hydrodynamic force on them

F h = ˆSp σ • n dS, (3.1)
where S p is the surface of the particle and n its outward unit normal. σ is the stress tensor, which includes the force caused by the pressure, p, and the viscous efforts:

σ ij = -pδ ij + µ ∂u i ∂x j + ∂u j ∂x i . (3.2)
The first moment M of the force distribution can be expressed as [START_REF] Guazzelli | One sphere in Stokes flow[END_REF])

M ij = ˆSp σ ik n k x j dS, (3.3)
which can be decomposed into a symmetric and antisymmetric terms: [START_REF] Greenshields | 2015 OpenFOAM[END_REF](Greenshields, -2015)).

M ij = S ij + A ij . The
symmetric portion S is called the stresslet [START_REF] Batchelor | The stress system in a suspension of force-free particles[END_REF] and is given by

S ij = 1 2 ˆSp (σ ik x j + σ jk x i )n k dS (3.4)
whereas the antisymmetric term A contains the same information as the torque T = ´Sp x × (σ • n)dS:

A ij = 1 2 ˆSp (σ ij x j -σ jk x i )n k dS = - 1 2 ijk T k . (3.5)
The torque plays the same role in generating angular momentum as the force does for linear momentum. The stresslet though, is a less-known quantity that is the result of the resistance of the rigid sphere to a straining motion. Unlike the hydrodynamic force and torque, it is not needed in the equations of motion for a particle [START_REF] Guazzelli | One sphere in Stokes flow[END_REF].

These equations will be used to simulate the flow and to compute the forces and their first moment felt on the particles. This will be useful for studying the closure problem, which includes the coupling between the fluid and the spheres and it is quantified by the forces. The methodology and workflow is explained in the remaining part of the chapter.

OpenFOAM

All the simulations are done using the open source CFD software OpenFOAM, which stands for Open Source Field Operation and Manipulation (Weller et al., 1998;[START_REF] Greenshields | 2015 OpenFOAM[END_REF]Greenshields, -2015;;[START_REF] Jasak | Openfoam: Open source cfd in research and industry[END_REF]. It is a C++ library for the development of customized numerical solvers, including pre-and post-processing utilities, and is mainly used for the solution of continuum mechanics problems including coupling with multiphysics, such as multiphase flows, chemical reactions, turbulence, heat and mass transfers, acoustics, solid mechanics (fluid structure interaction or particulate flows) and electromagnetics [START_REF] Mueller | Consistent closures for euler-lagrange models of bi-disperse gas-particle suspensions derived from particle-resolved direct numerical simulations[END_REF][START_REF] Soulaine | Pore-scale modelling of multiphase reactive flow: application to mineral dissolution with production of CO 2[END_REF][START_REF] Horgue | An open-source toolbox for multiphase flow in porous media[END_REF][START_REF] Petrazzuoli | Numerical prediction of peclet number in small-sized fixed bed reactors of spheres[END_REF].

The software includes certain applications that fall into two categories: solvers, that are designed to solve a specific problem, and utilities, which are designed for data manipulation. As it is an open source library, new solvers and utilities can be easily created by its users with some prerequisite knowledge of the underlying method, physics and programming techniques involved. The overall structure of OpenFOAM is shown in figure 3.1.

We will use OpenFOAM to perform Direct Numerical Simulations (DNS) on periodic fixed beds of spherical particles, which will be either periodic in one direction (mono-periodic) or in all three directions (tri-periodic). This means that the fluid equations are solved without any turbulence model and that all the scales are resolved.

The SIMPLE algorithm

The SIMPLE algorithm will be used for solving the equations governing the fluid motion, within the context of the simpleFoam steady-state solver, which means that ∂u/∂t = 0 in the NS equations (Eq. (1b))). SIMPLE is an acronym for Semi-Implicit Method for Pressure Linked Equations, and as its name implies, it uses the pressure to solve iteratively the steady state momentum equation. The sequence for each iteration n is as follows (Ferziger et al., 2002):

1. Advance to the next iteration t = t n+1 .

2. Initialize u n+1 and p n+1 using the latest available values of u n and p n .

Construct the momentum equations.

4. Under-relax the momentum matrix. This means that u n+1 = u n + α u ũ, with ũ a velocity correction.

5. Solve the momentum equation to obtain a prediction for u n+1 .

6. Construct the pressure equation. As evidenced by the algorithm, the solver solves the equations for each variable u and p sequentially and the solution of the preceding equation is inserted in the subsequent equation.

The non-linearity appearing in the momentum equation is resolved by computing it from the velocity and pressure values of the preceding iteration. When the momentum equation is solved, it delivers a velocity field which is generally not divergence-free (therefore, not strictly incompressible as equation (1a) states). Then, the continuity and momentum equations are used to construct an equation for the pressure, which, if solved correctly, delivers a divergence free velocity field when it is inserted into the momentum equation [START_REF] Hinch | Think Before You Compute: A Prelude to Computational Fluid Dynamics[END_REF]. This is repeated until the results have converged. I invite the reader to read section 7.3.4 of Ferziger et al. (2002) for a more detailed explanation of the algorithm.

As it was mentioned in the description of the algorithm, SIMPLE uses an under-relaxation method. The relaxation factors α u and α p are constants that are multiplying the algebraic equations corrections at each time step. They are used to alter the path towards convergence, preventing the solution from diverging, and they determine the rate of convergence for the simulations. The optimum value of the relaxation factor depends on the mesh and the chosen numerical schemes, and is specific to each particular problem. An underrelaxation factor (less than one) increases the stability of the simulation and an overrelaxation factor (greater than one) increases the rate of convergence at the risk of unstable or diverging solutions.

Residual control will be used to quantify the convergence of the simulations. Residuals show the difference between successive solutions of the equations to be solved. In the case of steady simulations, residuals should decrease as the run evolves, and reach a minimum value (ideally zero) once a steady-state solution is found.

As a last note on the subject, it was seen that the convergence of the solution strongly depends on the number of correctors nCorrectors for the SIMPLE algorithm, which had to be set greater than 2. This is the number of times the pressure equation is solved in a time step, and it is set to 3 in our computations. Besides, in order to reduce the influence of mesh non-orthogonality, the nNonOrthogonalCorrectors are set to 4; the mesh will be explained in section 3.6.

The PIMPLE algorithm

In the case of unsteady flows, where ∂/∂t = 0, the PIMPLE algorithm can be used. The main difference between the SIMPLE and PIMPLE algorithms is that in the latter there are more iterations and corrections when calculating the pressure and velocity fields. Because the problem is time-dependent, the CFL condition has to be established so as to keep a stable solution, which means that, for example, for a 1-dimensional grid with length intervals ∆x, a time step ∆t and velocity |u|, the Courant number C must be

C = |u|∆t ∆x ≤ C max ,
with C max a maximum value, generally chosen to be of the order of or less than 1.

Numerical schemes

The purpose of any discretization is to transform partial differential equations (such as equations ( 1b)) into a corresponding set of algebraic equations. OpenFOAM uses the Finite Volume Method (FVM) to do this, where the discretization of the solution domain produces the positions of points, cell faces and volumes in which the solution is sought and the description of the boundary. The space is divided into control volumes (CV) or cells, and equation discretization gives an appropriate transformation of terms of governing equations into algebraic expressions. Control volumes do not overlap and fill the whole computational domain. Particularly, in OpenFOAM the CV is bounded by a set of flat faces and each face is shared with only one neighbouring CV. A typical control volume is shown in FIGURE, which was taken from H. Jasak's thesis [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF]. The cell faces can be divided into two groups: internal faces, that are the ones between two control volumes, and boundary faces, that are the ones that coincide with the boundaries of the domain. The face area vector S f points outwards from the cell, is normal to the face and had a magnitude equal to the area of the face. On the other hand, boundary face area vectors point outwards from the computational domain [START_REF] Jasak | Error analysis and estimation for the finite volume method with applications to fluid flows[END_REF]. (1996) The fact that the CV can take any polyhedral shape allows us to create unstructured meshes, where the CV do not have a constant volume. We can also apply a local grid refinement, where CV are added in parts of the domain where high resolution is necessary (for example in a gap between two close spheres, or in regions of strong gradients), without affecting the rest of the mesh.

Once the mesh is constructed, we can discretize the integral form of the governing equations over each control volume. The basic quantities (e.g., mass and momentum) are conserved at the discrete level and the equations are solved in a fixed coordinate system (defined by the mesh) that does not change in time. For further reading on the finite volume method, I invite the reader to consult chapter 4 of Ferziger et al. (2002).

Second-order Gaussian integration schemes are used to compute the different terms of the NS equations. As grids are mainly composed of hexahedral cells, the Gauss linear scheme is used for the gradient operator, Gauss linear corrected for the Laplacian schemes, and bounded Gauss linearUpwind is used for the divergence operator (Jasak, 2015). For example, the Gauss gradient scheme calculates the gradient by using integrals over faces:

ˆVP ∇φ dV = ˛∂V P dS φ = f S f φ f , (3.6)
where P is the central point of the CV, φ is the fluid flux passing through a face f of the CV and S is its surface with its corresponding normal. The face value of φ is evaluated at the centers of cell faces, from the cell centers values. In our case, this is done by using a linear interpolation scheme, which uses central differencing for the interpolation.

Boundary conditions

Aside from the fluid equations, we also need to set boundary conditions on the solid objects. The no-slip condition u = 0 for the velocity and Neumann boundary conditions (∂p/∂n = 0) for the pressure are both used on the solid spheres and on the non-periodic reactor walls. Solid walls are so treated by a no-penetration condition, so that the fluid does not pass through them.

In all simulations, there is at least one periodic boundary, for which we will pay particular attention to. This is modeled by simulating an imposed flow rate, which is set by adding a volumetric source term into the stationary NS equations that compensates the pressure gradient and thus the global viscous friction forces on solid bodies. The fluid is so accelerated from rest by a constant body force that conducts to a prescribed average velocity (or equivalently flow rate) in the computational domain. In the particular case of a z-periodic bed with a period of length H the steady momentum equation reads:

ρ(u • ∇)u = -∇p + µ∇ 2 u -fẑ. (3.7)
We can define a modified pressure p as

p = p + f z, (3.8)
with p periodic (p(x, y, z = 0) = p(x, y, z = H)). The pressure gradient is then computed as ∇p = ∇p + f ∇z. When averaged over the whole volume this results as

∇p V = 1 V ˆ(∇p + fẑ)dV.
Because of Stokes's theorem, this becomes

∇p V = 1 V ‹ ∂V pndS + fẑ,
and because of periodicity the first term in the right hand side is zero, so then the averaged pressure gradient becomes

∇p V = -fẑ, (3.9)
which is different from zero. With this formulation, the output pressure that is calculated with the DNS is p. This was done in OpenFOAM by setting a pressure jump between the periodic boundaries ∆p, so that for a spacial period of length H, |∆p|/H = f .

Bed generation

Because we are interested in studying border effects on packed fixed beds of random arrangements of spheres, the main and largest simulations presented in this work are made up of randomly placed spheres in a reactor cylindrical geometry. All the random arrangements of beads in the cylindrical container are built using the Discrete Element Method (DEM) software Grains3D (Wachs et al., 2012). In general, a DEM method is used as a numerical model capable of describing the mechanical behavior of, for example, assemblies of spheres, where each sphere has its own equations of motion. In particular, they follow Newton's second law of motion and they are also described by force displacement laws, such as stiffness and friction. In this case, the contact forces involved are a Hookean elastic force, a viscous dynamic force and a tangential friction force for each sphere. The fixed beds arrangements are computed in three steps: first, the particles are dropped into a non-periodic cylindrical container (which has a bottom solid wall) from an initial point above and outside the container. The spherical particles are driven by a gravitational force until they are randomly settled at the bottom of the cylinder. Afterwards, the two periodic boundaries of the bed are set at the minimum and maximum positions of the settled particles (these will become the inlet and the outlet of the system), and another simulation is ran where the spheres are given an initial random velocity so that they can move inside the now periodic domain and reach a stable configuration. At the same time, their radius is increased at each time step, so that the porosity of the bed decreases at each temporal iteration. During expansion, particles experience multiple collisions before they reach their final diameter that satisfies the maximum solid concentration possible to mesh, and the simulation is stopped just before the particles come in contact with one another. Lastly, once we have the maximum possible radius so that there are no particles touching each other, a third run is done using the particles' last position as an initial condition so that the particles re-accommodate into their final position, i.e. towards a low-energy equilibrium configuration, which will then be the one used for meshing. These steps are illustrated in figure 3.4. The beds simulated in the present work contain between 72 and 2600 spheres each, and the porosity varies between [10 -50]%.

The resulting beds are characterized by calculating fluid volume fraction as a function of r, denoted as ε(r). This is shown in figure 3.5 for three particular cases with D/d = [5. 13, 8.03, 10.15]. The results are compared with those obtained by Mueller (1992); Benenati

Step 1: Random insertion of spheres

Step 2: Spheres are settled

Step 3: Insertion of periodic boundaries

Step 4: Spheres move randomly inside the periodic domain and increase their radius at the same time

Step 5: Spheres with the maximum possible radius are reaccommodated The results are compared with those obtained in Mueller (1992), Goodling et al. (1983) and Benenati & Brosilow (1962b).

& Brosilow (1962b); Goodling et al. (1983), where these authors measured experimentally the porosity profiles in beds whose sphere-to-reactor diameter aspect ratio are similar to ours. It is observed that the results match, even though we are not under the exact same conditions. Once the spheres arrangement is computed and characterized, we have to mesh it, as explained in the following subsection.

Meshing strategy

Overview

As already stated, the mesh plays a very important role in the convergence and accuracy of a simulation. A simulation can indeed diverge or give non-physical results if the mesh is not of a sufficient quality (e.g., highly skewed faces, high aspect ratio) or too coarse, because of the poor quality of the solution. Moreover, it is particularly challenging in the case of porous media as it is necessary to mesh the solid-fluid interface and resolve the small gaps between solids. We will construct unstructured meshes with OpenFOAM utilities as detailed in the next section. This means that the control volume have different volumes and that the mesh is not orthogonal nor Cartesian (hence the need for the nonOrthogonalCorrectors used in section 3.2). All the simulations will have at least one periodic direction. Non-matching grid periodicity is used and is handled by the cyclicAMI condition, as detailed in subsection 3.6.

The general workflow of the meshing process is summarized in figure 3.6. It consists of five steps. First, the blockMesh utility is used to generate a background grid. We then use snappyHexMesh utility, which generates automatically three-dimensional grids containing mostly hexahedra (hex) and split-hexahedra (split-hex) from triangulated surface geometries that in our case will be the spheres. The mesh approximately conforms to the surface by iteratively refining an initial background grid (in our case, made with blockMesh) and morphing the resulting split-hex mesh to the surface. The specification of mesh refinement levels is very flexible and the surface handling is robust with a pre-specified final mesh quality. The background grid defines the refinement level 0. The level 1 corresponds to a division of level 0 cells by a factor 2 in all directions, and so on for higher levels. In order to resolve better the small gaps, we can choose to add refinement levels specifically in these regions. This will be further detailed in section 3.6

The createPatch utility allows us to generate the periodic conditions. It creates cyclic patches out of selected boundary faces using the cyclic AMI (Arbitrary Mesh Interface) method. This method conducts the solvers to interpolate data on coupled periodic boundaries, even if the points and faces are non-matching. We can specify in its input dictionary the cyclic conditions along the translational vector.

The renumberMesh is then used. It renumbers the cell list in order to reduce the linear solver matrix bandwidth, which will improve simulation stability and speed. And lastly, we use checkMesh to check the quality of the obtained mesh. This raises warnings or critical errors when there are problems with the periodic boundaries, if there are highly skewed or non-convex cells and whatever issues might arise from a bad quality mesh. If this check fails, the simulations are likely to crash.

Periodic meshes

As it has been mentioned, we will use the cyclicAMI boundary condition, which couples a pair of periodic boundaries whose inner points and faces do not exactly match, by interpolating the results at each periodic cell within a set tolerance, which we set to be 0.0001. However, points and faces outer bounds should geometrically match. This makes it highly mesh-dependent, because the periodic boundaries have to be as similar as possible so as to generate a reliable periodic solution. This becomes more difficult when there are solid obstacles (the spheres) near or across the periodic borders. In this case, the spheres need to be duplicated in each periodic direction, so that the mesh can be generated identically (or as close to identically) on each side. That is to say that each periodic sphere has to have its (periodic) counterpart. An example of this is shown in figure 3.7, where the unit cell is a Face-Centered Cubic cell with periodic spheres on the walls, all with their periodic homologous counterparts.

Moreover, aside from repeating the spheres and their periodic homologous, explicit edge capturing is used in the snappyHexMeshDict dictionary by providing eMesh files. These files represent discretized circles that result from the intersection between a sphere and the border of the domain, which is then reproduced in the periodic direction. These discretized circles provide snappyHexMesh the points/segments that the generated mesh is required to match. Besides, refinement levels are also explicitly specified on those edges. This is done so that the periodic boundaries can match as close as possible, and so that the interface between the spheres and the walls can be well resolved. A visualization of a sample eMesh file is shown on figure 3.8a). The eMesh is shown in solid red, and it can be noted that there is one more level of refinement along it.

Another thing that has to be taken into account for a good AMI interpolation is to verify that there are no tangent spheres near the periodic boundaries so as to avoid any strongly skewed cells that could affect the convergence of the results. We do this by manually setting the periodic boundaries at a place where the spheres are neither too close nor too far away a) b) Snapshot of a sphere in the border of a periodic boundary and its eMesh file, which represents the circle that results from the intersection of the sphere and the wall of the domain. b) Superposition of the grid in a border of the domain and its periodic homologous counterpart from the borders. Finally, in order to check the quality of the interpolation, the AMIWeight quantity is shown in the checkMesh when the mesh is evaluated. This interpolating weight quantifies the similarity between a pair of periodic boundaries. It is equal to 1 for a perfect match and 0 in the opposite case. In all the simulations done in this work, its value is 0.99 in average. Figure 3.8b) shows two superposed periodic grids: one in red and the other in blue. The lines are not easily distinguishable because the grids match almost perfectly.

Reactors

Cylindrical reactors are meshed with the geometry obtained by Grains3D (section 3.5). Reactors are periodic in the main flow direction because of two main reasons: the first one is to avoid any inlet/outlet boundary condition effects, and the second one is to avoid large meshing and computational times. Grids have typically tens of millions of cells and this number increases as 2 n (n is the level of refinement in the snappyHexMeshDict ) within the small gaps.

The background mesh is computed with blockMesh using a butterfly O-H topology grid, as shown in figure 3.9 along with the coordinate system. This is based on creating the H square-like shape in the middle of the circular section surrounded by O topological blocks that connect the diamond with the borders of the circular cross section. The main issue with the butterfly structured mesh topology is that the cell sizes cannot be strictly uniform everywhere, in particular along the radial and azimuthal directions. This can easily be observed in the left snapshot in figure 3.9. Nevertheless, it was ensured that the maximal cell sizes of the background grid were below a specified resolution, so that everything is well resolved in all parts of the domain. Surface refinements are specified in the snappyHexMeshDict, and are set to be between levels 1 and 2 on the sphere surfaces, and to 2 on the edges at the intersections of spheres and periodic boundaries, specified as eMesh files entries. The number of cells between each refinement level (entry nCellsBetweenLevels in the snappyHexMeshDict dictionary) is set to 4. The additional refinement level in very small gaps, namely the gapLevelIncrement keyword in the dictionary, is set to 2. These different refinements can be appreciated in figure 3.10. 3.10a) and b) show two slices of a standard fixed bed, where we can observe the level 1 refinement at the spheres' surfaces and the level 2 in the small gaps between them. 3.10c) shows a snapshot of a periodic boundary, with the reactor cut in half. The edge refinement is at level 2 as well, and with this taken into account the AMI weights are 0.99 in average for all cases. Finally, figure 3.10d) shows a generic transversal cut of the grid, where the different levels of refinement can be appreciated as well as the packing concentration. The spheres below the cut are shown for illustrative reasons.

In average, 85% of all the cells are hexahedra, which makes the grid suitable to use a gaussLinear gradient scheme. There are 0.07% of under determined cells in average, and 8% of concave cells. The percentage of highly skewed faces is negligible, being only 5.3e -5% in average. The meshes have between 45 and 90 million cells, and they have a background grid resolution of at least 30 cells per sphere diameter, varying on the case. The relevant remaining parameters (number of cells, resolution) will be detailed once the results are presented.

Validation

Before going on to more complex simulations, it is necessary to validate the whole workflow, i.e. meshing and resolution of the Navier-Stokes equations. The different steps in which we did this are detailed in the following subsections.

Force over a single sphere

As a first test case, we computed the force f D over an ordered array of spheres and compared it to the results obtained by Hill et al. (2001). The system consists of a Simple Cubic (SC) ordered array of spheres, with a porosity ε = 0.592. This is done by creating a unit cell with a sphere in its center. The boundaries are all periodic, so that we used the cyclicAMI condition for all variables on domain's limits. A no-slip and Neumann boundary conditions were imposed on the sphere wall for the velocity and pressure, respectively, and the Reynolds number explored varies between 0 and 501 .

Problem description

The pressure and fluid velocity are governed by the incompressible Navier-Stokes equations, introduced in section 3.1, with the no-slip boundary condition at the sphere surfaces and periodic conditions at the bounds of the computational domain.

We recall that the Reynolds number is defined in terms of the sphere diameter d, the superficial velocity U and the kinematic fluid viscosity ν as

Re = U d ν .
The stram-wise component of the non-dimensional drag force F is defined as

F = f D 3πµdU . (3.10)
The denominator on the right-hand side of equation 3.10 is the Stokes drag force on a single sphere in an unbounded fluid, and hence the positive deviations of F from unity indicate the contribution of hydrodynamic interactions (finite solid volume fraction) and fluid inertia (finite Reynolds number) to the magnitude of the average drag force. At moderate Re, the spatial structure of the velocity field depends strongly on the solid volume fraction, the Reynolds number and, for ordered arrays, the direction of the average pressure gradient relative to the axes of the arrays (Hill et al., 2001). The results obtained in Hill et al. (2001) are noted as F H and ours as F .

Mesh convergence 20 40 60 80 cpd Because the mesh plays such an important (if not the most) role in the simulations, it is imperative that it is resolved well enough so that it is able to capture the physics of the problem that is being simulated. In this case we will quantify this by studying how much the results vary for different grid resolutions, which are represented in cells per particle diameter d (cpd from here on) and the results do not have to depend on the mesh resolution. The resolution of the simulations were done as detailed in section 3.2. For a uni-dimensional mesh with length L and N cells, the cpd is calculated as cpd = d×N L . We compare the results for 17 different background grid resolutions with the previous reference, for the case of the sphere in the middle of the domain and Re ≈ 50 and then compute the drag force felt by the sphere. The results of the mesh convergence are shown in figure 3.11, represented by the error between the force calculated by OpenFOAM and the results obtained by Hill, as a function of the cpd. The results and the error stabilize at ∼ 2% when the mesh resolution is at least 25 cells per sphere diameter. So it can be concluded that we have reached a mesh-independent physical result, that the simulations are well resolved, and the results are accurate and reliable.

10 -2 10 -1 10 0 |F -F H | F H [%]
It is worth noting though, that the mesh resolution in the reactor simulations will be higher as it will vary between the resolution set by the background mesh and the one when all the small gaps are taken into account at different parts of the bed. That is, the cpd will be between

d × N L < cpd < d × N × 2 n L ,
where L is a characteristic length of the reactor (it can either be its diameter D or height H), N the total number of cells in a given direction (r, θ or ẑ), and n the refinement level.

Tri-periodic sphere U We will now concentrate on the case of a single mesh resolution, with 36 cells per sphere diameter, and the surface refinement around the sphere is set to 2. The comparison between our results and those obtained by Hill et al. (2001) are shown in figure 3.12, where it is evident that there are no significant differences between both results. This is further illustrated by the relative error, defined in terms of the force calculated using OpenFOAM F and the one calculated by Hill F H as, Error = |F -F Hill |/F . This is shown in the inset, where we obtain a maximum relative error of 1.8%.

n = r n ∥ U n ⊥ U 10 20 30 40 50 Re 1.0 1.5 2.0 |F °FH | F H [%] n k U n ? U n = r a) b) (2) (3) (1) (1) (2) (3)
Moreover, in order to assess the behaviour of the cyclicAMI condition, we computed three different SC arrangements: the one already presented, with the sphere at the middle of the domain (a), another one with a periodic sphere cut into two in the stream-wise direction direction (b), and a third one with its periodicity perpendicular to the flow (c). These three cases are illustrated in figure 3.13a). The relative errors are shown in figure 3.13b). Errors are never higher that 2%. As expected, the error increases a little bit when the spheres are in the border, because of the AMI interpolation. It is interesting to note that the case with the highest error (albeit small) is when the spheres are periodic in the direction that is perpendicular to the flow.

Poiseuille flow

Another test was done to validate the methodology, this time by calculating the flow with dynamic viscosity µ through a mono-periodic pipe of radius R and height H with no beads. We expect to recover the results from the Hagen-Poiseuille equation for the axial component of the velocity [START_REF] Sutera | The history of poiseuille's law[END_REF], which reads .11) where ∆p is the pressure difference between the inlet and the outlet. A no-slip condition (u = 0) is imposed at the walls, along with Neumann boundary conditions for the pressure. The velocity at the inlet and outlet of the tube are periodic in the axial (z) direction, and a constant body force is applied in the z direction. The Reynolds number is varied between 10 and 150, and the length of the tube is 20 times its diameter, so that H = 40R. This is shown in the left snapshot of figure 3.14.

u z (r) = ∆p (R 2 -r 2 ) 8µH , ( 3
The velocity profiles obtained are shown in solid lines in figure 3.15 for Re = [10, 50, 100, 150], and compared with the analytical result, equation (3.11) shown in the red star markers. It is observed that the results match almost perfectly. The average relative error is 2%. The radial coordinate is normalized by the pipe's radius R and the velocity profile u z is normalized by u max = u z (r = 0) so that all the curves collapse. In the right panel of figure 3.14 a solution for the velocity field can be observed for Re = 10.

Spurious periodic effects

In order to verify that there are no spurious effects generated by the length of the grid (for example a coupling from the walls), we simulated two reactors with 2H length for D/d = [5, 8], calculated the pressure drop, and compared the results against a case where the period is H. Given that the physical problem is the same, there should not be a difference in the pressure drop between the two cases, as it does not depend on the periodicity of the system.

The averaged pressure drop is calculated as explained in section 3.4, following equation (3.9). There are no significant differences between the pressure drop results: for the D/d = 5 case we have that (∇p 2H -∇p H ) ∇p H = 0.23%, whereas for D/d = 8, (∇p 2H -∇p H ) ∇p H = 0.20%, with both cases computed for Re = 60. 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 steady state solver, the solutions once the problem starts becoming transient do not differ that much from their average value [START_REF] Finn | Characteristics of Porescale Vortical Structures in Random and[END_REF][START_REF] Hill | The transition from steady to weakly turbulent flow in a close-packed ordered array of spheres[END_REF]. We verified whether the non-stationary modes are weak by comparing results obtained with SIMPLE to those given by the non-stationary solver PIMPLE, for a particular case with Re = 200 and D/d = 8. The PIMPLE simulations are computed with a variable time step in such a way that the maximum Courant number is 0.5. The final residuals are 10 -8 for both the velocity and pressure, whereas they are of O(10 -5 ) for the SIMPLE solution. An Euler-implicit time scheme is used to resolve the time derivative, ∂u/∂t, which is a first order transient scheme.

SIMPLE vs. PIMPLE for high-Re flows

1 2 3 4 t/T H 1.3E -04 1.5E -04 |f P IM P LE -f SIM P LE | f SIM P LE [%]
We compare the pressure gradient obtained at each time step with the one obtained from the stationary simulation. The time in the simulations is normalized by an integral time which is defined as the time that the fluid needs in average to travel through the whole domain of length H: T H = εH/U . As can be observed from figure 3.16 all the results match down to an error of O(10 -4 ), thus showing that a global quantity such as the pressure drop is accurately resolved by the SIMPLE algorithm.

So as to compare local quantities, we computed the fluid velocity's probability density function (pdf) at different time iterations t i up until 4 integral time units, shown in colors in figure 3.17, and compared it to the one obtained using SIMPLE, which is shown as a black solid line. Once again, there are no significant differences between the results. Moreover, we can also compare the total forces on the spheres calculated by both algorithms. Figure 3.18 (Top) shows the relative error defined as E = |F pimple -F simple |/F simple for 13 random spheres of the bed as a function of time, and F simple and F pimple are the forces computed by SIMPLE and PIMPLE respectively. It can be seen that the error remains relatively low, and that there are peaks where it can be up to 7% for some spheres. We tried to compute the Fourier spectra of the forces but it could not be exploited, as it seems that the signal is chaotic, or cannot at least be expressed as a convolution of different harmonics. Figure 3.18 (Bottom) shows the time-averaged error for the 13 spheres. Its value is 2.5% at the most, indicating that the SIMPLE algorithm provides an accurate enough solution.

Experimental methodology

In this chapter, the experimental techniques and methodology will be presented. It is divided in two parts: section 4.1 will be dedicated to the global measurements and section 4.2 to the local measurements. In both, the experimental setups are presented, along with the acquisition techniques, characterization of the instruments, signal processing, postprocessing analysis, and the shortcomings encountered.

Global measurements: pressure drop in fixed beds

The main objective of this experimental campaign is to measure the pressure gradient in confined fixed beds of spherical particles. To this aim, we will measure the pressure drop in the experimental device described in the following section, for different particle sizes and different Reynolds numbers. Therefore, the parameters that will be varied are the particle Reynolds number Re p , which we recall to be

Re p = U d ν ,
and the particle-to-reactor ratio D/d, with d the particle diameter, D the reactor diameter, U the velocity of the flow as if the spheres were not present (that is, the superficial velocity), and ν its viscosity. The superficial velocity is related to the flow rate Q as U = 4Q/(πD 2 ) , and it will be the parameter that we will vary in the experiment. The variation of the Reynolds number will allow us to study how the pressure drop depends on Re in transitional regimes, whereas the variation of the particle size will allow us to see how the effects of the borders and particle distribution come into play.

In the following, the experimental setup will be introduced and in the following subsections all the instrumentation that was developed for the characterization of the instruments used will be presented. It was necessary to design an experimental apparatus from scratch where we could measure the pressure drop, which is a global (or macroscopic) quantity. To do this, we designed a setup as the one schematically shown in figure 4.1. This was constructed by Olivier Razebassia from the Mechanic Shop at the Laboratoire de Physique of ENS Lyon. It consists in a cylindrical test section made of a clear plexiglass column of diameter D = 0.04m that has several piezoresistive pressure sensors distributed vertically along its walls, each at a distance ∆h from each other, so that we will be able to measure a pressure gradient ∇p = ∆p/∆h, with ∆p the pressure difference between two sensors.

The main priority that we had when designing the setup was that it should allow us to be flexible while carrying out the experiments. That is, we did not know a priori if there were going to be any nozzle effects such as an undesired acceleration of the fluid [START_REF] Singh | Effect of nozzle geometry on critical-subcritical flow transitions[END_REF], how high from the inlet the bed should be located, if we needed to add another porous medium to laminarize the flow between the inlet and the entrance to the porous medium, and so on. This was achieved by constructing several sections of the cylinder that could be screwed and unscrewed, with different heights and several grids that would be useful for different purposes. Several sections of the experiment are shown in figure 4.2. This also makes it easier to change the bed at each realization of the experiments, as it is not necessary to disassemble the whole setup. The test section containing the bed is part of a closed upstream-flow water loop that also includes a pump and a flow meter, and has a height of L exp = 0.4m and a diameter D = 0.04m. Filtered water flows upstream through the circuit (the inlet is at the bottom of the reactor). A thermal bath is used to control the temperature, thus allowing us to have a fixed and controlled viscosity. Given that each time that the particle distribution is changed, some air is inserted into the circuit, particular attention must be paid to the removal of bubbles. This is a typical issue encountered in these kinds of setups when we have to constantly assemble and disassemble each time that we do a new experiment run. The bubbles are not only bothersome, but they can also change the physics of the problem, given that they add a third phase into a two-phase (solid and liquid) system, and could consequently change the physical results.

To this end, as a first step there is a reservoir on top of the reactor so that the bubbles can escape the closed water loop when reaching the free surface of the tank. Additionally, the filtered water is degassed by leaving it in circulation at 55 °C for nearly 12 hours and then cooled down to 20 °C, the temperature at which the experiments are carried out. This facilitates the suppression of bubbles and it was corroborated visually before starting each acquisition. In total, each experiment took about 24 hours each time the bed is modified.

In order to make sure that the flow was laminar before entering the bed and that there were no nozzle effects, there is a honeycomb before the bed, between two empty tube sections. The honeycomb is 0.05m high, whereas the two empty sections are 0.2m high, placed before and after the honeycomb. This was done after observing a jump in the pressure curves at a certain Reynolds numbers, and we believe that it is because there was an instability coming from the nozzle. Because the nozzle is 4 times smaller than the cylindrical test section, from the conservation of mass we have that U nozzle = 16U cylinder , and so we can have two very different regimes near the nozzle and near the bed. Once the honeycomb was added, this jump in pressure disappeared. As already mentioned, this setup is intended to be used to measure global quantities, such as porosity and pressure profiles. The pressure sensors allow us to have ∆p/∆h, and the flow meter allows us to have the velocity U = Q/A, with Q the volumetric flow rate and A = π/4D 2 the area of the cross section of the (empty) bed. Figure 4.3 shows the test section of the setup with a partially filled bed, where the pressure sensors attached to the cables can be seen. The specifics of the sensors are detailed in section 4.1.2.

Both the pressure sensors and flow meter have an analogue output, and all of the signals were measured using a National Instruments PXI-1002 card which provides 8 analogic input channels, has a 23-bit resolution and measures at a frequency of up to 100kHz. The acquisition is done using the LABView software (which is developed by National Instruments), and the post processing is done both in MATLAB and Python.

The flow meter is connected to two other components. The first one is a signal attenuator, and the other is a low-pass filter. The attenuator is used because the voltage range that the analogue output from the flow meter was in a different voltage range than the one accepted in the acquisition card, and the low-pass filter was used to reduce ambient electromagnetic noise from the signal. This will be detailed in section 4.1.2.

Stainless steel spherical beads which have a density of 8g/cm 3 from Marteau et Lemarie (France) are used. They were chosen so as not to have any fluidization. This material was chosen because of its high density, as previous tests were carried out with glass beads and they all got fluidized at the Re numbers that we are interested in. 1 In addition, a grid on The flow rate is one of the main control parameters which has been used to tune the Reynolds number. It is therefore crucial to have an accurate measurement. A magnetic flow meter was used to measure the volumetric flow rate, model MAG-VIEW MVM-020-QA from Bronkhorst, which is an electromagnetic inductive flow meter for electrically conducting liquids with a minimum conductivity of 20µS/cm, such as water. Its operating range is between 1 and 20 l/min, its response time is smaller than 100 ms, and it has an operating temperature between -10 to 60℃. It also has an analog output which delivers a square wave with its main frequency proportional to Q, with a factor F * of 1000 pulses per liters/s.

A raw the signal obtained with the N I -Dacq is shown in figure 4. 6 (left). As the signal is noisy and as was mentioned in the previous section, a commercial low-pass filter Standford Research (SR640) was applied so that only the first harmonic of the signal is obtained, along with a signal attenuator. This was done so that the signal was between 0 and 3 Volts. The which the bed starts getting fluidized. The FFT of the spectre measured by the pressure sensors (which captured the oscillating behavior of the bed) was also computed to get the principal frequency of the kinetic wave generated by the movement of the beads. The flow rate is then proportional to the signal frequency f max and is given in cm 3 /s by

Q = 1000 × f max F * .
f max is the amount of pulses per second of the flow meter. This is computed by making a power spectral density (PSD) estimate of the flow meter signal. The PSD describes how the power of a signal in time is distributed over frequency. We use a Welch algorithm with several estimates of PSD in different windows of the signal which are then averaged to reduce the variability of the spectral estimate. In figure 4.7, the spectrum of the flow meter signal is shown, which has only one mode due to the usage of the low-pass filter and its maximum gives f max . Once the flow rate is computed, it is straightforward to calculate the velocity, given that

U = Q/A
(for incompressible flows). A is the cross-sectional area of the (empty) reactor. The resulting velocity is the so-called superficial velocity (the velocity as if the spheres were not present) and is the one that is used in the definition of the Reynolds number. For each velocity measurement, the flow rate was determined from a 60s reading from the flow meter.

The flow rate is fine-tuned by a solenoid control valve Bürkert-Type 2875 which has an operating pressure range between 0 and 25 bar and a 24V DC operating voltage. It is a normally-closed valve and its opening is proportional to the voltage. It is controlled by a Bürkert Type 8605 controller that converts an external standard signal into a pulse-width modulated signal, which enables delicate adjustments of the opening of the valve.

(Global) Porosity measurements

The porosity of the bed, which is defined as the ratio between the volume occupied by the solid (the beads) and the fluid in the test section,

ε = V solid V f luid ,
is measured by knowing the weight of the bed in air. In this case,

V air = V empty reactor = πL exp (D/2) 2 ,
and the weight of the solid bed in air is straightforward because we know the beads' density, so

V solid = V beads = m beads /ρ beads , hence ε = m beads πρ beads L exp (D/2) 2
Pressure sensors and signals We used high-precision piezo-resistive sensors model MS5401-BM from TE Connectivity to measure the pressure. They have a sensitivity of 150mV/bar with a linear deviation of 0.05% with a 1 bar range.

The sensor includes a silicon membrane where borosilicate glass wafers are vacuum glued to the back for a reference pressure. All of this is installed in the center of a metal capsule filled with a silicone gel and directly mounted onto a printed circuit where the acquisition cables are plugged into. Beause of this coating the sensors break if they come into contact with the beads, so we had to isolate them from the bed with the use of washers that were added onto the setup where the captors are installed. The electronic components were mounted by Frank Ropars from the Electronic Shop at LPENSL in the context of J.Vessaire's PhD thesis (Vessaire, 2019).

The advantage of piezo-resistive sensors compared to piezo-capacitive is that they allow the user to measure both mean pressure and fluctuations. To this end, a voltage amplifier was used with a gain of G = 10 so it holds now that 1 bar corresponds to 1.5V. The signal is digitized and acquired by the N I -Dacq.

The calibration of the sensors (allowing to convert the voltage signals into bars) was done several times for all of the sensors and then those values were averaged. An example of the calibration of one sensor is shown in figure 4.8. This was done by measuring the hydrostatic pressure p in a column of water with variable height z, so in general we have that p = ρgz+p 0 . In this case, we should have a measured pressure p that goes as p = γρgz +δ, with γ a scaling factor and δ an offset. We obtained a calibration factor of 1.49e-5 V/Pa = 1.49V/bar, which is consistent with the sensor specifications and the amplifier.

In figure 4.9, a typical pressure signal p is shown (in Volts). The aspect ratio between its standard deviation σ and its mean µ is σ/µ = 0.03, so there is a variation of about 3% around the mean value. As to see how well converged the first and second moments of the signal are, in figure 4 where P is the probability density function. This is another way of measuring the moments of the distribution and as the both reach a plateau, it is safe to assume that both the mean value and standard deviation of the pressure are statistically well converged. The pressure gradient ∇p is then calculated as the slope of the p(z) curve. 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. An example of such curve is shown in figure 4.11. The pressure acquisition was done for 120s and time averaged for each point.

Summary

Having presented the setup and the characterization of each of the instruments, the summary of a typical run is as follows:

1. Insertion of the spherical beads into the test section. We ended up doing experimental campaigns for three monodispersed arrangements with three aspect ratios: [5,8,10]. The beads were poured into the section and we then shook the container.

D/d =
2. The porosity of the bed is measured by weighting it on a scale.

3. We let the filtered water circulate at 55 °C and then cooled it down to 20 °C so as to degas the system.

4. We vary Re with the use of the solenoid valve and we measure the velocity with the flowmeter, which is connected to the NI acquisition card.

5. The pressure gradient is computed by measuring the pressure with the high-precision piezo-electric sensors at different heights of the bed, and we only used the sensors that were far away from the inlet and outlet. Their signal is recored by the NI-Dacq as well.

6. The post-processing of the data is done in MATLAB and Python.

Local measurements: pore-scale hydrodynamics

The aim of this experimental campaign is to study the local velocity field at the pore-scale.

In order to do this, we used Particle Tracking Velocimetry (PTV) using index-matched hydrogel beads to make the porous bed and we studied the velocity fluctuations, which we recall to be defined as

u = u -u , (4.1)
where u is the velocity of the fluid that saturates the bed.

By looking at the fluctuations, we will have a thorough description of the hydrodynamics at the local scale, and we will be able to calculate one-and two-point statistics, such as the velocity's probability density functions (pdfs), correlations and structure functions.

In the following sections the experimental device will be described, followed by the characterization of hydrogel beads that were used as the solid phase, a first characterization of the bed, the details of the PTV algorithm and a trivial case that was used to validate the data acquisition and post processing.

Experimental device

As was the case for the other experimental device, it was also necessary to design and build a setup from scratch. This was done with the help of Olivier Razebassia from ENS Lyon. The base of the set up is similar to the one described before: it consists of a cylindrical test section of diameter D = 0.09m, with an closed upstream-flow water loop, a pump, a solenoid valve, a flow meter and a reservoir at the top to help to get rid of the bubbles, with the same instruments that were used for the other setup, and the characterization that was done previously holds. A schematic representation of this section of the device is shown in figure 4.12.

Because we are going to carry out optical measurements, we have to pay special attention to the possible reflections and optical index variations. As we have a circular cross section, any incoming light from a laser sheet will be reflected and refracted by the uneven surface. In order to solve this, we added a rectangular vessel surrounding the cylinder where the bed is put so that we can fill with water and drain it independently from the closed loop of the primary section of the device. This solves the issue of having optical distortions due to the curvature of the test section, impairing both illumination and visualization.

The illumination system consists of a 5W Quantum Opus 532 laser with a 532nm wavelength (green). The beam enters through an light-guiding arm (Dantec dynamics) that has an internal alignment mirror system that maintains the coherence of the laser beam, and can be moved in several degrees of freedom without affecting the beam. This is highly practical for PTV measurements because it allows us to choose the place where we want to put the laser sheet.

At the end of the arm, there is a cylindrical lens so that the beam is expanded and a vertical laser sheet is generated, illuminating a plane that is parallel to the flow rate. This sheet, of approximately (5 × 0.4)cm (length×width), illuminates into a test section of the fixed bed. The lens is attached to a stable, linear translating stage Zaber that has one degree of freedom that allow us to move the sheet over 75cm perpendicular to the flow rate plane (depth) in order to scan different zones of the bed.

The water is seeded with micrometric tracer particles whose motion is recorded by two Phantom v12 high-speed cameras which are connected to an external trigger that syncronizes them. They are in turn connected to the computer which controls the acquisition sequences (this will be detailed in the following sections). The cameras can reach up to 6242 fps and have a 12-bit, 1280 × 800px resolution.

A schematic representation is shown in figures 4.13 a) and b).

Hydrogeal beads

In order to look inside the bed we need to have "transparent" solid beads. That is, we need the beads to match the refractive index of the fluid. Refractive index matching consists on pairing a solid material and its carrier fluid, each with a similar index of refraction, so that the solid is virtually transparent.

Our particles of choice were hydrogel beads, which are superabsorbent materials that are expanded 200 times their initial weight by absorbing water, and are thus made up of ∼ 99% water, so their refractive indices are almost the same, with studies reporting a 0.1% difference with water [START_REF] Byron | Refractive-index-matched hydrogel materials for measuring flow-structure interactions[END_REF]Dijksman et al., 2012;[START_REF] Harshani | Experimental study of porous media flow using hydro-gel beads and LED based PIV[END_REF]. This can be observed qualitatively in figure 4.14, where the"Université de Lyon" badge is perfectly readable when the hydrogel beads are submerged in water. This choice of hydrogel beads was mainly motivated by the fact that there is no need to change the fluid (in our case, water) to index match the particles. If we were to use other materials they would be very sensitive to the mixtures of the fluid needed to index match the solid, and these can change their optical properties with time (e.g., due to changes in temperature) and this would add an additional experimental constraint [START_REF] Wiederseiner | Refractive-index and density matching in concentrated particle suspensions: A review[END_REF]. Most of all, other choices are either with viscous fluids (e.g. PMMA and glycerol) or not so simple to handle (ammonium thiocyanate which is corrosive or DMSO). Hydrogel beads have been used in previous Particle Imaging Velocimetry (PIV) studies, which further encouraged us to use commercially available hydrogels [START_REF] Harshani | Experimental study of porous media flow using hydro-gel beads and LED based PIV[END_REF]Weitzman et al., 2014;[START_REF] Byron | Refractive-index-matched hydrogel materials for measuring flow-structure interactions[END_REF][START_REF] Jensen | Sap flow and sugar transport in plants[END_REF].

A drawback of hydrogel beads however is that they are fragile and easily deformable when a mechanical force is applied to them, so before each experimental run we checked that they maintained their their spherical shape and that they were not broken, this can be seen in figures 4.17, 4.19 and 4.20. This visual appreciation is good enough because the shape of the spheres is not of the outmost importance in these experiments, but the pores are. The beads used have a mean diameter d = 14mm after being submerged in water for at least 36h.

Moreover, the removal of the bubbles was not done with the same protocol as in the previous setup (i.e., heating the water to 55℃), but by letting the water flow through the circuit for as long as it took without changing its temperature (usually about 16h). This was done so as not to change the properties of both the water (e.g., its viscosity) and the beads. Bottom: the container with hydrogel beads is filled with water. The beads index-match the water, as the badge can be clearly read even when they are submerged in water.

The water flowed at the lowest flow rate so as not to deform the beads. After letting the water circulate, we visually corroborated that there were no bubbles in the system, which is easy because the bed is transparent.

PTV Calibration

Because we are going to work with data obtained by digital cameras, the positions obtained by them is represented in pixels. It is then very important to have a very neat and correct calibration so that we can have the right transformation from digital units into physical ones (px→mm). Not only that, but we also need a transformation that allows us to transform the 2D images into the 3D space. This is done with the method developed in [START_REF] Machicoane | A simplified and versatile calibration method for multi-camera optical systems in 3d particle imaging[END_REF], where the calibration is done in at least two different 2D planes at different distances from the cameras. Each camera only requires two (or more) plane-to-plane transformations to be determined between the sensor pixel plane and two (or more) planes in the real world, by linear or cubic transformations. The target shown in figure 4.15 is used to calibrate the system. It has two different planes in depth, each consisting of 72 equidistant points separated by 5mm. on the PTV data.

Particle Tracking Velocimetry

Particle Tracking Velocimetry (PTV) consists in an optical lagrangian technique where tracer particles are tracked with fast cameras, and once their trajectories are acquired their instantaneous velocity fields and other hydrodynamic magnitudes of interest can be computed.

The tracer particles have to be small enough so that the smallest scales of the system can be resolved and they have to be neutrally buoyant so as to be treated as fluid elements. Their three-dimensional position can be measured by two -or more-fast cameras that record the region of interest frame by frame.

The algorithm for tracking the particles was developed in MATLAB and adapted from previous codes used in the Turbulence group at ENS Lyon (Dumont, 2021; Laplace, 2022), and was done with the help of Alexandre Ponomarenko, a postdoc candidate at the laboratory at the time. The PTV is done in four different steps which will be detailed below. x x This step consists of removing the background noise by setting a threshold value and by identifying the tracer particles from their intensity and size. As a first guess possible particle candidates at pixels of higher intensities. Once those point-like particles are identified, their sub-pixel position is estimated by a two-dimensional gaussian, with its center (its maximum) representing the particles' positions (Ouellette et al., 2006a). From this, the position in each frame is saved and we will be working these hereafter. two raw time-superposed images taken with both cameras in a generic experiment, where the centers/particles are identified as specified before.

Tracking

Once the centers are found, the particle trajectories are reconstructed for each of the cameras. This is achieved using a predictive tracking algorithm similar to the 4-frames best-estimate algorithm in Ouellette et al. (2006a).

Matching

In a first step, trajectories in pixels are projected onto the "real world" (RW) coordinates, i.e., physical coordinates, by using the calibration transformations previously described in section 4.2.1, momentarily neglecting the 3D nature of the actual tracks. This preliminary step allows us to match corresponding trajectories seen by both cameras. Once the trajectories for each camera are both in the same coordinate system (in mm) we identify the ones that were measured simultaneously by both cameras independently. That is, we have to match the trajectories. Because each camera measures a slightly different region, the same trajectory is going to be measured at a slightly different relative position that varies, so in order to match two trajectories it is necessary to specify the maximum distance δ that there can be between the same trajectory measured by both cameras. That is to say that the trajectories X 1 for camera 1 and X 2 for camera 2 have to be the same, and we impose that ||X 1 -X 2 || < δ ∀t. The last thing that we set is that both trajectories last for the same number of frames (i.e., the same time).

Ray crossing

Lastly to retrieve the actual 3D trajectories, the "rays" are crossed. By "rays" we mean that every pixel of the camera sensor can be associated to an imaginary ray so that any particle along that line will be imaged on the same pixel. This line can be obtained by projecting each pixel coordinates for each particle detected on each frame of each camera to the two calibration planes. The intersection of the lines corresponding to matched trajectories then allows to retrieve the 3D tracks in the real world, as schematized in figure 4.18. We then have to cross the rays from each camera and their intersection is going to provide us with the 3D positions of the particles. This is done by using the data obtained from the calibration, that transforms px→mm. In order to choose the tracer particles and to choose our region of interest, a first visualization of the bed was done by bathing the hydrogel beads in a fluorescent dye (figure 4.19). We then illuminated the bed with the laser and took a picture. It is worth noting though that this was not a packed bed, as there was no grid on top. This was done to have a rough estimate of the size of the pores, which is of the order of 1-3.5mm, and we expect this value to decrease when the bed is fully contracted.

The measured size of the pores leads us to use fluorescent polyethylene particles as tracers, which have diameters ranging between 27 to 35µm (∼ 100× smaller than the pores). They have a 0.996g/cc density, so they can be treated as neutrally buoyant. Because of this, they can be treated as elements of the fluid.

Acquisition

The cameras are piloted by National Instruments' software LabVIEW, where we can set the cameras parameters (acquisition frequency, exposure time, spacial resolution, etc), and we can control the translating stage. Several things have to be kept in mind before starting each experiment.

On the pictures acquired by the camera, each tracer particle has to be at least 3px in diameter to be properly detected by the sub-pixel gaussian fit center finding algorithm. If they were less than 3px we would have a pixel-locking (or peak-locking) effect, which is the tendency for the measured location of a particle image to be biased towards integer values (Michaelis et al., 2016). They also have to have an intensity high enough to be detected by the code, and all values below a certain intensity will be ignored and not treated as tracer particles. One of the main aspects that needs to be taken care of is the sampling frequency, or the frames per second fps, that needs to be set, and it is determined by the spatial resolution of the cameras and the tracers velocity.

Videos need to be recorded in such a way that we can be sure that we are tracking the same particle over two successive frames. For example, if we have two nearby particles in a frame, and in the following one they moved 10px each, a priori there is no way to know which particle is which. On the other hand, if for example there are 1 or 2 pixels per particle between two successive frames, it is easier to determine which particle is which at the next frame.

With that being said, we can roughly estimate the number of fps needed. To choose it, we need to estimate the maximum velocity that we would expect in our setup, which in our case would be a the velocity in the pores, u p . We can estimate this by saying u p ∼ U/ε, with U the superficial velocity and ε ≈ 0.4.

Next, we will need the spatial resolution of the cameras a, which has units of L/px. Noting f ≡ 1/dt the framerate, we can compute the particle displacement Disp between two successive frames as Disp = u p dt/a, where dt is the time between frames. Replacing this into the framerate, we get

f = 1/dt = u p a/Dist
With all of this taken into account, the f ps was set at 2200Hz, with a space resolution of 896 × 800px, and we recorded 3 seconds (a total of 6600 frames) for each acquisition. This was done several times for several planes of the bed, and for several Re.

The coordinate system is defined as (x, y) on the plane of acquisition, with y the flow direction, and z the other direction perpendicular to the flow rate, which is the depth of the bed. This is shown in figures 4.20 and 4.13.

Measurements are done in the course of hours (between 3 to 6 hours for each run) at night when there is less light pollution. The laser sheet is moved in the z-direction so that we have measurements in different depths of the bed, and at different times. We then compute the the velocity fields for all the data collected. This is possible because the system we are studying is stationary, so it does not change with time. This was verified by looking at the trajectories (there are certain ones that easy to identify, such as the ones on recirculation points), and by verifying that the dynamics do not change with time.

Once the measurements are taken, the data is saved. Each measurement for each camera takes up to 6Go of space, and we can save up to 1To of data for each experiment. The raw data is then processed by the PTV algorithm and is then analyzed as detailed below.

U

PTV Errors

As illustrated in figure 4.21a), whenever the cameras' rays are crossed, there is an error committed between the point found by their intersection and the real point in the physical space. As a matter of fact, rays never exactly cross, both because of uncertainties in centerfinding precision and calibration accuracy. These uncertainties are more pronounced here than in classical PTV in homogeneous fluids due to slightly imperfections of the index-matching. We define the error as the minimal distance between the two matching rays, ||r ray 1r ray 2 ||. Its probability density function (estimated over all the detected and matched particles) is shown in figure 4.21b) for an arbitrary experiment. Its mean value, which we take to be the mean error is 0.1mm for all of the experiments, which is approximately 20 times smaller than the typical pore size (see figure 4.19 for reference). This determines the raw 3D resolution of our PTV data and the level of noise. Histogram of the error, defined as the distance between the two rays in the intersecting points.

Postprocessing

In order to smooth the trajectories and compare the signal-to-noise ratio on velocity estimates, the trajectories obtained from the PTV algorithm are fitted by polynomials of order 3, which as can be observed in figure 4.22a), it is a good approximation, as trajectories are short. The mean difference between the "real" trajectories and the fitted polynomial d = (X real -X f it ) 2 depends on each coordinate, and is

d x = (x real -x f it ) 2 = 8e-3mm, d y = (y real -y f it ) 2 = 9e-3mm, and
d z = (z real -z f it ) 2 = 6e-2mm
The error in the z-direction is higher because it is highly sensitive to the calibration in the z-direction, which has a larger associated error. Indeed, considering the optical arrangement of the cameras, the x and y coordinates are highly redundant from each 2D view, while z is fully retrived by . The probability density function of the errors in the trajectories is shown in figure 4.22b), and we choose to work with the trajectories that have an error d i < 20µm in the x and y directions. From the fitted trajectories we can then compute the velocity and acceleration fields by differentiating the positions. As a result of the added temporal coherence along trajectories, the error from fitting the trajectories is about 10 times smaller than the error committed from individual ray crossing.

Validation

Before passing onto the "real" experiments and in order to validate global workflow of the code, we carried out an experiment with no hydrogels, and injected tracer particles into the flowing water at different Re. We expect the tracer particles to follow trajectories parallel to the flow and to have the same velocity as the one measured by the flowmeter because as there is no porosity u = U . The difference between the velocity measured with the PTV algorithm and the one measured with the flowmeter was of 0.1%, thus validating the workflow described above.

Experiments meet simulations

The experiments allow us to study the system at a macro and a micro scale. From the first set of experiments we measure the pressure and porosity of the bed, both global quantities, to compute the presure drop and see how it variates with Re and D/d. From the second we measure the velocity field and acceleration fields, and as it will be detailed in chapter 6, we will compute correlations and structure functions.

The experiments will not only be used for the "calibration" of the numerical simulations, but in fact the experiments and the simulations will complement each other. Whereas the experiments will provide us the real-life physics, the simulations will allow us to go even further. For example, in the simulations we will have access to the local pressure field, which is not available experimentally with the current setups. They also allow us to study the forces involved in the system, which we can not measure experimentally.

With this interplay between experiments and simulations we will have a thorough description of the beds in different scales, which can be extrapolated to other kinds of porous media. Chapter 5

Global confinement effects

This chapter is presented in the form of a paper that is going to be submitted between January and February 2023. This study is oriented towards macro-scale hydrodynamics in fixed beds, and it presents both numerical and experimental results (the experiments were conducted in the first set-up presented in chapter 4). In order to quantify the border effects, we explore three different configurations with D/d = [5, 8, 10], fixed porosity and we explore a Renynolds-numer range that varies between 200 and 1000 in the experiments and 60 and 300 in the simulations/

As it has been mentioned in the introduction, there have been plenty of studies that evidence how the confinement effects of packed bed is reflected in, for example, the mean porosity of the medium, but we are taking a different approach than the ones found in the literature. The study is mainly divided into two sections: the first one is dedicated to the variability of the pressure field and pressure gradient, and how they depend on the random arrangement of the particles. In the second part we explore the physical mechanisms responsible for how the pressure drop changes with D/d.

We first want to explore the wall effects from a different perspective: how the variability of the results is affected by the presence of the walls. This is a valid question, given the fact that the arrangement of particles is random and there must by an intrinsic variability linked to that. We study this experimentally by computing the pressure drop for 9 different configurations (the only variable that changed was the distribution of the spheres) and calculate how much the results vary from their mean (i.e., their standard deviation). We observe both a variability in the pressure curves and the pressure drop obtained for each configuration.

The other question that we want to answer is if the measurement depends on the placement of the captors at the walls of the bed: is it the same if we measure at θ 1 or θ 2 (two generic angles)? We explore this numerically by obtaining the pressure along a vertical line at different angles of the bed, and in order to be as close as possible to the experimental measurement, we do this at the wall of the reactor. Effectively, what we observe is that the pressure varies with the placements of the "sensors" (our numerical probing lines) and that in fact this variability can be associated to the one observed for different ensembles: it seems that measuring at, say, 9 different angles of the bed is comparable to measuring 9 different ensembles with the same D/d, ε and Re. This is an interesting result, given that we could do a sort of ensemble averaging by using a single numerical bed and measuring at different angles. Moreover, this variability increases with D/d, making the wall effects more evident.

The last question we ask ourselves in this first part is whether it is the same to measure the pressure drop in the walls (as it is usually done in experiments) to calculating the averaged pressure gradient over the whole volume. This is of interest because while experiments are usually done in such a way that the sensors are placed at the walls of the reactor, the theoretical developments are usually done by averaging over the fluid volume.

In the second part of the article we present how the pressure gradient varies with D/d and Re both experimentally and numerically, and we recover results that are consistent with previous studies. We then try to quantify the influence of the walls in the pressure drop with the aid of the simulations. We do this by calculating the total particle-fluid force of all of the spheres and the wall. We observe that in fact, the pressure drop is not affected by the walls per se, but by how the walls affect the particle distribution. Using experiments and numerical simulations, this study examines the intrinsic variability in confined fixed beds of randomly arranged spherical particles, specifically in relation to the pressure field and pressure gradient within the bed when there are wall effects. We observe that as the sphere-to-reactor diameter ratio decreases, the pressure gradient presents a stronger variability, particularly in relation to where measurements are taken within the reactor. The study also quantifies the difference between measuring the pressure drop at the wall of the reactor 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 walls of the reactor, finding that the pressure drop follows a consistent 1/Re scaling regardless of the confinement due to changes in the mean porosity of the bed. We also find that the vast majority of the pressure gradient is due to the distribution of the solid spheres within the bed, rather than the walls themselves.

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); [START_REF] Elouali | Physical models for packed bed: Sensible heat storage systems[END_REF]; [START_REF] Barker | Heat transfer in packed beds[END_REF]. 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 and local porosity variation near the reactor wall Mueller ( 2019 For a fixed bed of spherical particles the parameters involved are the Reynolds number Re = U d/ν based on the superficial velocity (i.e., the velocity as if the spheres were not present), the bead diameter d 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 f luid /V total ; and D/d, which is the ratio between the particle and reactor diameter, d and D respectively. When border effects are negligible, i.e. when D/d → ∞, the pressure drop depends only on Re and ε. In this case, the Ergun model [START_REF] Ergun | Fluid Flow Through Packed Columns[END_REF] is generally used to describe the pressure drop of random arrangements of spherical fixed particles. Nevertheless, this model is not suitable for cases where the border effects are evident, particularly when D/d < 10 Eisfeld & Schnitzlein (2001); [START_REF] Erdim | A revisit of pressure drop-flow rate correlations for packed beds of spheres[END_REF]; Hill et al. (2001); Clavier et al. (2015); [START_REF] Foumeny | Correlations of pressure drop in packed beds taking into account the effect of confining wall[END_REF]. It has also been reported that, in the laminar regime, 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 & Schnitzlein (2001).

In Eisfeld and Schnitzlein Eisfeld & Schnitzlein (2001), the variability of the coefficients of equation (5.3) was studied by analysing more than 2300 data points with different border 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 & Gibilaro (2004) where the authors define a bulk velocity in terms of D/d. In De Klerk De Klerk (2003) and Zou et al. Zou & 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 more or less 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? 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. [START_REF] Erdim | A revisit of pressure drop-flow rate correlations for packed beds of spheres[END_REF]), 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.

The questions that we want to address in this work are: Do the Ergun scaling laws are valid when D/d < 10? If so, what are the physical mechanisms involved? Does the randomness of the arrangements affect the results? In order to answer all these questions, pressure drops in fixed beds are studied both experimentally and numerically, for three different configurations in the inertial regime. The simulations will 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 calculating it like it is typically done in the experiments. The results are compared with those obtained in Mueller (1992) and Benenati & Brosilow (1962).

Computational approach

The same problem is addressed numerically in 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 software Grains3D Wachs et al. (2012). The fixed beds are done in three successive steps: first, the particles are dropped into a nonperiodic 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 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 so 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 in figure 1, where the results are consistent with those obtained by Mueller (1992); Benenati & 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 Weller et al. (1998). For each bed, the flow domain in between the particles is meshed with the OpenFOAM F. Falkinhoff, J.L Pierson, L. Gamet, M. Bourgoin and R. Volk 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 maximum 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 analyzing 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 [START_REF] Delaney | Defining random loose packing for nonspherical grains[END_REF]; De Klerk (2003). Moreover, the system was made periodic along the main flow direction (zdirection), so as to avoid any inlet/outlet boundary conditions effects and to compare it to the experimental measurements. The number of cells ranges between 45 to 70 millions, depending on the case, the axial periodicity being 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]. 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 ∇ • u = 0, (2.1)

ρu • ∇u = -∇p + µ∇ 2 u, (2.2)
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 equation ( 2.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 equations (2.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 differencing for the interpolation.

In order to impose a mean flow with periodic boundary conditions, a forcing term is added to the incompressible momentum equation such as

ρu • ∇u = -∇p + µ∇ 2 u -fẑ, (2.3)
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 simple relation

p = p + f z, (2.4)
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 boundaries (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 summed over the entire flow domain.

We also used a residual control, down to 10 -6 for both the velocity and pressure fields. A similar condition was also used in [START_REF] Magnico | Hydrodynamic and transport properties of packed beds in small tube-tosphere diameter ratio: pore scale simulation using an eulerian and a lagrangian approach[END_REF]. 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 computations performed by Hill et al. (2001) with a Lattice Boltzmann method as a benchmark case. Ordered arrangements of spheres are computed by placing a sphere in the middle of the three-periodic unit cell (i.e., a single-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 figure 2, this quantity was calculated and compared to the one obtained by Hill et al. (2001). 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 not larger than 2%.

F * = F 6πµRU ,
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 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 solver 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.

Experimental set up

The pressure drop measurements are conducted in the experimental device shown schematically in figure 3. 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 plexiglas column with a diameter D = 0.04m, 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 cm high. Two additional 0.2 m empty sections are added so as to separate the bed from the nozzle and a 0.05 m high honeycomb is placed so as to suppress velocity fluctuation in the incoming flow.

The test section is first filled to its full height with monodispersed stainless steel 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 ̸ = 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 dp dz , is then calculated by computing the slope of the p(z) curve as shown in figure 4 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.

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 figure 5 which displays the normalized pressure, p(z)/ρU 2 , measured experimentally at the surface along the tube for 9 independent 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 ∼ O(10).

Investigation of the local pressure field is impossible in the experiments, so we shall 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 = cte. It reads:

p ′ (r, θ, z) = p -2/R 2 R 0 p(r, θ, z)rdr = p(r, θ, z) -⟨p⟩ r,θ (z), (4.1)
where ⟨•⟩ r,θ denotes a spatial average over variables r and θ while ⟨•⟩ θ is an average over θ only and ⟨•⟩ V an average over the whole volume. Figure 6 shows the pressure fluctuations field p ′ /ρU 2 , at two planes at two different heights of 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 F. Falkinhoff, J.L Pierson, L. Gamet, M. Bourgoin and R. Volk 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 a certain 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 figure 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 = 1 ρU 2
dp dz that we can plot against 1/Re.

The results are displayed in figure 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 toward a constant at high Re which corresponds to the inertial regime beyond the purely viscous (Darcy) regime.

As can be seen in figure 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 values of α and β, drawn as a solid line in figure 7. We then obtain the standard deviation from the ensemble average,

σ G p = G p 1 Re -G p 1 Re 2 1 2
, and find that the standard deviation averaged over all the values of the Reynolds number is σ 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 (figure 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 figure 8a) for a particular case of D/d = 10 at Re = 200. Each point is calculated by averaging the pressure over θ ± dθ at different z, which is analogous to having finite-sized pressure sensors at the wall (i.e. they span an angular length at the wall of the reactor θ ± dθ) at different heights: ⟨p(r = R, θ, z)⟩ θ (z). 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, θ + dθ, z)⟩ is shown in figure 8b) for 4 different values of the Reynolds numbers. This was calculated for different θ ± dθ by a linear fit of the curves as is done experimentally. 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, it is striking that 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 organisation 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 figure 9. This was calculated by computing the standard deviation of the superficial pressure drop described above and its mean value. 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: σ G p (r=R,θ,z) /⟨G p (r = R, θ, z)⟩ θ 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 O(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 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. Nevertheless, this variation is found of the order of 1%, which is less than the error that one commits when measuring different beds or at a particular angle.

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.

F. Falkinhoff, J.L Pierson, L. Gamet, M. Bourgoin and R. Volk

The pressure loss is a function of 7 physical parameters so that 8 dimensional quantities are involved:

∆p, ρ, U, µ, 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 = ∆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 ≡ d∆p/L ρU 2 , Re, ε, D/d, L/D,
that can be written as:

G p = f (Re, ε, D/d, L/D).
(5.1)

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, ε, D/d).
(5.2)

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 (figure 5). We verified that L = H was not involved as well in the DNS by checking that G p remains the same when doubling the height of the bed.

The functional form of equation 5.2 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 [START_REF] Idelchik | Handbook of hydraulic resistance, 2nd edition[END_REF]. 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 = α(D/d, ε) 1 Re (Re ≪ 1),
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 [START_REF] Idelchik | Handbook of hydraulic resistance, 2nd edition[END_REF], which allows us to write

G p = β(D/d, ε) (Re ≫ 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 = α/Re + β. There is no a priori reason for which α and β should not depend on the Reynolds number in the inertial regime Re = O(100) when the scale separations is small, as there could be boundary layer or transitional effects. 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 [START_REF] Ergun | Fluid Flow Through Packed Columns[END_REF], which reads:

G p = (1 -ε) ε 3 A 1 -ε Re + B .
(5.

3)

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 [START_REF] Ergun | Fluid Flow Through Packed Columns[END_REF]. This equation has been observed to be valid for D/d > 20 [START_REF] Erdim | A revisit of pressure drop-flow rate correlations for packed beds of spheres[END_REF], and it 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 , and the numerical one is presented as the gradient averaged over the wall surface, that is, the average of the gradient shown in figure 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 cases follow a linear evolution of G p with 1/Re: G p = α/Re + β which 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 → ∞. 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, and thus generates a lower pressure gradient for lower scale separations.

The slope and intercept of the linear fits are shown in table 3, multiplied by a factor F. Falkinhoff, J.L Pierson, L. Gamet, M. Bourgoin and R. Volk depending on the porosity so as to have a direct comparison with A and B from equation (5.3), noted as A ergun and B ergun in the table.

D/d A Ergun ε 3 (1-ε) α exp ε 3 (1-ε) α dns B Ergun ε 3 (1-ε) β exp ε 3 (1-ε) β
Both the numerical and experimental values are within the values reported in Eisfeld & 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 previous works De Klerk (2003); Zou & Yu (1995) both variables depend on the other, so their contributions might not be separable. That is to say, that one might not be able to do, e.g.

α(ε, D/d) = α ε (ε) × α 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 figure 12 Gibilaro (2004). Whether this due to some friction at the outer boundary or because of 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. It 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 calculated 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 (5.4) This means that the pressure gradient is a result of the different contributions of the forces of the system, and we can separate those two contributions in order to see which one is more influential. Equation (5.4) 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 (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 F. Falkinhoff, J.L Pierson, L. Gamet, M. Bourgoin and R. Volk is the same normalization factor used in G p . As was mentioned, the total forces by unit volume are 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 W alls , and F T otal = 1.05 × F Spheres . Not only that, but both F W alls and F Spheres follow the same 1/Re scaling, which validates the robustness of the scaling (the walls contribution was replotted 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.

⟨∇p⟩ V = 1 V F Spheres + F Walls .
Figure 13b) shows F T otal 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).

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. That is, 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.

Conclusions

Fixed beds of randomly arranged particles present an intrinsic variability that is linked to the arrangements themselves, and 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 the scale separations get weaker (i.e., D/d decreases) the pressure gradient presents a stronger variability depending on where we measure 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. This is of particular interest because the theoretical developments are usually derived for bulk-averaged averaged quantities. We observed that there is a small difference between measuring at the wall of the reactor and averaging over the whole volume, of 2.5% at the most.

Our study showed that the variability has an intrinsic error that can reach 10% when the scale separations are smaller, and this can be of interest for models regarding the pressure drop in confined beds, as there can be two separate models that account for a 10% difference that might be inevitable, as it is intrinsically related to the wall effects on the packing.

On the other hand, we studied how the mean pressure gradient is affected when the walls play an important role. As a first observation, all pressure drops follow the same 1/Re scaling, which is consistent with the one proposed by the celebrated Ergun correlation. It is also noted that they decrease with D/d, given that the average porosity of the bed changes with D/d, and thus the pressure drop is affected. We quantified the contribution of the walls and the solid spheres by doing a simple force balance of the system, which satisfies that the pressure drop of the system is almost 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 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, whose force contributes to the pressure gradient a small percentage that is within the variability error of the measurement.

Chapter 6

Micro-scale hydrodynamics: a turbulent perspective Figure 6.1 -Time-lapse of a generic experiment at a generic plane of the bed. Even though we will be dealing with flows that are in the inertial or transitional regime, the Lagrangian dynamics is not trivial: we can observe different behaviors in different pores, including vortices in some pores and more ordered flows in other pores. This behavior is certainly caused by the randomness of the arrangement.

Experimental results

We will study the local multi-scale hydrodynamics of a fixed bed made up of hydrogel beads of diameter D h = 14mm with a particle to bed diameter ratio D/D h ≈ 6.43, by using indexmatching PTV as detailed in chapter 4. Its porosity is estimated to be ε ≈ 0.4 (De Klerk, 2003). As it has been explained in previous chapters, the confinement of the bed generated by the walls can have an effect on global quantities such as pressure. Here, we will explore whether the wall effects and the geometry are significant at the pore-scale.

As explained in chapter 4, we measure at a fixed z-plane of 4mm width, x being the other transversal direction and y the direction of the flow. We study the flow at four different Reynolds numbers based on the hydrogel beads' mean diameter and the superficial velocity 124,169,203,211], which is a regime that ranges between the inertial and the unsteady laminar (transitional) regime [START_REF] Wood | Modeling Turbulent Flows in Porous Media[END_REF]. Figure 6.1 shows a time-lapse image of a generic experiment in this regime, with the tracer particles that behave as elements of fluid shown in white. It can be qualitatively observed that the Lagrangian dynamics are not trivial: there are some re-circulation zones and the random arrangement generates pores where there can be different behaviors such as recirculation and stagnation points, channellike pores and so on (see, for example [START_REF] Patil | Turbulent flow characteristics in a randomly packed porous bed based on particle image velocimetry measurements[END_REF]). The flow shows a clear spatial complexity, with shear zones and strain, and rotation dominated structures, which recall the characteristic multi-scale structures of a turbulent flow. However, one should keep in mind that these structures are steady and stable in time. One of the goals of this chapter is to explore to which extent the spatial complexity of this steady laminar flow does resemble to turbulence. This will be done by deploying the same arsenal of statistical diagnoses to characterize the spatial fluctuations of velocity and their correlations as used for the study of turbulence. It can be seen that the Reynolds number based on the measured local fluid velocity u y can be several times higher than the "macroscopic" Reynolds number Re(U Q ) based on the superficial velocity U Q . This evidences the fact that not only there are very high local values of velocity, but this also probably generates high values of local acceleration as well, corresponding to more complex hydrodynamics at the pore-scale compared to the global scale of the bed. With all of this taken into account, we carried out a statistical analysis of the velocity and acceleration fields: one-point statistics including probability density functions and twopoint statistics including two-point correlations and structure functions. This will all be detailed in the following sections. Because of the symmetry of the system, we assume that the two perpendicular components x and z are equivalent. The system can thus be described by studying only one of the components. This will be further explored with the numerical results. As it has been mentioned, the velocity varies quite significantly within the bed. In order to quantify this variability, we first calculate the velocity component (transversal u x and streamwise u y ) probability density function (pdf), shown in figure 6.3. Gaussian distributions are plotted in red-dashed lines for comparison, as they are typically observed in fully-developed turbulent flows (Ouellette et al., 2006b).

U Q : Re = [

One-point statistics

The u x distribution is symmetrical around its mean value. It is quite remarkable that it is not that far off from a Gaussian, even though it has been established that these are not actual turbulent flows. On the other hand, the stream-wise pdf is asymmetrical around its mean value, with lower values of probability for negative velocities when compared to the Gaussian distribution. This is consistent with the fact that there is a fixed direction of the velocity, which is driven by the pressure drop that has been thoroughly studied in the previous chapter. The difference between both distributions shows that the flow is not perfectly isotropic in the x and y directions, at least at the scales that are set by the PTV measurement window, which approximately spans about (4 × 4)D h . This will be further explored later with the numerical results.

The acceleration pdfs are shown next in figure 6.4, and compared once again to a Gaussian distribution. The first observation is that they are remarkably similar to the pdfs observed in fully turbulent flows [START_REF] Porta | Fluid particle accelerations in fully developed turbulence[END_REF][START_REF] Voth | Measurement of particle accelerations in fully developed turbulence[END_REF], where extreme values of acceleration can be reached, as it is evidenced by the comparison with the Gaussian pdf and its stretched tails. This might point towards a kind of turbulent "signature" that is present in the system, which would be consistent with the fact that the tortuosity of the flow generates regions of relatively high and low porosity values and therefore the local velocity values change as well, generating extreme values of local acceleration. We have already stated that our flow is steady and is not actually turbulent. Our observation and that in [START_REF] Holzner | Intermittent lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF] at lower Re show that a turbulent-like pdf of acceleration can be retrieved without a complex temporal dynamics of the global flow, but just with a rich Lagrangian dynamics in a complex steady flow. Overall, the one-point statistics show a similarity between the velocity and acceleration probability density functions and those found in turbulent flows. The qualitative observations also show the presence of recirculation points and changes in local velocities. At this point, it is therefore tempting to push forward the comparison and explore two-point statistics seeking for possible multi-scale correlations, which are one of the most emblematic signatures of turbulence. Is it possible that the spatial structures of the flow present some turbulent signature in its correlations? In order to answer this question, we pass on to a multi-scale analysis, borrowing tools from the turbulence formalism, such as structure function analysis to identify the presence of an inertial range, an energy cascade and the formalism developed by Kolmogorov in 1941 (K41 theory). See, for example chapter 6 of reference [START_REF] Pope | Turbulent Flows[END_REF] or reference [START_REF] Davidson | Turbulence in Rotating, Stratified and Electrically Conducting Fluids[END_REF] chapter 9 for the multi-scale statistics in turbulence.

Two-points statistics

In the following sections we will be working with the velocity fluctuations, defined as u i = u iu i , as we are interested in the small-scale hydrodynamics and the role of the fluctuations at this scale.

Correlation functions and large scales Figure 6.5 -Left: Two-point auto-correlations for u x and u y . The velocity fields are correlated at the smallest distances before becoming null. Afterwards, they tend towards an oscillating pattern, which reflects the geometry of the bed. Right: Crossed auto-correlation function. The two components of the velocity are never fully correlated (that is, they never reach a value of 1).

In order to observe how two elements of fluid are correlated, we calculated the autocorrelation function of the x-and y-components of the velocity fluctuations u i = u iu i :

C ij (dr) = u i (r + dr) u j (r) i = x, y, (6.1)
with dr the relative distance between two elements of fluid (i.e., two tracer particles) in an r-sphere: dr = ∆x2 + ∆y 2 + ∆z 2 , and the mean is done for all the separations or increments1 . The averaging • is done over an ensemble of pairs of points with a given separation dr, and evenly sampling the entire flow. The results of C ii (dr) are presented in figure 6.5 (Left). It is observed that they both lose correlation at dr 0.2D h , and then continue to oscillate around zero for higher separations, which is due to the presence of the hydrogel beads. Considering that a typical pore size ranges between 1 and 5 millimeters 2 , this reveals the presence of correlation which develop typically at the pore-scale of the bed, and it shows that the velocity fluctuations are strongly correlated at this scale. On the other hand, the crossed-correlation C xy oscillates around zero for all separations, but it reaches a value of approximately 0.3 at the smallest scales. This is evidence of the non-gaussianity of the velocity distributions (see for example Ex.3.19 of Pope (2000)); this is shown in figure 6 The correlation length is calculated when this quantity reaches a plateau. Bottom: Correlation lengths for the transversal and axial components. They vary with the Reynolds number and seem to converge towards a constant value so as to reach a certain level of isotropy.

In order to see if there is a characteristic correlation length involved, we calculate the area under the curve of equation (6.1):

L i = dr 0 C ii (dr) σ u i dr i = x, y, (6.2)
which typically reaches a plateau, defining what is known in turbulence as the correlation or integral length L (Pope (2000),section 6.3). The results are shown in figure 6.6 (a, where the curve oscillate around a constant value, which is due to the geometry of the bed. The correlation lengths are calculated at their maxima (which are where C ii = 0) and are shown in figure 6.6 (b). It is interesting to note that L y increases monotonously with Re whereas it is the opposite for L x . Moreover, their ratio tends to one (figure 6.6 (c)). This points toward the isotropy of the system: whereas for lower Re there seems to be a higher correlation length in the transversal (x) direction, it is the opposite for the higher-Re case, where they both tend towards the same value. This means that as the Re number increases the large-scale fluctuations of the flow become more isotropic. This is consistent with the asymptotic regime for highly turbulent flows in porous media as explained in [START_REF] Wood | Modeling Turbulent Flows in Porous Media[END_REF], although it is remarkable that this isotropy already appears at a lower Re at the pore-scale, which is reflected by the fact that L x /L y ≈ 1.

Inertial range of scales

In order to zoom-in into the multi-scale pore-scale hydrodynamics, we calculate the secondorder structure function S 2 , which contains the same information as the auto-correlation function with a focus on the smaller scales. We first calculate it in the r-sphere for the xand y-components of the velocity, so that S 2 is defined as

S 2 i (dr) = δu 2 i (dr), (6.3)
with δu i = u i (r + dr)u i (r) and dr = r 2r 1 = (x 2x 1 , y 2y 1 , z 2z 1 ) for two different points r 1 and r 2 .. This S 2 tends to 2σ 2 u i when the velocity field is no longer correlated. We will also focus on the longitudinal second-order structure function, defined as

S 2 = (δu • δr) 2 |dr| 2 . ( 6.4)
The second-order structure function is of particular interest in turbulence because in homogeneous and isotropic turbulence (HIT) it carries one of the most celebrated signatures of the multi-scale nature of turbulence and exhibits the following scaling [START_REF] Pope | Turbulent Flows[END_REF], eq. 6.30): (6.5) in the inertial regime, where the scales are much larger than the Kolmogorov length scale η (which corresponds to scales below which the dynamics is dissipated by molecular viscosity ν such that η = (ν/ )3/2 ) and much smaller than the integral length scale L. is the energy injection rate and C 2 is a universal constant, generally referred to as the Kolmogorov constant 3 . In HIT C 2 ≈ 2 [START_REF] Pope | Turbulent Flows[END_REF][START_REF] Sreenivasan | Multiphase flows: Rich physics, challenging theory, and big simulations[END_REF]. The component-wise secondorder structure functions also follow the same two-thirds scaling, although with different constants. Not only that, but it comes from the self-similar K41 theory that the structure function of order p follows:

S 2 (dr) = C 2 ( dr) 2/3 ,
S p (dr) = C p ( dr) p/3 . (

Figure 6.7 -a): Second-order structure functions for each velocity component as a function of the relative distance in a r-sphere, which depends on all three spatial coordinates. A two-thirds scaling is evident at the small scales dr < 0.2D h , which points towards the presence of an inertial range. b) Second-order longitudinal structure function, claculated with all three components of the velocity and of the relative distances. It shows once again an inertial range at the smaller scales.

It is known that K41 fails at predicting actual scalings for high-order structure functions (p ≥ 4) (Frisch, 1995). This reveals the intermittency phenomena and the fact that turbulence is not scale invariant -in particular small-scale energy dissipation is unevenly distributed in space. The dataset of the present work does not allow to reach the statistical convergence required to explore such high order statistics. We will therefore focus on low-order statistics (p 3) seeking for possible scalings in the present system, possibly reminiscent of K41 dynamics. If these scaling are maintained, we would have strong evidence that the spatial organization of these laminar flows shares similarities with a turbulent flow. Given its temporally steady nature, it would then be a practical example of actual "frozen turbulence", which is a classical approximation when analyzing wind tunnel experiments using the so-called Taylor hypothesis [START_REF] Pope | Turbulent Flows[END_REF].

Up until now we cannot fully link the previous results to a turbulent-like formalism, because we still have to identify an inertial regime, which would appear in the second-order structure functions. Figure 6.7a) shows S 2x and S 2y respectively. It shows a clear 2/3 spectrum that continues on until dr ≈ 0.2D h for all the Reynolds numbers studied. This points to the presence of an inertial range at scales dr < 0.2D h . The oscillating pattern for dr > 0.2D h is the long-scale correlation which oscillates because of the presence of the beads, and is in average 2σ 2 u i as expected. We can then verify that the uncorrelated plateau for S 2 where the inertial range ends, is reached at approximately L ≈ 0.2D h . Figure 6.7b) shows S 2 , which once again supports the presence of an inertial range at dr < 0.2D h and the large-scale decorrelation. We can eventually calculate the analogous of the Kolmogorov constant of the system, but we first need to determine the equivalent of , the turbulent energy dissipation rate. It is worth noting that at smaller scales (smaller than η) there should be a ballistic regime (that goes as dr 2 ) characteristic of the smooth dissipative dynamics of small scales and this should provide us the energy-dissipation rate, but we do not recover that result here. This might be because we do not fully resolve the finest length scales of the flow (the seeding density of tracer particles does not give access to statistically relevant information for separations below dr ∼ 10 -2 D h ), and also because the measurements become noisier at the smaller scales, given that it becomes difficult to match the particles within PTV errors (primarily related to imperfect index-matching effects).

Energy cascade

The presence of an inertial range does not suffice alone to have a full turbulent-like description, which is incomplete without calculating the energy cascade, if there is one. In order to do this, let us now turn to the third-order longitudinal structure function S 3 defined as

S 3 = (δu • δr) 3 |dr| 3 ,
which satisfies within the inertial range

S 3 = - 4 5 dr (6.7)
for a direct cascade. Contrary to S 2 , these inertial-scaling laws can be derived from the Navier-Stokes equations [START_REF] Lindborg | A note on kolmogorov's third-order structure-function law, the local isotropy hypothesis and the pressure-velocity correlation[END_REF]. From the physical point of view, while S 2 reflects the scale by scale kinetic energy distribution of the system, S 3 carries the information on how energy flows across scales. The fact that in 3D HIT S 3 < 0 reflects the direct nature of the energy cascade: energy is injected at large scales (∼ L) and flows through inertial scales down to the small scales (∼ η) where it is dissipated by viscosity. Some systems (such as 2D turbulence) are known to exhibit a range of scales with the opposite trend, energy flowing from the smaller scales to the bigger ones, and then S 3 > 0. The third-order structure function is then directly linked to the inertial-range energy cascade present in turbulent flows by its dr relation and sign. Unfortunately S 3 needs third-order velocity statistics which need more data to converge, so it was not possible to calculate it accurately. However, this can be overcome by considering a derived second order relation.

Indeed, for a locally homogeneous and isotropic flow [START_REF] Mann | Experimental study of relative, turbulent diffusion[END_REF][START_REF] Ott | An experimental investigation of the relative diffusion of particle pairs in three-dimensional turbulent flow[END_REF][START_REF] Hill | Opportunities for use of exact statistical equations[END_REF] δa

• δu = -2 , (6.8)
with the sign depending on whether there is a direct (< 0) or inverse (> 0) energy cascade. This quantity yields information about the energy cascade by having a constant value across scales in the inertial range and its sign, and does not involve any empirical unknown constant so it allows to estimate directly (as with S 3 ). It has a better convergence rate as it requires second-order statistics, although at the cost of a high-order derivative order, as it requires acceleration to be known besides of the velocity. The lagrangian framework of our diagnosis is particularly suited to accessing this acceleration. We calculated this quantity by only using the y and x components, because the z component is much noisier than the other two. Moreover, because of the symmetry of the system, the z component should behave as u x and a x , so we define δa • δu = (2δa x δu x + δa y δu y ). Figure 6.8 shows the results, where it is first observed that it remains relatively constant for the range of scales dr, where the inertial range is observed for S 2 . Not only that, but it remains negative in that range, especially in the higher-Re flows, which indicates the presence of a (local) direct energy cascade at the inertial range.

Figure 6.8 -The acceleration-velocity structure function, which is supposed to remain constant in a homogeneous and isotropic turbulent flow, and its negative sign indicates the direction of the energy cascade. In this case it remains negative for almost all distances at where the inertial range has been identified.

We can then calculate the value of the energy injection rate by averaging over dr < 0.2D h :

= | δa • δu | 2 . (6.9)
The values are shown in table 6.1 and vary between [2 × 10 -3 -2.7 × 10 -2 ]m 2 /s 3 .

The calculated energy dissipation rate reflects the energy transfer across the inertial scales, and it should be equal to the energy injection rate at the larger scales. Dimensionally in turbulence, the energy injection rate is related to the velocity standard deviation σ u by

= C σ 3 u L , (6.10)
where C is a constant and L is the correlation length. Not that in turbulence C is known not to be a universal constant and highly dependent on the geometrical properties and the forcing of the turbulence, in particular in the vicinity of production regions of turbulence (such as wakes) [START_REF] Vassilicos | Dissipation in turbulent flows[END_REF]. Figure 6.9 shows | δa • δu | as defined in the previous subsection, which in the framework of the present analogy with turbulence equates to 2 vs. σ 3 u . It shows reasonable agreement with a linear fit, with C ≈ 0.14, hence validating the scaling law (6.10).

Figure 6.9 -The energy dissipation rate (minus a factor of 2) as a function of σ 3 u . Is is observed that they follow a linear relation.

Re

[m 2 /s 3 ] η [µ m] τ η [s] C 2 124 (1.90± 0.50)e-03
151± 97 (2.29± 0.068)e-02 (1.20± 0.27)e-03 169 (6.70± 1.00)e-03 110± 33 (1.22± 0.015)e-02 (1.00 ± 0.13)e-03 203 (2.77± 0.14)e-02 77.5± 0.98 (6.00± 0.018)e-03 (5.92± 0.25)e-04 211 (2.71± 0.13)e-02 77.9± 9.34 (6.07± 0.17)e-03 (7.13± 0.27)e-04 Table 6.1 -Turbulent parameters calculated from the results obtained by calculating δa • δu = ∓2 , including the energy dissipation rate, the Kolmogorov length-and time-scales, and the Kolmogorov constant.

Full characterization of the turbulent-like dynamics

In order to have a full turbulent-like characterization, we can also calculate other turbulent parameters, such as the Kolmogorov length-and time-scales η and τ η respectively, defined as .11) The values are shown in table 6.1. η ranges between 150 and 78µm4 , and it decreases with Re, which indicates that the inertial range increases with increasing Re: compared to L = 2.8mm, we have that L η = [18. 6 -36] from the lower to the higher Re cases. τ η on the other hand ranges between [2.3e-2 -6e-3]s and decreases with increasing Re, which if we make the analogy with a turbulent flow, it would be the eddy-turnovertime at η.

η = ν 3 1/4 and τ η = ν 1/2 . ( 6
We can also calculate C 2 from equation (6.5). As shown in table 6.1, the values range between 10 -3 to 10 -4 , which not only it is not constant, but it is significantly lower than the value found in HIT, which is close to 2. This is evidence of the small "signature" that is present in the system. With all these quantities at hand we can calculate several Reynolds numbers based on different scales, and they are shown in table 6.2. We first define a "turbulent" Reynolds number that depends on the velocity standard deviation and the correlation length scale L = 0.2D h : .12) Because of equation ( 6.10), we can also write 6.13) which involves both the dissipation and correlation scales. We can also define a Reynolds number based on the Taylor microscale λ (see [START_REF] Pope | Turbulent Flows[END_REF], ch. 6), that by making use of equations (6.12), (6.10) and the definition of λ is: This R λ is found to range between 120 and 160 (see table 6.2). It is interesting to note that this is comparable to the turbulent Reynolds number of actual turbulent studies as performed in grid-turbulence wind tunnels [START_REF] Obligado | Preferential concentration of heavy particles in turbulence[END_REF][START_REF] Mora | Energy cascades in active-grid-generated turbulent flows[END_REF]Bourgoin et al., 2011). Even though we are dealing with a system that is not macroscopically or globally turbulent and more notably is a steady laminar flow (though spatially complex), at the pore-scale there is a global multi-scale dynamics that is not that spatially different from that found in turbulent wind tunnel experiments for example. This once again reflects the multi-scale nature of the system: even though at the global scale of the bed the flow is not turbulent, we find that the combined pore-scale hydrodynamics are compared to those of a turbulent flow, at least at the spatial sense. Simulations were computed for a bed with D/D h = 6.45 and height H/d = 4.84; both parameters are similar to those from the experiments. The grid is made up of 50.4 × 10 6 cells, and its resolution varies from 19 cells per bead diameter (based on the background grid) to 75 when the refinement is taken into account. For the sake of comparison we will define the y coordinate as the stream-wise direction and we will note the particle diameter as D h .

Re L = σ uy L ν . ( 6
Re L = C -1/3 L η 4/3 , (
R λ = C -1/2 15Re L . (6.14) Re Re L = σ u L/ν Re L = C -1/3 (L/η) 4/3 R λ = C -1/2 √

Comparison with numerical results: porosity effects

The Reynolds numbers based on the superficial velocity are Re = 160 and 211, and the porosity is ε = 0.49, which is the most different parameter from those of the experiments5 . We will use this to our advantage to study how the local hydrodynamics varies when the porosity changes. A snapshot of the velocity field of a simulation for Re = 160 is shown in figure 6.10, where the colorbar shows the Reynolds number based on the module of the local velocity |u|. As with the experiments, there are regions with very high local velocities when compared to the superficial one.

One-point statistics

The velocity probability density functions are shown in figure 6.11. The results are separated into the different components, as the simulations allow us to study the z-component of the velocity field without any measurement error. We can further compute the velocity pdf of the radial velocity, u r = u 2

x + u 2 z in order to study the transversal isotropy. As a first observation, unlike in the experiments, all the distributions have significant spikes near zero. This is very likely due to the fact that the simulations resolve all the scales, in particular the thin boundary layers near the beads, which are out of reach of the present experimental PTV (aiming at a multi-pore characterization of the flow) and which have a no-slip boundary condition. Regarding the transverse (x, z) and radial (r) components, they all share the same distribution and the same symmetry found with the experimental results. This validates our assumption that the the two transversal components are analogous, and it is further evidenced by the u r distribution. Moreover u x = u y = u r = 0 for both of the Reynolds numbers explores. This distribution differs further from a normal distribution than the experimental results, which is expected due to the marked peak around zero already discussed, which comes with a relative reduction of fluctuations resulting in sub-gaussian fluctuations in the ±2.5σ u range. As for the stream-wise component, the same skewness as in the experiments is observed, where negative values have once again a lower probability, and it is even lower than the one calculated on the experiments. This can be interpreted as a possible role of porosity effects: as the medium is less porous, recirculation points might become more rare, thus generating a mainly upwards flow.

The acceleration is calculated by the material derivative, because ∂/∂t = 0 in the stationary simulations: (6.15) and the corresponding probability density functions are shown in figure 6.12. They share the same shape as the ones obtained experimentally, and extreme events are once again found to be orders of magnitude more probable than in a system that follows a gaussian distribution (above ±5σ a ). It is worth noting that all the components share (almost) the exact same distribution, evidencing a clear isotropy of acceleration fluctuations.

a = Du Dt = (u • ∇)u,

Two-point statistics

The two-point statistics are calculated in a cartesian interpolated mesh generated using MATLAB in order to reduce computational costs. This is done by first generating a cartesian mesh, finding the nearest values in the OpenFOAM non-cartesian mesh and then averaging the fields into the closest cartesian grid point. A snapshot of the interpolated velocity field in an arbitrary x plane is shown in figure 6.13, which can be compared to the one shown in figure 6.10 in the original grid. In order to further corroborate that we get the same results, the mean values of the velocity and acceleration are compared to those obtained with OpenFOAM and there is only a 2% difference between both values.

The second-order structure functions are calculated with the increments along cartesian directions only. That is, we calculate

S 2 ij = [u i (j + dj) -u i (j)] 2
i, j = x, y, z, (6.16) This is done so that we can discriminate the contribution coming from each component and to further study the (an-)isotropy of the flow. Note that this does also allows us to separate longitudinal structure functions S 2 i = S 2 ii and S ⊥ 2 ij = S 2 ij , with i = j. Let us first study the different contributions for a particular Re and then compare the different Re cases.

We first consider structure functions (longitudinal and transverse) for the stream-wise components of the velocity (u y ). Figure 6.14 shows S 2y (dy) (Right) and S ⊥ 2y (dx) and S2 y (dz) (Left) for Re = 160, that is the second-order structure function for the stream-wise velocity with relative separations in the stream-wise direction and the two transversal ones. As a first observation, all the curves tend towards 2σ 2 u i , which is the asymptotic limit in (6.16) for large separation, when the velocity fluctuations are no longer correlated. Second, there is no clear two-thirds scaling in S 2y (dy), while a small inertial range appears for the transverse structure functions when the relative separations are in the x and z directions. This shows that the separations in the transversal components play a higher role than the stream-wise ones for building the inertial range. We also note that in S ⊥ 2y (dz) and S ⊥ 2y (dx) the fields become uncorrelated at dx ≈ dz ≈ 0.2D h , which is consistent with what was found experimentally, and is once again at the pore-scale. The correlation length in the streamwise direction seems larger and of the order of D h . This could again be related to the higher porosity of the bed. Figure 6.14 -Second-order structure function of u y compared to the -2/3 spectrum (red dashed line), with increments in the transversal and axial directions. A small inertial can be appreciated when the increments are in the plane perpendicular to the flow. Figure 6.15 -Second-order structure functions for the transversal components of the velocity, shown for all increments. The isotropy of the components is evident by comparing the relevant curves. The inertial range is more significant in these cases.

To further explore the difference between the dependencies at all directions, figure 6.15 shows the longitudinal and transverse structure functions for u x and u z . In these cases we do observe an inertial range, although it is over a range of scales significantly smaller than in the experimental case. This can be explained by the the high porosity, resulting in a "less turbulent" system, where fluctuations are less probable and influential. It is worth noting that S 2x (dx) is almost equal to S 2z (dz), and S ⊥ 2z (dx) to S ⊥ 2x (dz), validating once again the isotropy in the xz plane, and the equivalence between the perpendicular directions Moreover, even though it was not present in S 2y (dy), there appears to be an inertial range for S ⊥ 2x (dy) and S ⊥ 2z (dy), which points towards the fact that most of the "turbulent" contribution comes from increments and velocity components in the directions perpendicular to the flow. Figure 6.16 shows some of the previous curves for Re = 160 and 211. As with the experimental results, curves collapse at different Reynolds numbers with the right normalization, and the transversal contributions to the emergence of an inertial range are once again evidenced in the plot.

The fact that the inertial range is not as significantly large as in the experiments is a possible indication on how it highly depends on the porosity of the bed. This opens perspectives in order to perform further experiments and/or simulations of different porosities. Simulations are currently being run at higher porosity in order to get a deeper and more conclusive insight on this point. This makes sense: a lower porosity means that there are smaller regions for the flow to pass through, thus having higher local velocities and accelerations promoting locally the emergence of recirculations and shear and globally a spacially more complex and rough velocity field likely to exhibit turbulent-like statistics. The present numerical study allows a partial turbulent-like description, as only multi-scale statistics involving transverse dynamics show some analogy with turbulence. Nevertheless, since an inertial range is still detectable then we can compute δa • δu for its different components and different contributions6 Figure 6.17 shows the y and x components and relative separations (the z components is not shown as the equivalence between the x and z components has already been established). Even though δa i • δu i is not constant as observed in turbulence where an inertial range in scales exists, the difference in sign is quite impressive. Besides a clear qualitative difference can be seen between the stream-wise longitudinal velocity-acceleration structure function δu y (dy)δa y (dy) and the other contributions.

Figure 6.17 -Velocity-acceleration structure function for different components of velocity and acceleration as a function of increments in the axial and perpendicular components. It is negative at the small scales on all cases except for δa y δu y , where it is ≥ 0.

Let us recall that as it has been explained, the sign in this quantity is very important in turbulence, as it dictates the direction of the energy cascade (if it is negative the energy transfer is form the larger to the smaller structures and vice-versa). As seen in figure 6.17 when the relative separation is in the x direction, the transverse and longitudinal velocityacceleration structure-functions are entirely negative at the smaller scales (dx 0.3) before oscillating at the larger scales, which is due do the macro-scale geometry. δa x δu x (dy) is negative as well, although not as much as the two curves on the right. On the other hand, when both the correlation and the relative distances are in the stream-wise direction, the value is either zero or positive. These observations corroborate the previous hints that the turbulent-like contribution comes mostly from the dynamics perpendicular to the mean stream, which are isotropic and homogeneous. Figure 6.18 shows the total acceleration-velocity structure function (that is, when all three components of velocity and acceleration are taken into account) as a function of the dx and dy increments. In the case where the relative distances are in ŷ, there is a small range where δa • δu(dy) is negative, which is consistent with the previous results. On the other hand, when we calculate it as a function of the increments in x, the structure functions is entirely negative for the smaller scales dx < 0.2D h . This along with the results obtained with the second-order structure functions shows that the HIT "turbulent signature" comes mostly from the strucutres of the flow in the perpendicular plane, and not from the stream-wise direction. A possible explanation for those differences between stream-wise and perpendicular dynamics could be related to the hypothesis of local homogeneity generally considered in turbulence and which may better apply to the perpendicular than the stream-wise dynamics in the present case. As discussed in [START_REF] Hill | Opportunities for use of exact statistical equations[END_REF] section 4.2, inhomogeneity corrections (in particular regarding crossed velocityacceleration structure functions) come mostly from the pressure gradient term in the NS equations. In the present case, the stream-wise motion is entirely driven by the pressure gradient ∇p, and consequently the pressure plays a major role in this direction, possibly leading to deviations with respect to the simple prediction δa • δu = -2 for HIT.

Summary and final discussion

This chapter is divided in two sections: experimental and numerical results. We used the refractive index-matching particle tracking velocimetry technique to measure the local velocity field in the experiments. This was done in a system where the porosity is estimated to be ε ≈ 0.4. These experiments show that there is a "turbulent signature", where we recover results from homogeneous and isotropic turbulence, including an inertial range and the presence of a direct energy cascade at dr < 0.2D h , which appears from a multi-scale/multipore analysis, and all the different structures and contributions of the pores that, combined, lead to a "frozen turbulence" analogy, where all the turbulent-like quantities can be defined and estimated. First, based on velocity-acceleration structure functions we can extract and equivalent of the turbulent energy dissipation rate, which is the mean parameter to estimate all the "turbulent-like" quantities. We can for instance calculate the Kolmogorov lengthand time-scales and a Reynolds based on the Taylor microscale, R λ , which shows that the flow at the pore-scale is comparable to other flows that are considered to be turbulent. This is a result that reflects the multi-scale nature of porous media: even though the system is globally in a transitional (or inertial) regime that is neither laminar nor turbulent, the hydrodynamics at the pore-scale resembles that of a developed turbulent flow.

The simulations on the other hand allow us to study complementary aspects: the porosity effects and the local isotropy and homogeneity by discriminating the different contributions of a and u. We verified that the dynamics in the plane that is perpendicular to the flow is isotropic and homogeneous. It has been argued that the porosity plays a significant role as the inertial range found is not as large as the one calculated in the experiments, and we note that most of the contribution for the inertial range comes from the transversal components. The inertial range is more significant when it is calculated with the transversal components of the velocity, evidencing that most of the "turbulent signature" comes from the xz plane. Most of the negative contribution, which accounts for a direct energy cascade, comes once again from the perpendicular directions, whereas it is positive when it is calculated parallel to the flow. This last result can be related to the inhomogeneity caused by the pressure that drives the whole stream-flow.

These are all interesting results for the volume averaging theory, which is discussed in the following chapter, where we can safely assume that the system is isotropic and homogeneous in the direction perpendicular to the flow. It also shows the limits of the V.A.T, given that we could be filtering out these pore-scale effects when averaging, although this will not affect the macroscopic results given that, for example, the Darcy and Darcy-Forchheimer laws can be derived using this formalism [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF].

Chapter 7

Some results on the closure problem in the averaged equations

In the following chapter we are going to explore some aspects of the closure of the volumeaveraged equations in porous media by studying the non-closed terms and how they are affected by the finite wall effects.

Averaged equations

As it has been investigated in the two previous chapters, the hydrodynamics at small scale is very different than at global scale in fixed beds. In particular, we have seen that, though stationary, the velocity fluctuations present a turbulent-like behavior at the smaller scales, even though the flow is laminar. We wonder then whether the impact of the small scale turbulent-like behavior might be present in the large-scale dynamics. Can a sub-grid homogenized description of this complex motion be derived to model the large scale dynamics? Indeed, when modeling porous media at an engineering scale (i.e., global) we need a framework to average the equations of motion so as to avoid solving all the smallest scale of the process since the time required to run the simulations might be impractical, or because the physical processes occurring at small scales are not of technical relevance when constructing a reactor for instance.

We thus need a mathematical tool that can include both the micro-and macro-scale aspects of the equations of motion in porous media, where the system is intrinsically multiscale. Volume-averaged techniques are usually used in these kinds of problems, where we relate both the fluid and solid phase contributions. This averaging technique generates more unknown than equations and requires closure models in order to have a full description of the problem. We provide an overview of the method proposed in [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF] and [START_REF] Jackson | The Dynamics of Fluidized Particles[END_REF], where the reader will find further details if needed.

The averaging is developed in terms of locally spaced averages that are defined in terms of a weighting function (or a filter function) g in a way that the overall average value at position x for any point function q is defined as1 :

q f (x) = Vf q(y)g(|x -y|)dV y ,
where the integration is done over the whole fluid volume V f . It is worth noting that in other works, notably in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF], a representative elementary volume (R.E.V) is chosen and the averaging is done over that control volume. This is a particular case of the formalism proposed by Jackson, given that g can be for example a heavyside step function that defines the R.E.V. This filter approach in contrast is similar to that used in subgrid L.E.S models (M. [START_REF] Kuerten | Point-Particle DNS and LES of Particle-Laden Turbulent flowa state-of-the-art review[END_REF]. In order to make the averaging process independent of the specific filtering function g, the following scale separation needs to be verified:

d r L,
where d is the particle diameter, L the macroscopic length scale and r the filtering length scale. If we average the mass conservation and stationary Navier-Stokes equations for the (incompressible) fluid phase with this operator we obtain the following expressions [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF])

∂ε ∂t + ∇ • (ε u f ) = 0, and (7.1) ρε∇ • ( uu f ) = ∇ • (ε σ f ) - p Ŝp σ(y) • n g(|x -y|)dS y , (7.2)
where ρ is the fluid density, ε its porosity and σ the fluid stress tensor2 . The last term on the right hand side in equation ( 7.2) is integrated over a particle surface S p , and it is summed over all particles indexed by p.

Up until here the formalism is equivalent to the one proposed in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] for homogeneous R.E.Vs, and consecutive developments for specific cases such as in turbulent flows in porous media (Soulaine & Quintard, 2014) and in moving porous media [START_REF] Wang | Volume-averaged macroscopic equation for fluid flow in moving porous media[END_REF]. The issue with this expression is that this term cannot be directly interpreted like a "lagrangian" force [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF]. That is, we want the force to be described from the spheres' "point of view". In order to do this, it is necessary to describe the force at the position of the particles' mass center, x p . If we do the following expansion 7.3) This series can be truncated with acceptably small error [START_REF] Jackson | Locally averaged equations of motion for a mixture of identical spherical particles and a newtonian fluid[END_REF]. The first term in the right hand side of equation ( 7.3) can now be interpreted as a force term that represents the force exerted on the particles by the fluid, as felt by the spheres at their center of mass. We express the first term in the RHS as n F p , where the average is done over the particles and n = N/V is the density of particles per volume (N is the total number of spheres in the volume V). On the other hand, the second term in equation ( 7.3) involves the first moment of the loads on the particle, namely the torque (anti-symmetric component) and the stresslet (symmetric component) [START_REF] Guazzelli | One sphere in Stokes flow[END_REF]. We will express this term as M p , where M represents the first-order moments of the loads and the average is once again done over all the solid particles. It is worth noting that this term does not explicitly appear in the development proposed in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF], where there is only a force term. This is because so far we have not done any stronger assumptions than the separation of scales, whereas in [START_REF] Whitaker | The Forchheimer equation: A theoretical development[END_REF] we also need local homogeneity. It is not evident whether the forces alone suffice to close the equations in the case where there are wall effects present, and where there might be inhomogeneities in radial direction. We decompose the velocity product as

g(|x -y|) = g(|x -x p |) -(y -x p ) • ∇g(|x -x p |) + ... we have that p Ŝp σ(y) • n g(|x -y|)dS y = p g(|x -x p |) Ŝp σ(y) • ndS y -∇ • g(|x -x p |) Ŝp [σ(y) • (y -x p )]ndS y + ... (
u i u k f = u i f u k f + u i u k f
, where u i = u iu i f are the velocity fluctuations.

Finally, using equation (7.3), equation ( 7.2) can be then rewritten as: 7.4) where Σ is an effective stress tensor defined as

ρε∇ • ( u f u f ) = ∇ • Σ -n F p , (
Σ = ε σ f -u u f + M p . (7.5)
We can observe from equations (7.4) and ( 7.5) the solid-fluid interactions due to the force and the first moment, that is also averaged over the particle phase. It is also worth noting that we have a Reynolds tensor-like term that involves the velocity fluctuations one-point correlations. These equations can get closed by studying ∇ • Σ and F p as a function of the different parameters involved in the system.

Because of the geometry that we are working in (axisymmetric), we are particularly interested in studying the z-component of equation (7.4). In cylindrical coordinates it reads:

ρ u r f ∂ u z f ∂r = ∂Σ zz ∂z + 1 r ∂(rΣ zr ) ∂r -n F z p (7.6)
If we neglect the variation of the averaged quantity with the radial direction one gets the same balance as presented in Chapter 4. Once again, there are three control parameters: the porosity ε, the Reynolds number Re and the separation of scales, quantified by D/d. We will explore the effects on the different terms of equation ( 7.4) by using several ensembles. This is done so as to have better statistics on the averaging in the particle phase, as more than 1000 particles is necessary. Simulations are still ongoing at this point so as to achieve better statistical convergence3

This set of simulations will allow us to study the porosity effects by varying ε between 0.9 and 0.5; the Reynolds effects by setting Re = 100 and 200, and finally the border effects by comparing D/d = 10 and 5 while keeping all the other parameters fixed4 . We will also be reanalyze some other aspects of the simulations presented in chapter 5.

Bulk averages

In this section, we investigate the average of several quantities over the whole control volume. One of the main results from chapter 4 is that

dp dz = 1 V p F spheres z + F wall z , (7.7)
where the sum is done over all the solid spheres and F wall is the force at the lateral wall of the bed. This result can be proved using equation 7.6, neglecting the variation of the averaged quantity in the radial direction, and it is valid when averaging over the computational domain.

Additionally, we can calculate the force probability density function in the bulk, and even though the results might not have fully converged, they still give us some clues of the overall force distribution of the system. Figure 7.1(Left) shows these results for the force calculated for 1000 spheres normalized by their mean value for D/d = 10 and ε = 0.5 and two different Reynolds numbers Re = 200 and 100. There is a peak at F z / F z ≈ 1.35, which means that there is a good amount of spheres that have a solid-fluid force 20% higher than the average. This is different from the force distribution presented in [START_REF] Hardy | Simulation of the flow past random arrays of spherical particles: Microstructure-based tensor quantities as a tool to predict fluid-particle forces[END_REF] for tri-periodic beds. This difference of the force distribution when compared to its mean force can consequently be attributed to finite wall effects.

To further explore the wall effects, figure 7.1(Right) shows the same results for Re = 200 and D/d = 10 and 5 for comparison. The average of the latter case was done for 360 particles. We can see that the peak remains almost at the same place, so there are still a significant amount of spheres that feel a force higher than the average, and there is a higher probability of finding spheres with a higher F z / F z when compared to D/d = 10.

Moreover, in order to study porosity effects using an extreme case, we calculated the force distribution for 9 different fixed beds (with a total of 900 spherical particles) with D/d = 10 and porosity ε = 0.9. The pdf is shown in figure 7.2 for Re = 200 and the ε = 0.5 case is shown for comparison. The difference between the pdfs is remarkable, as apart from a small peak, there is no main peak for ε = 0.9 when compared to the less porous case. This can be explained by a local almost-homogeneity, as will be explored in the next section, so that the force is (almost) evenly distributed among the spheres. It is worth noting however that the peak is at F z / F z ≈ 2.8. As a first step, we calculated the different components of the Reynolds stress tensor. All the crossed components are near zero

u x u y f ≈ u x u z f ≈ u y u z f ≈ 0,
when compared to the other components, so we are interested in the diagonal components of the stress tensor. Moreover, we verified that

u x u x f ≈ u y u y f ≈ u r u r f ≈ u θ u θ f ,
which is consistent with what was found in the previous chapter at the pore-scale, where we observed that the transversal components are comparable. Figure 7.3 shows u r u r f and u z u z f normalized by the square of the averaged streamwise velocity in the fluid phase, ( u z f ) 2 for two confined cases with scale separations D/d = 5 and 10.

As a first observation, u r u r f increases with the Reynolds number whereas the opposite is true for u z u z f , and the latter is dominant over the contribution of the transversal component, which points to the fact that the Reynolds tensor is axisymmetric [START_REF] Pope | Turbulent Flows[END_REF]. Moreover, there is an anisotropy between the r and z components, and these results agree with what has been observed in [START_REF] Mehrabadi | Pseudo-turbulent gas-phase velocity fluctuations in homogeneous gas-solid flow: fixed particle assemblies and freely evolving suspensions[END_REF] for tri-periodic arrangements and in [START_REF] Uhlmann | Interface-resolved direct numerical simulation of vertical particulate channel flow in the turbulent regime[END_REF] in particulate channel flows. This anisotropic behaviour is also consistent with what was found in previous chapters at micro-scale: as was concluded using the numerical results, the main contribution to the inertial range observed in the smaller scales in the second-order structure function (the 2/3 scaling) seems to come from the transversal components of the velocity for high Reynolds numbers. Indeed, what we observe here is that as the Reynolds number increases the radial component has a tendency to become more significant, whereas the opposite is true for the axial one, leading towards a more isotropic state where both contributions become more and more comparable. We also quantified this by computing the ratio between the longitudinal and transversal components of the integral lengths, which tends towards unity with increasing the Reynolds number.

We observe that the Reynolds stress tensor is more significant (when compared to the averaged velocity) when the D/d increases. In This might be due to the fact that for higher D/d the average porosity is lower (so there is less fluid volume) (De Klerk, 2003). Because the pores are smaller, more fluctuations can be generated and consequently have a higher contribution. If we think qualitatively of the different passages of the fluid flowing through the medium, we can expect to have less velocity fluctuations where the paths are less tortuous and there is a high local porosity (such a region could be for intance near the reactor wall), whereas at places of low porosity we might have more extreme events, such as the ones presented in chapter 6, where we observe that high values of acceleration are probable (also reported in [START_REF] Holzner | Intermittent lagrangian velocities and accelerations in three-dimensional porous medium flow[END_REF]). All of this could be interpreted as a porosity effect, but in this case, the porosity itself is due to the confinement6 . We will now consider a fixed bed with D/d = 10 and porosity ε = 0.9. Figure 7.4 shows the longitudinal and transversal components of the Reynolds tensor for both beds with D/d = 10 and porosity ε = 0.5 and 0.9 for comparison. As it can be observed, they decrease with increasing fluid volume fraction, as expected from the previous discussion. In this case, the change in porosity is not caused by the walls but by the imposed loose packing. The variation with Re of the Reynolds tensor is less evident when ε = 0.9, and the results do not seem to tend to isotropy as in the case with ε = 0.5, or at least not as "fast" as in that case. In fact, the slope in the curve for ε = 0.9 is slightly smaller than the ε = 0.5 case (also observed in [START_REF] Mehrabadi | Pseudo-turbulent gas-phase velocity fluctuations in homogeneous gas-solid flow: fixed particle assemblies and freely evolving suspensions[END_REF]). This means that the symmetry is already well defined: there are less variations in the radial direction with the Reynolds number and the same happens for the stream-wise direction. If we say that the flow is almost Poiseuille-like, the axial symmetry is well defined. The difference in mean value (lower for ε = 0.9) is consistent with what is observed in [START_REF] Mehrabadi | Pseudo-turbulent gas-phase velocity fluctuations in homogeneous gas-solid flow: fixed particle assemblies and freely evolving suspensions[END_REF] and [START_REF] Reddy | Recent developments in experimental (piv) and numerical (dns) investigation of solidliquid fluidized beds[END_REF]: the increase of the solid volume fraction leads to a decrease in anisotropy, as lower porosities do not allow the formation of distinct wake structures.

Radial inhomogeneities

Let us now go into a smaller set of scales, say meso-scales, which are in between the pore scale and the global scale. That is, we are going to take into account the local inhomogeneities and study the different closure terms filtered at different scales.

Given the fact that the ∂ r operator is present in equation ( 7.4), we are interested in the radial inhomogeneities. We explore this by averaging different quantities over cylindrical shells ∆r i = [r i , r i + δr i ], where r i is the considered position of the filtered quantity when plotting the results.

The average is done over the other two components z and θ as well. For instance, the radial profile for a variable k is

k f r (r) = 1 V i r i +δr i ri 2π 0 H 0 k(r, θ, z) dr dz dθ, (7.8)
where V i is the fluid volume made up by the two concentrinc rings delimiting the integrational domain. Even though we are doing a radial average, for the sake of clarity we will leave out the • f r operator when talking about locally averaged quantities.

First glance: velocity and porosity profiles

Before going onto studying the different terms of equation ( 7.4), we shall first calculate the radial profile of the velocity u z (r). Figure 7.5a) shows the velocity profile normalized by the superficial velocity U for a particular case of D/d = 10 and Re = 20, 100 and 200. All the curves collapse, indicating that the velocity profile does not depend on its mean superficial velocity. This is not the case with the Reynolds stress, as evidenced in figure 7.3 where they are normalized by the average velocity in the fluid phase, which marks one difference between the behavior of the mean flow and its fluctuating component.

In contrast, if the porosity of the medium is changed and there is a higher fluid volume fraction both the velocity and porosity profiles change significantly, as shown in figures 7.6 a) and b), where the curves for ε = 0.5 are shown for comparison in a lighter shade. The porosity profiles are normalized by the mean porosity ε for visualization purposes, and no clear oscillating pattern is observed, which is due to the small number of solid spheres which are randomly inserted inside the cylindrical domain. The velocity profile in constrast, shows oscillations with lower amplitudes, with a higher value near the center instead of near the walls (unlike in the ε = 0.5 case). This may be reminiscent of the beginnings of a Poiseuille profile, shown in chapter 3, which reaches a maximum value near r = 0. We will calculate u z u z (r) by computing its mean over concentric rings between two radii, ∆r i = [r i , r i + d/2], with varying r i . This choice of averaging is done so that we cover half a period of the mean velocity profile. As expected from figure 7.3, its mean value decreases with the Reynolds number. The border effects are evidenced by the increase in value near the wall, which is 33% larger than the bulk-averaged value u z u z f . This might be due to the fact that near the walls not only the fluid volume fraction is larger, but there are also fluid-spheres and fluid-walls interactions which might enhance velocity fluctuations (instead of just having spheres-fluid interactions at the center). It does not show the same defined oscillatory pattern found in the velocity radial profile, but there are still inhomogeneities present in the radial direction, so that the ∂ r u z u z term in equation (7.6) might be different than zero at these scales. Figure 7.8a) shows the same quantity but averaged for concentric cylinders in d/2 intervals, and the overall behavior is captured by averaging using this resolution. Figure 7.10 shows the averaged z-component of the forces in the radial direction F z (r), once again calculated using d/2 bins for all of the cases studied, normalized by the bulk average F z . Let us first analyze the case with D/d = 10 and two different Reynolds numbers, Re = 100 and 200. As reflected by the force profile, there is an almost oscillatory behaviour and the force seems to be evenly distributed on the spheres. Moreover, the mean radial force decreases as the spheres approach the walls. This might be explained qualitatively using the following arguments: the force is affected by the velocity gradient, and even though the porosity near the walls is higher and the velocity as well, the flow could be smoother there and the velocity gradient is weaker. On the other hand, in regions of smaller porosity velocity (near the center of the bed, at r = 0) gradients are higher and have a higher effect on the particle-fluid force.

We also observe that the radial profile of the forces does not change with the Reynolds number under this normalization, which follows the same trend observed in the velocity profile, where all curves collapse when normalized by the superficial velocity.

Figure 7.10 shows the same profile for a fixed Reynolds number (Re = 200) and porosity ε = 0.5 and two different particle-to-tube ratios D/d = 10 and 5, this time normalized by the tube radius R for visualization and comparison purposes. In the D/d = 5 case the forces seem to fluctuate less, and they decrease less when the spheres are approaching the walls. This might be due to the fact that there is not enough place for the forces to "relax" towards a lower value: there is a large effect from the sphere diameter d so that oscillations with D/d = 5 have a bigger period.

Finally, we study a third case at fixed Reynolds number Re = 200 and D/d = 10 and different porosities ε = 0.5 and 0.9, shown in figure 7.11. In this case the ε = 0.9 profile is higher than the one for ε = 0.5, which is in agreement with the peak found in the bulk pdfs (see figure 7.2). This seems to show that in fact the force might not be affected so much by the velocity fluctuations (which are higher for ε = 0.5), but by the local mean velocity: in the case of ε = 0.9 the local dynamics are not driven by the velocity fluctuations with the same intensity as for ε = 0.5.

Summary

The wall effects are all evident in the packed beds with ε = 0.5 and are reflected in the velocity and porosity profiles and in a lesser extent in the Reynolds tensor. This might lead to the fact that it is necessary to account for the wall-fluid interactions when there is a high solid concentration, which add inhomogeneities in the radial direction at the meso-scale. We also noted that the force is evenly distributed among the spheres, even when they are confined, and this can be qualitatively explained by considering the behavior of the local velocity field.

In contrast, even though they are confined, these effects are not seen when the fluid volume fraction is higher, in particular for ε = 0.9, which might point towards the fact that there is a critial porosity at where the border effects are no longer evident, although further studies need to be done to reach this conclusion, for example by studying a loose bed with D/d < 5.

Macroscale

Chapter 5 explored the global border effects in packed beds where the sphere-to-reactor diameter ratio is lower than ten (D/d < 10). We did this both experimentally and numerically, and both results complemented each other, as they allow us to study different aspects of the same problem.

The first aspect we studied is the variability of the pressure field and pressure gradient, and how it is affected by D/d and the Reynolds number. This variability is associated to the random arrangement of the spheres and in order to quantify it we measured experimentally the pressure drop for nine different packed beds with D/d = 10 under the same conditions. The standard deviation of the results (compared to the ensemble average) was found to be approximately 4%.

We further related this variability of the different configurations with the one found within a single packed bed, but by measuring at different angles of the reactor wall. This is of interest because typically pressure sensors are placed along the walls at a certain angle, but there is a source of error that is linked to the local arrangement of the spheres, and how this measurement would vary if we measured at a different angle. We explored this numerically and we found that the standard deviation for D/d = 10 is close to 4%, which point towards the conclusion that there could be an equivalence between doing several different realisations of the experiment and placing the sensors at different places of a particular packed bed.

We also observed that this variability increases with decreasing D/d, which is due to wall effects on the pressure drop. If the observed trend is correct, this variability could go up to 10% for D/d < 5 and could have an effect on the several models proposed in the literature.

As it has been said, there is no general agreement upon a correlation that describes the pressure drop in terms of ε, Re and D/d. However, if there are two models that are within the, say, 5% errorbar, they should be equivalent, given that the variability should be taken into account.

The second global aspect that we studied was how the pressure gradient is affected with D/d. We found that it follows an Ergun-like relation for all Reynolds numbers for both the experiments and simulations, and that it decreases with decreasing D/d. This can be explained by the fact that when D/d decreases the porosity increases: this means that there is less hydrodynamic resistance for the flow to pass through the bed so that the pressure drop is weaker when the porosity is higher.

Using conservation arguments we found that the pressure drop is proportional to the sum of the spheres and wall forces when we consider the whole computational domain. This allowed us to quantify the contribution from the spheres and the walls. We observed that for D/d = 10, 98% of the pressure gradient comes from the contribution of the spheres, whereas it is 95% for D/d = 5. Because of this we conlude that the change in the pressure drop when compared to the Ergun relation for non-confined beds is due to how the particles are arranged and how their average porosity changes with D/d, rather than by the presence of the walls themselves.

Microscale

The first part of chapter 6 was dedicated to the experimental study of the local hydrodynamics of the flow with Reynolds number ranging between 120 and 210 passing through a fixed bed with D/d ≈ 6. We observed qualitatively that the Lagrangian dynamics of the tracer particles is not trivial at all: although stationary, the flow shows a clear spatial complexity, with shear zones and rotation-dominated structures, with recirculation and stagnation points.

This complex behavior motivated us to compute one-and two-point statistics of the fluid velocity and acceleration fields. We found that high values acceleration are probable, as already observed in other studies, and that the acceleration probability density function presents exponential tails similar to the ones observed in fully turbulent flows. We then computed the two-point eulerian correlation function and we found that two elemens of fluid are spatially correlated at lengths smaller than approximately d/5.

We also computed the streamwise and transversal correlation (pore-scale) lengths as a function of the Reynolds number, and their ratio tends towards unity as Re increases. This seems to point towards local isotropy: as the Reynolds number increases both correlation lengths become more and more comparable, so that there is less of a difference between the two directions of the flow.

The spatial complexity of the flow, the acceleration pdfs and the small-scale correlations inspired us to use the tools typically used to study and characterize turbulent flows. We first computed the second-order structure function, where an inertial range is observed in homogeneous isotropic turbulence (HIT) by its two-thirds scaling. Surprisingly, we observe this turbulent(-like) scaling at scales smaller than the correlation length.

If there are scales that are similar to the inertial range observed in HIT, there has to exist an energy cascade at these scales. This is quantified by the correlation of velocityacceleration spatial increments, which is directly linked to the energy dissipation rate and is constant at the inertial range. Effectively, we observed that δu • δa remains relatively constant at this range of scales, and that it has a negative value. Following turbulent-like arguments, we thus say that there is a direct energy cascade at these scales, where the energy goes from the larger scales to the smaller ones.

These results show that at the smallest scales we have a stationary turbulent-like system, even though the flow is laminar. We then computed and estimate of the energy dissipation vs. σ 3 u , and found that they follow are reasonably linear relation, which is what is observed in fully developed HIT.

In order to have a full turbulent-like description we used to calculate the Kolmogorov length and time scales, which we used to calculate the Reynolds number associated to the taylor microscale in turbulence, R λ . We found that it is highly comparable to the ones found in wind tunnel experiments where fully turbulent flows are studied. This validates once again the argument that there is a "real" turbulent-like behavior at the smallest scales of the medium.

In the second part of the chapter we presented the results for simulations computed at two of the Reynolds numbers studied experimentally, similar D/d and with a relatively lower porosity. This allowed us to study the influence of the porosity on the local hydrodynamics. As a first result, the velocity pdfs differ slightly from the experimental results, and we attribute the difference to the resolution of the simulations, as for instance we cannot experimentally resolve the scale near the beads. We also observed that the plane perpendicular to the flow is isotropic, as the pdfs of the two transversal components are almost equals. The acceleration pdfs share the same exponential tails as the ones observed experimentally.

We took advantage of the additional data provided by the simulations and we decomposed the experimental analysis by computing magnitudes relative to the transversal and streamwise components. We computed the transversal and longitudinal (cartesian) secondorder structure functions, and we found a smaller "inertial range" (that is, the range of distances that it spans is significantly smaller). This might be qualitatively explained using the following argument: as the porosity is larger and the pore-scale larger, there are less velocity fluctuations and/or extreme events (such as extreme acceleration values) that can be responsible for a turbulent-like behavior, which is driven by velocity fluctuations.

Moreover, there was almost no inertial range present in the streamwise longitudinal structure function, S yy (dy) (y is the direction of the flow). On the other hand, when either one of the transversal components are involved, there is a small two-thirds scaling. This seems to indicate that in fact the turbulent "signature" comes from the planes that are perpendicular to the direction of the flow (that is, the rθ or xz plane).

We then computed the velocity-acceleration structure function for the different components. We did not observe a constant value, but this can be explained using two different arguments. The first one is that we are not calculating exactly the same quantity that we calculate in the experiments, as in the experiments the increments (or relative distances) are calculated in spherical coordinates, wheareas here we compute it in only one cartesian direction. The other argument is that it might also be a porosity effect. As it has been specutaled with the second-order structure function, the turbulent-like behavior is less significant when the porosity is higher.

Nevertheless, we were able to extract information from these quantities. As we observed in experiments, δu • δa is a negative quantity at small scales. In fact, what we observed in the simulations is that once again, the negative contribution comes from the trasnversal components, which validates our point that most of the turbulent signature comes from the structures of the flow in transversal planes. The streamwise velocity-acceleration longitudinal structure function remains positive or is equal to zero at all scales.

The derivation of the δu • δa -relation was done for homogeneous flows. In fact, what our data seems to show is that the transversal components are homogenous, whereas the streamwise direction is not, which is consistent with the fact that most of the inhomogeneity corrections come from the pressure gradient term in the NS equations, which is what drives the entire system in the streamwise direction.

Mesoscale

Finally, in chapter 7 we presented some preliminary/exploratory results on the closure problem in the averaged equations. We studied the different terms involved in the non-closed equations by quantifying how they vary with D/d, Re and ε.

We found that when bulk averaged the only non-zero components of the Reynolds stress tensor are the ones in its diagonal, and we paid attention to the streamwise and radial components. We observed that as the Reynolds number increases velocity fluctuations tend to become isotropic, reflected by the fact that u z u z decreases with Re while u r u r increases. This is in agreement with what was found at the pore-scale in chapter 6 with the longitudinal and trasnversal correlation lengths. We found these results for beds with ε = 0.5 and D/d = 5 and 10. Moreover, confinement effects are observed in the mean value of the components of the Reynolds stresses, as its streamwise component can be up to 4 times larger in a bed with D/d = 10 when compared to one with D/d = 5. This might be explained by the fact that due to the presence of the wall, the porosity increases so that there are less velocity fluctuations. It is thus a porosity effect (as was also speculated and explored in chapter 6), but the porosity itself changes because of the presence of the walls (as was observed with the pressure drop in chapter 5).

In order to study porosity effects for a case where the other parameters remain unchanged, we studied how the Reynolds stresses vary for D/d = 10, Re = 200 and porosities equal to 0.5 and 0.9. We recovered the previous results, where the averaged Reynolds stresses are lower for the loose bed when compared to the packed one.

As for the forces, we observed that the force balance found in the global effects chapter (Ch. 5) could be interpreted as a particular case of the bulk-averaged equations. It is interesting to note that the walls account for 2% to 5% of the balance, while the fluid-wall forces do not appear in the volume-averaged equations.

We then studied the local radial inhomogeneities. We observed the wall effects in the porosity and velocity radial profiles: ε(r) tends to unity near the walls, which offers a less resistant path for the flow to pass through and consequently there are higher velocities near the walls as well. This was observed for ε = 0.5. On the other hand, the profiles remain almost constant for high porosities (ε = 0.9), which means that the wall effects are significantly stronger at lower porosities. The inhomogeneity generated by the wall presence is also observed in the Reynolds tensor, which increases near the wall for ε = 0.5 and D/d = 10 and 5. In contrast, it does not significantly change when ε = 0.9. These results show that radial inhomogeneities should be taken into account when doing averaging procedures at the meso-scale, in particular for beds with low porosity.

We finally studied the forces radial distribution. For D/d = 10 and ε = 0.5 we observe a slight decrease in the sphere-fluid forces as the spheres approach the wall, but there are no border effects as siginficant as the ones observed in the velocity profiles or in the Reynolds tensor, as they are almost evenly distributed. This might be explained qualitatively using the following arguments: the force is affected by the velocity gradient, and even though the porosity near the walls is higher and the velocity as well, the flow could be smoother there and the velocity gradient is weaker. On the other hand, in regions of smaller porosity velocity (near the center of the bed) gradients are higher and have a higher effect on the particle-fluid force.

In contrast, when we consider a loose bed ε = 0.9 the forces are once again evenly distributed and the curve is not as different as the one for ε = 0.5. In particular, the ε = 0.9 force profile is higher (in average) than the one for ε = 0.5, which is in agreement with the peak found in the bulk pdfs. This seems to show that in fact the force might not be affected so much by the velocity fluctuations, but by the local mean velocity.

All of these results give clues as to how the wall effects are present and reflected in the radial inhomogeneities at the mesoscale. In particular, the wall forces have not been accounted for in the volume-averaged model and they should be, given the fact that they change the mesoscale hydrodynamics. In this chapter we were able to validate, corroborate and use results that were obtained both in the global and micro scale flows, which reflects the multi-scale nature of the problem.

Chapter 9

Perspectives

The results obtained during this thesis open many perspectives for future studies. We will candidly point out a few of the questions that arised and remained open during this threeyear work, and aspects that we would like to continue studying and exploring.

At the global scales

Measuring the pressure drop on fixed beds can be seen as a relatively simple and straightforward experiment and it has been widely studied for the last seventy years or so, mainly due to its numerous applications, but despite its apparent simplicity a full model that describes the pressure drop when border effects are significant has not yet been found. Nevertheless, in light of the results obtained regarding the variability of the pressure drop (and how it increases with decreasing D/d) we consider that the variability should be taken into account when proposing new models, and when measuring the pressure drop.

New experiments could be done with D/d < 5 to verify the trend that we found for the variability. This could also be repeated for beds made up of different particle shapes, and used as a proxy or guide for the error estimation.

At the pore scale

There are several aspects that could be explored to further study the turbulent-like behavior at the pore-scale:

• A conditioned statistical study could be done on the two-point statistics by conditioning the computation on pores with similar flow patterns. Would be obtain the same scalings if we only took into account pores where the flow is channel-like? What about the pores where we only observe circulation points? This would help to quantify the interplay between the different kinds of local behaviors.

• In K41 theory intermittency appears when S p = C p (ε dr) p/3 , which is at p = 4. This itermittency could be studied at the pore-scale by computing more statistics on the tracers trajectories, either by doing more experiments or filming more time (seeding more tracers into the setup would present further experimental complications, especially with tracking algorithms).

• It would be interesting to observe at which Reynolds number this behavior starts. Is it valid for all flow regimes?

• Another interesting experiment using refractive-index matching PTV would be the injection of a colorant into the flow and to observe the overall dissipation and transport, especially when the confinement is strong. We would expect to be able to observe the preferential sampling of the flow, and how it is affected by the wall presence. This could be then compared to simulations by the injection of a passive scalar.

The closure problem

The need for the consideration of the wall effects has been evidenced at the meso-scale.

In fact, it would be interesting to develop a formalism that takes into account the walls themselves (a priori, this term would appear in the fluid-solid forces term).

More computations would be needed to have more statistics, and we would like to explore a wider parameter range. We could further explore the tendency towards isotropy and its relation to porosity, the particle-to-tube diameter ratio and Reynolds number.

Figure 1 . 1 -

 11 Figure 1.1 -Wind and solar power are intermittent in nature, as it is usually the case with most of the renewable energy sources. Figure taken from Lawson (2019).

Figure 1 . 2 -

 12 Figure 1.2 -A schematic of an AA-CAES system. Figure taken from Mozayeni et al. (2017)

Figure 2 . 1 -

 21 Figure 2.1 -Figure taken from Wood et al. (2020) to show the multi-scale hierarchy present in porous

Figure 2 . 2 -

 22 Figure 2.2 -Simulated (left) and real-life (right) fixed beds of randomly-arranged spherical particles.

Figure 2 . 3 -

 23 Figure 2.3 -Predictions of mean porosity values for packed beds with D/d < 30.Figure taken from Guo et al. (2019). The works referred to in the legend correspond respectively to Sato et al. (1973); Dixon (1988); Fand & Thinakaran (1990); Foumeny et al. (1996); Zou & Yu (1995); De Klerk (2003); Benyahia & O'Neill (2005); Ribeiro et al. (2010); Cheng (2011)

Figure 2 .

 2 Figure 2.4 -Radial porosity variations for randomly packed beds with different D/d ratios. Figure taken from Mueller (1992). The works referred to in the legend correspond respectively to Roblee et al. (1958); Benenati & Brosilow (1962a); Ridgway & Tarbuck (1968); Goodling et al. (1983); Mueller (1992); Giese et al. (1998)

Figure 2 .

 2 Figure 2.5 -General outline of the results presented in chapters 4 through 7, along with the relevant parameter ranges.

Figure 3 . 1 -

 31 Figure 3.1 -Overview of OpenFOAM structure, taken from the OpenFOAM user manual(Greenshields,

Figure 3 . 2 -

 32 Figure 3.2 -A typical control volume used for the Finite Volume Method. Figure taken from Jasak(1996)

Figure 3 Figure 3 . 3 -

 333 Figure 3.3 -Illustration of the DEM method for a generic pair of spheres.

Figure 3 . 4 -Figure 3 . 5 -

 3435 Figure 3.4 -Process of constructing a random periodic arrangement of solid spheres

Figure 3 . 6 -

 36 Figure 3.6 -Workflow of the mesh generation process.

Figure 3 . 7 -

 37 Figure 3.7 -Snapshot of a sphere in the border of a periodic domain in the particular case of a Face-Centered Cubic unit cell. The domain was repeated three more times, so that the periodic geometry can be appreciated.

Figure 3 Figure 3 .

 33 Figure 3.9 -A generic background mesh generated with the blockMesh utility, composed by an O-H butterfly topology. Left: the xy plane. The center square H block can be appreciated, along with the variation in resolution. Center: Radial coordinates. Right: Relevant lengths: the cylinder radius R and its height H.

Figure 3 .

 3 Figure 3.11 -Force over the sphere for different grid resolutions, measured in terms of cells per particle diameter (cpd), compared with the force measured by Hill et al. (2001).

FFigure 3 .

 3 Figure 3.12 -Non-dimensional force on a simple cubic array of ε = 0.592 as a function of Re obtained by OpenFOAM and compared to the one calculated by Hill et al. (2001). Inset: relative error between both results.

Figure 3 .

 3 Figure 3.13 -a) Different cases of a simple-cubic ordered array. b) Relative error for the non-dimensional force on a simple cubic array of ε = 0.592 as a function of Re obtained by OpenFOAM, compared to the one calculated by Hill et al. (2001) for the SC cases shown in a).

Figure 3 .Figure 3 .

 33 Figure 3.14 -Left: Pipe of radius R and height H = 40 × R. Right: Velocity field for a case with Re = 10.

Figure 3 .

 3 Figure 3.16 -Comparison between the pressure drop calculated using the SIMPLE and PIMPLE algorithms for a fixed bed of D/d = 8 and Re = 200

Figure 3 .

 3 Figure 3.17 -Comparison between the velocity probability density functions of the fluid phase calculated using PIMPLE and SIMPLE algorithms. The colored curves represent different times and the black solid curve represents the SIMPLE solution.

Figure 3 .

 3 Figure 3.18 -Top: Relative error between the force over a single sphere calculated by SIMPLE and PIMPLE for thirteen random spheres of the bed, each in a different color. Bottom: Time-averaged relative error for 13 random spheres of the bed.

Figure 4 . 1 -

 41 Figure 4.1 -Schematic representation of the experimental setup.

Figure 4 . 2 -

 42 Figure 4.2 -Different sections of the setup. They allow flexibility at the assembling of the experiment

Figure 4 . 3 -

 43 Figure 4.3 -Test section of the reactor with a partially filled bed. Pressure sensors are attached to the cables, and are a ∆h distance apart from one another.

Figure 4 . 4 -

 44 Figure 4.4 -Stainless steel particles in different sizes were used.

Figure 4 . 7 -

 47 Figure 4.6 -(Left): Raw signal from the flowmeter. It has overshoots in the extremes of the square signals.(Right): Flow meter signal after being filtered by a low-pass filter and a signal attenuator.

  Figure 4.8 -Calibration curve for one of the pressure sensors.

Figure 4

 4 Figure 4.9 -A typical pressure signal in Volts captured by the sensors.

2 ) 2 )

 22 Figure 4.10 -Cumulants of the mean value µ and the square of the standard deviation σ respectively.

Figure 4

 4 Figure 4.12 -Schematic representation of a section of the experimental set up used to do Particle Tracking Velocimetry.

  Figure 4.13 -Schematic representation of the optical setup viewed from two different angles: a) from the front and b) from above.

Figure 4 .

 4 Figure 4.14 -Top: empty container filled with air. Middle: the container is filled with hydrogel beads.

Figure 4 .

 4 Figure 4.17 -Post-processed images taken by the two cameras of a generic experiment, the frames are all superposed.

  Figure 4.17a) shows

Figure 4 .

 4 Figure 4.18 -Scheme of the Center to Rays process. Red circles show the particles that are seen in the pictures. This step allows us to draw lines between the cameras and the observed particles on each frame.Figure adapted from Dumont (2021)

  Figure adapted from Dumont (2021)

  Figure 4.19 -Hydrogel beads were bathed in colorant to as to have an estimate of the size of the pores.

  Figure 4.20 -Coordinate system used in these experiments. The xy plane is the one recorded by the cameras and the flow is in the y direction. The depth is in the z-coordinate, and it is perpendicular to the flow, along with x.

  Figure 4.21 -a) Comparison between a real point (•) and the point found by the PTV algorithm ( ). b)

  Figure 4.22 -a) Generic trajectories fitted by an order 3 polynomial. The dots represent the real data and the solid lines are the fits. b) Probability density functions (pdf) of the error between the real trajectories and the fitted ones, d = (X real -X f it ) 2

  Figure 4.23 -Snapshots of different acquisition of the tracer particles.

  ); 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 & 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 border effects become significant.

Figure 1 .

 1 Figure1. The fluid volume fraction (i.e., the porosity) as a function of the radius of the bed. The results are compared with those obtained inMueller (1992) andBenenati & Brosilow (1962).

Figure 2 .

 2 Figure 2. The 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. al Hill et al. (2001). Inset: the relative error shown as a percentage.

F 2 )Figure 4 .

 24 Figure 4. Example of the pressure obtained at different heights z i 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.

Figure 5 .

 5 Figure 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.

Figure 6 .

 6 Figure 6. Pressure field at two different height of the bed, for D/d = 10, ε = 0.518 and Re = 200. Results from Direct Numerical Simulations.

Figure 7 .

 7 Figure 7. Pressure drop measured experimentally for nine different random arrangements of spheres with D/d = 10. 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.

Figure 8 .

 8 Figure 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 DNSs.

Figure 9 .Figure 10 .

 910 Figure 9. Standard deviations of ⟨G p (r = R)⟩(θ) compared to its mean over all angles, ⟨G p (r = R)⟩ θ .

Figure 11 .

 11 Figure 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.

Figure 12 .

 12 Figure 12. Pressure drop for cases with D/d = [4.25, 5, 5.04], compared with the Ergun model for (D/d → ∞).

Figure 13

 13 Figure 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 sphere for D/d = 5, 8 and 10.

Figure 6 . 2 -

 62 Figure 6.2 -Velocity map of an experiment with the Reynolds number based on the superficial velocity Re(U Q ) = 211. The local velocity u y though can however become much higher in some pores.

Figure 6 .

 6 Figure 6.2 shows the average velocity map for a particular experiment with Re = 211.It can be seen that the Reynolds number based on the measured local fluid velocity u y can be several times higher than the "macroscopic" Reynolds number Re(U Q ) based on the superficial velocity U Q . This evidences the fact that not only there are very high local values of velocity, but this also probably generates high values of local acceleration as well, corresponding to more complex hydrodynamics at the pore-scale compared to the global scale of the bed.

Figure 6 . 3 -

 63  of the transversal and stream-wise components of the velocity, respectively. They are expressed in such a way so that they can be compared to a Gaussian distribution of mean zero and standard deviation equal to one, which is typically observed in turbulent flows.

Figure 6 .

 6  of the transversal and stream-wise components of the acceleration, respectively. They are expressed in such a way so that they can be compared to a Gaussian distribution of mean zero and standard deviation equal to one. Both distributions present exponential tails, which are observed in turbulent flows where extreme values of acceleration are probable.

Figure 6 . 6 -

 66 Figure 6.6 -Top: The area under the curve of the auto-correlation function as a function of dr.The correlation length is calculated when this quantity reaches a plateau. Bottom: Correlation lengths for the transversal and axial components. They vary with the Reynolds number and seem to converge towards a constant value so as to reach a certain level of isotropy.

Figure 6 .

 6 Figure 6.10 -Velocity map obtained from the simulations for D/D h = 6.45 for a Reynolds number based on the superficial velocity Re(U Q ) = 160. As seen in the experiments, local values of the velocity can get much higher than U Q .

Figure 6 .

 6 Figure 6.11 -Left: Probability density function (pdf) of the transversal components of the velocity u x and u z , and its radial component u r = u 2 x + u 2 z , for two different Reynolds numbers, and they are all compared to a gaussian distribution. Right: Pdf of the stream-wise component of the velocity for two different Reynolds numbers.

Figure 6 .

 6 Figure 6.12 -Probability density functions for all components of the acceleration field for two different Reynolds numbers, compared to a gaussian distribution. They all have exponential tails, indicating that extreme values of accelerations are more probable.

Figure 6 .

 6 Figure 6.13 -The same velocity field shown in figure 6.10 interpolated into a cartesian grid.

Figure 6 .

 6 Figure 6.16 -Second-order structure functions for two different Reynolds numbers. The inertialrange is (almost) non-existent for S2 y (dy), whereas it appears fro the other component.

Figure 6 .

 6 Figure 6.18 -Total velocity-acceleration structure function (all three components of the velocity and acceleration are taken into account) as a function of two different increments: one at the stream-wise direction and the other at the transversal one. It is evident that most of the negative contribution comes from the perpendicular component.

Figure 7 . 1 -

 71 Figure 7.1 -Force probability density function for D/d = 10, ε = 0.5, Re = 100 and 200 (Left), and ε = 0.5, D/d = 5 and Re = 200 (Right).

Figure 7 . 2 -Figure 7 . 3 -

 7273 Figure 7.2 -Force probability density function for D/d = 10, Re = 200 and ε = 0.5 and 0.9.

Figure 7 .

 7 Figure 7.4 -Diagonal components of the Reynolds stress tensor as a function of the Reynolds number for two beds with scale separation D/d = 10 and porosities ε = 0.5 and 0.9.

Figure 7 . 7 -

 77 Figure 7.7 -The stream-wise component of the Reynolds stress tensor u z u z calculated as a function of r for fixed D/d = 10 and two different Reynolds numbers (a)), and for fixed Reynolds number (Re = 200) and D/d = 10 and 5 (b)).

Figure 7 .

 7 7 shows the stream-wise component of the Reynolds stress tensor over three different ensembles with scale separation D/d = 10 and porosity ε = 0.5 for Re = 100 and 200 averaged for different concentric cylinders with intervals d/10, u z u z (r).

Figure 7 .

 7 Figure 7.11 -Forces radial distribution for different porosities ε = 0.9 and 0.5 at fixed D/d = 10 and Re = 200.

  table 1 they also present comprehensive list summarizing the different models for different D/d ratios represented in the figure.

  .10 the cumulants K of p and p 2 are shown respectively, which are by definition p = ˆp p

min P (p )p dp and p 2 = ˆp p min P (p )p 2 dp ,

  Revisiting the influence of confining wall on the pressure drop in fixed beds Florencia Falkinhoff 12 , Jean-Lou Pierson 1 , Lionel Gamet 1 , Mickael Bourgoin 2 and Romain Volk 2 1 IFP Energies Nouvelles, 69360, Solaize, France 2 Ens de Lyon, CNRS, Laboratoire de physique, F-69342 Lyon, France (Received xx; revised xx; accepted xx)

	To be submitted in Jan-Feb 2023	1

Table 1 .

 1 Parameters of the simulations. ε = V f luid /V total : fluid volume fraction. Re = U d/ν Reynolds number based on the superficial velocity U , the diameter of the particles d and the kinematic viscosity ν.

	D/d H/d 10.15 6.35 0.485 ε 8.03 5.02 0.490 5.13 3.85 0.518	Number of cells ρ 68.8M 1 0.53-7.89 1e-4 20 -300 U ν Re 58.9M 1 0.42 -6.20 1e-4 20 -300 45.3M 1 0.26-3.83 1e-4 20 -300

Table 2 .

 2 Range of parameters studied experimentally and numerically.

Table 3 .

 3 Experimental and numerical values of the slope and intersection of the G p (1/Re) linear fits, α and β respectively.

	dns

Table 6 .

 6 2 -Different Reynolds numbers based on different scales. We pay special attention to R λ , which is based on the inertial range as its characteristic scale.

	15Re L

Table 7 .

 7 1 -Simulations carried out in this chapter. The number of ensembles corresponds to the different arrangements studied for identical D/d, H/d and Re conditions, and the Total number of particles corresponds to the sum of all the spherical beads in all the ensembles.

	D/d H/d	ε	Re	# of ensembles Total # of particles
	10.15 6.35 0.90	200	9	900
	10.15 6.35 0.485 {20;60;100;150;200}	3	1500
	5.13 3.85 0.518	200	5	360

Our definition of the Reynolds number is equivalent to ReD in[START_REF] Wood | Modeling Turbulent Flows in Porous Media[END_REF].

The Reynolds number inHill et al. (2001) is defined in terms of the particle radius, whereas we define it in terms of the particle diameter here.

This problem will not be addressed in this thesis, but it is highly interesting to observe the frequency modes of the bed, how the fluidization depends on a critical Re number and on the density of the beads. A few simple tests were carried out to this end, out of scientific curiosity, and a jump at the pressure drop was observed at the Re at

It is worth noting that even though we are measuring at a fixed z-plane the width of the laser sheet is 4mm, which is comparable to the typical size of the pores, which ranges between 1 and 5mm. We are thus computing the 3D trajectories within the sheet, which has a width that is adequate to resolve the pore-scale.

Following an approximation done in[START_REF] Andreotti | Granular media: between fluid and solid[END_REF], the characteristic pore size s can be calculated by the following relation: φ/φmax = 1/(1 + s/D h )

, where φmax is the maximum solid fraction. Using this, the characteristic pore size is approximately 16mm in our case.

C2 is usually defined in terms of the Fourier transform of the second-order structure function, which is the energy spectrum E(k) = C k 2/3 k -5/3 , and C2 and C k are analytically related[START_REF] Sreenivasan | Multiphase flows: Rich physics, challenging theory, and big simulations[END_REF], thus C2 is also called the Kolmogorov constant.

It is worth noting that these scales validate our choice of tracer particles, which should typically be smaller than the Kolmogorov scale in order to be treated as tracer particles, and their diameter is 30µm. This also explains why we do not see the dr 2 regime in figure6.7: the smallest scales that we study are higher than η.

The difficulties of achieving meshes with lower porosities have already been discussed in chapter 3.

It is worth noting though that we are not expecting to get the same results as before, first of all because of the previous discussion and secondly because equation(6.8) is derived for an r-sphere[START_REF] Hill | Opportunities for use of exact statistical equations[END_REF], while in this case we are looking at the different components

We will not take into account the time-variable t because we are dealing with steady-state flows. The timedependency in the equation is trivial by considering the function q(y, t)

[START_REF] Singh | Effect of nozzle geometry on critical-subcritical flow transitions[END_REF] The notation presented here is such that uu f = u ⊗ u f is a tensorial product, and thus ∇ • uu f = is a vector

As it has been mentioned in chapter 3 meshing the arrangements is not trivial and the time that it takes to generate a particular mesh can take weeks.

We could further explore the border effects by studying for instance D/d = 20 or a tri-periodic arrangement of spheres where no border effects are present.

Although are studying it, the term u r u r f does not appear in equation(7.6)

Here we have labelled the "fixed" porosity at ε = 0.5, but in fact it slightly varies with D/d, as it is expected from the confinement effects (see table

7.1).

Acknowledgements

The calibration is done at different distances from the cameras, so that we have enough points to do the ray crossing in the PTV (this will be detailed in the next section). For each picture of the planes taken, we identify the square (marked in red in fig. 4.15) and the circles in the target. Using the square as reference, we define a coordinate system that starts at (0, 0) for each plane, and each circle is then labelled at each plane. This is shown in figure 4.16. Because we know the physical units (5mm between each point), we can compute the transformation necessary to convert the pixels into meters. The transformation is then used spherical beads with diameter d = [7.938, 4.762, 3.969] mm and density 8 g 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 whatever the 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 hours 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 MVM-020-QA from Bronkhorst (Netherlands). This provides a direct measurement of the mean velocity U = 4Q/(πD 2 ) that is used to define the Reynolds number Re = U d/ν 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 can be 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 highprecision piezo-resistive sensors with a sensibility of 150 mV/bar and a linear deviation of 0.05%/bar. In order to increase the signal-to-noise ratio, the signal is amplified The oscillating pattern is due to the geometry of the bed and shows a period of approximately one particle diameter d, and higher velocity values are observed near the walls. This can be explained by looking at the porosity profile shown in figure 7.5 b) 7 : the presence of the walls generates a higher fluid fraction near them, thus offering a less-resistant path for the fluid to flow through, thus generating higher values of local velocity. The same type of behaviour is observed for D/d = 5 (not shown here). These velocity profiles show that the term involving ∂ r u z f is not zero when we calculate it at these scales, as its local average varies with r as opposed to the global average. In order to consider the wall effects, we now compare the same terms with two different scale separations D/d = 10 and 5, and fixed porosity ε = 0.5 and Reynolds number Re = 200. The average for the beds with D/d = 5 is done over 5 different ensembles of particles and the results are shown in figure 7.9b), this time normalized by the reactor radius R so that they can be easily visualized. The border effects are once again evidenced by the increase of the streamwise velocity fluctuations, with once again a 33% increase near the walls when compared to the bulk average and showing the inhomogeneities when the there is a higher fluid volume fraction.

Lastly, we study the porosity effects for fixed scale separation D/d = 10 and Reynolds number Re = 200, and two different porosities ε = 0.5 and 0.9. The profile for ε = 0.9 was averaged over 9 different configurations. As shown in figure 7.9, the profile decreases with increasing porosity, which is expected from the bulk-averaged values presented in figure 7.4. There are less inhomogeneities when ε = 0.9, which is in agreement with what was observed in the velocity profile shown in figure 7.6. There are also less solid-fluid interactions and low-porosity regions for fluctuations to develop, and thus the system is more homogeneous, and we do not observe the same increase near the walls.

Force distribution

The other remaining closure term is the one that accounts for the sphere-fluid inter-phase force (we have already seen that the first moments of the forces do not affect the dynamics). Even though the statistics might not be fully converged (especially the case with D/d = 5), we can get some clues for the force distribution.

Summary and conclusions

Porous media are present in a wide range of phenomena, from industrial applications such as packed bed reactors to the medicine field. We are particularly motivated by the use of fixed beds in the context of AA-CAES technologies as a thermal energy storage unit, and our objective was to study the hydrodynamics in fixed beds of spherical particles that are randomly distributed inside a cylindrical container (e.g. a reactor), so that border effects are not negligible in this case.

Aside from our current motivation, fixed beds, as porous media, are a multi-scale system. That means that the hydrodynamics at the smallest scales of the system can be very different from that at the largest scales. In order to bridge this scale hierarchy upscale techniques are usually used such as volume averaging, that link the small-scale fluctuations with the large-scale dynamics. These techniques are typically done at a mesoscale, which should in principle be much larger than the pore scale and much smaller than the global scale.

The spheres distribution along with the wall presence produce tortuous paths and local changes in porosity, making the flows that saturate fixed beds incredibly complex, even if the flow is laminar. Moreover, there are fluid-solid interactions that can also affect the dynamics of the flow. Because of the complexity of the flows and the need to have all the information possible on the hydrodynamics, we studied the behavior of the confined flow at all three scales: the macro, the micro and the meso.

We did this by computing numerical simulations and doing two different experimental campaigns. Both the experimental and numerical methods allowed us to study hydrodynamical effects of the three non-dimensional parameters involved in the system: the Reynolds number Re, the porosity ε and the sphere-to-reactor diameter ratio D/d. In order to do this we had to design and construct two different experimental setups from scratch, to characterize the instruments used and we developed a Particle Tracking Velocimetry (PTV) code adapted for refractive index-matching PTV.

The simulations were all done using the open-source CFD software OpenFOAM. The packed beds were generated using the discrete element method software Grains3D, that allowed us to reach porosities down to 50%. After validating the OpenFOAM methodology we computed periodic arrangements of spherical particles with meshes that have between 40 to 90 million cells so that we could resolve all the small gaps. In total 18 different configurations with different D/d and ε where used to obtain physical results, and for each bed we studied the flow for different Reynolds numbers. The fluid solid forces seem to beevenly distributed at these scales