E Butaye

A Toutant
email: adrien.toutant@univ-perp.fr

S Mer

F Bataille

Development of Particle Resolved -Subgrid Corrected Simulations: Hydrodynamic force calculation and flow sub-resolution corrections

Keywords: Fluid -particles flows, DEM, Particle Resolved -Direct Numerical Simulation, Hydrodynamic force calculation, Fully resolved particles, Viscous penalization

published or not. The documents may come

Introduction

In the context of solar power generation, fluidized beds are used as a heat transfer fluid within verticals tubular receivers in solar towers power plants.

Under certain conditions of velocity and pressure, the flow attributes properties to solid particles that are specific to fluids (mixing viscosity, bubbly flow, waves. . . ). Over the last decades, numerous experimental studies in the literature have described This is an author's version, please cite: https://doi.org/10.1016/j.compfluid.2023.106071 fluidized beds behavior [START_REF] Kunii | Chemical reaction engineering and research and development of gas solid systems[END_REF][START_REF] Bi | Flow regime diagrams for gas-solid fluidization and upward transport[END_REF][START_REF] Gauthier | Influence of the particle size distribution of powders on the velocities of minimum and complete fluidization[END_REF][START_REF] Boissiere | Experimental study of gas-particle dense suspension upward flow for application as a new heat transfer and storage fluid[END_REF][START_REF] Gueguen | Gas-solid flow in a fluidizedparticle tubular solar receiver: Off-sun experimental flow regimes characterization[END_REF]. Flow regimes transitions are now well characterized in terms of particles shape [START_REF] Geldart | Types of gas fluidization[END_REF] and hydrodynamics [START_REF] Gueguen | Gas-solid flow in a fluidizedparticle tubular solar receiver: Off-sun experimental flow regimes characterization[END_REF]. However, local measurements of flow properties remain very difficult without invasive methods. Numerical studies appear therefore essential to fully characterize the flow, see [START_REF] Deen | Review of discrete particle modeling of fluidized beds[END_REF][START_REF] Van Der Hoef | Multiscale modeling of gas-fluidized beds[END_REF][START_REF] Alobaid | Progress in cfd simulations of fluidized beds for chemical and energy process engineering[END_REF] for a review of numerical methods and CFD simulations of fluidized beds.

Depending on the precision required, numerical modeling can be conducted at different scales. On large industrial geometry, Euler-Euler methods can be found. Averaged Navier-Stokes equations are solved in the fluid phase, and solid particles are represented with a locally averaged transport equation that ensures the solid behavior. A solid volume fraction is defined in each cell of the computa-tional domain and is transported at each time step through a continuity equation [START_REF] Ishii | Thermo-fluid dynamics of two-phase flow[END_REF][START_REF] Benoit | Three-dimensional numerical simulation of upflow bubbling fluidized bed in opaque tube under high flux solar heating[END_REF]. At laboratory scale, Euler-Lagrange methods represent individually each particle [START_REF] Zhou | Cfd-dem study of gas-solid flow characteristics in a fluidized bed with different diameter of coarse particles[END_REF]. Interactions between the flow and particles are not resolved but considered by means of drag, lubrication, added mass, history, and lift forces. The expression of such forces rely on correlation well referenced in the literature [START_REF] Tenneti | Drag law for monodisperse gas-solid systems using particle-resolved direct numerical simulation of flow past fixed assemblies of spheres[END_REF][START_REF] Zastawany | Derivation of drag and lift force and torque coefficients for non-spherical particles in flows[END_REF][START_REF] Beestra | Drag force of intermediate reynolds number flow past mono-and bidisperse arrays of spheres[END_REF][START_REF] Simonin | Drag force modelling in dilute to dense particle-laden flows with mono-disperse or binary mixture of solid particles[END_REF]. Particles are either simply advected by the fluid (one-way coupling) or collision laws are used to model interactions between them (four-way coupling) [START_REF] Vowinckel | Incorporating grain-scale processes in macroscopic sediment transport models[END_REF]. In this method, the solving of the fluid flow and the treatment of particles is known as Computational Fluid Dynamics -Discrete Element Method (CFD-DEM). Navier-Stokes equations are solved for the fluid phase and Newton-Euler equations coupled with collision models [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] are solved for the solid phase [START_REF] Tsuji | Direct particle sim-ulation of two-dimensional fluidized bed[END_REF]. On a smaller scale, Particle Resolved -Direct Numerical Simulations (PR-DNS) intend to fully characterize interactions between the fluid and particles. In this framework, the solving of the flow is identical to CFD-DEM, but forces exerted by the particles to the fluid are fully resolved. PR-DNS is widely used in the literature to study local interactions between the fluid and the particles, and especially to develop precise drag force models to inform Euler-Lagrange methods [START_REF] Buettner | Fluid-particle drag force measurements from particle-resolved cfd simulations of flow past random arrays of ellipsoidal particles[END_REF][START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF][START_REF] Duan | Fluid-particle drag and particle-particle drag in low-reynolds-number bidisperse gas-solid suspensions[END_REF][START_REF] Tenneti | Particle-resolved direct numerical simulation for gas-solid flow model development[END_REF].

At least 5 Eulerian meshes in the boundary layer are required to fully solve fluid-particle interactions in PR-DNS [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF]. The higher the Reynolds number, the smaller the dynamic boundary layer and thus the higher the grid resolution must be. However, the term PR-DNS is commonly employed when about 10 Eulerian meshes represent the particle. In most cases, the boundary layer is only partially resolved and the use of the term PR-DNS seems no longer relevant. PR-DNS are thus, with fully resolved boundary layers, limited to reproduce simple experimental case [START_REF] Cate | Particle imaging velocimetry experiments and latticeboltzmann simulations on a single sphere settling under gravity[END_REF][START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF]. Nonetheless, these PR-DNS simulations are used to provide models for Euler-Euler simulations. A gap lies between the very simple cases in PR-DNS and the actual application in industrial areas. The present article aims to address part of the aforementioned problem by enlightening fluid-particle interactions. As a first step, this paper focuses on a single particle in sedimentation. Prospects for using this methodology in more complex flows, such as fluidized beds, are drawn in the conclusion.

One propose to distinguish two cases in Particle Resolved Simulations (PRS). Firstly, simulations where boundary layers are fully resolved, known as PR-DNS. In this case, velocity gradients and pressure, interpolated from the Eulerian mesh, are fully resolved, resulting in an accurate computation of hydrodynamic forces exerted by the fluid on the particles. Then, an intermediate scale, named Particle Resolved -Subgrid Corrected Simulations (PR-SCS) where boundary layers are not resolved (see 1). In the latter, the Eulerian mesh is too coarse to correctly capture the velocity gradients and pressure at the interface between the fluid and the solid particles. PR-SCS scale is adapted for particle resolved simulations with a dozen meshes per diameter, as in [START_REF] Ozel | Particle resolved direct numerical simulation of a liquid-solid fluidized bed: Comparison with experimental data[END_REF][START_REF] Hamidi | Assessment of a coupled vof-front tracking/dem method for simulating fluid-particles flows[END_REF][START_REF] Luo | Particle-resolved direct numerical simulation of gas-solid dynamics in experimental fluidized beds[END_REF]. The unresolved part of fluid-particle interactions is modeled by a mesh-dependent correction developed with PR-DNS. Hereinafter, the term Particle Resolved Simulations (PRS) will refer to simulations where boundary layers, and thus the velocity gradients and the pressure at the interface between the fluid and the particles, are not fully resolved, and no correction is used to model the unresolved part of fluid-particle interactions. A first step to this work consists of studying the flow around a single particle at low Reynolds numbers, as in Stokes configuration [START_REF] Stokes | On the effect of the internal friction of fluids on the motion of pendulums[END_REF]. Thus, the analytical solution of the flow is well known and commonly referred to in the literature for validation. The coarser the mesh, the less captured the velocity gradients and the pressure, leading to sub-resolved hydrodynamic forces. For such resolution, the terminal velocity of the particle is no more in agreement with the theory. To overcome this, a correction term F c is added to the Navier-Stokes equation and discretized over the particle volume to correct the sub-resolved friction F P R-SCS . This procedure is very similar to the law of the wall in LES methods for turbulence modeling. However, one should point out that the proposed correction is essentially numerical, induced by the poorly velocity gradients and pressure capture. In this way, it differs from LES methods which aim to model subgrid interactions. The hydrodynamic force can thus be expressed as

F h = F P R-SCS + F c .
The paper is organized as follows. First, the modeling approach is presented in 2. The case of Stokes sedimentation is investigated in 3. In 4 a mesh dependent correction is built for PR-SCS. The strengths and limitations of the approach are assessed in 5. Finally, in 6, conclusions and perspectives to the present work are drawn.

• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • •
• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • •• • • •• • •• • •• •• • •• •• •• • ••• • • • •• • •• • • •• • • •• •• • •• • •• ••• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • •• • • • • • • • • • • • • • • • • • • • • • • • • •• • • • • •• • • • • • • • • • • • • • • • • Euler -

Modeling

Numerical modeling of PR-DNS flows has been extensively investigated. Two main approaches can be identified: body-fitted mesh [START_REF] Hu | Direct simulation of fluid particle motions[END_REF][START_REF] Hu | Direct numerical simulations of fluid-solid systems using the arbitrary lagrangian-eulerian technique[END_REF][START_REF] Nomura | An arbitrary lagrangianeulerian finite element method for interaction of fluid and a rigid body[END_REF] and fictitious domain [START_REF] Glowinsky | A fictitious domain approach to the direct numerical simulation of incompressible viscous flow past moving rigid bodies: Application to particulate flow[END_REF][START_REF] Khadra | Fictitious domain approach for numerical modelling of navier-stokes equation[END_REF][START_REF] Scardovelli | Direct numerical simulation of free-surface and interfacial flow[END_REF]. While the body-fitted method describes interfaces accurately, it requires remeshing while moving the interface, which results in a very high computational cost. This method is therefore not adapted for particle resolved simulation of fluidized beds, but very interesting for fixed arrays of spheres [START_REF] Massol | Simulations numériques d'écoulements à travers des réseaux fixes de spheres monodisperses et bidisperses, pour des nombres de reynolds modérés[END_REF]. The difficulty of representing the interface can be overcome through the use of the fictitious domain method. The mesh of interest is embedded in a mesh of simpler geometry. The inside of the particles is now represented, hence the name fictitious domain [START_REF] Khadra | Fictitious domain approach for numerical modelling of navier-stokes equation[END_REF]. Fictitious domain does not require remeshing. The grid used is fixed, and the interfaces are represented using different techniques, developed hereinafter. Navier-Stokes equations are solved on all the computational domain, including the inside of particles. A numerical treatment is required to enforce the solid body motion in this part of the domain. Two main methods are encountered in the literature to do so, namely the Immersed Boundary Methods (IBM) [START_REF] Peskin | Flow patterns around heart valves: a numerical method[END_REF][START_REF] Uhlmann | An immersed boundary method with di-rect forcing for the simulation of particulate flows[END_REF] and viscous penalization [START_REF] Ritz | A numerical continuous model for the hydrodynamics of fluid particle systems[END_REF][START_REF] Vincent | A lagrangian vof tensorial penalty method for the dns of resolved particleladen flows[END_REF]. In IBM, boundary conditions are explicitly imposed over the particle surface, whereas viscous penalization implicitly enforce fluid constraints on the interface. One should note that in IBM, even though the inside of the particle is represented by a fictitious fluid, it is not considered for the solving of the flow. The location of the interface can be achieved with different techniques. The rate of presence of solid in each Eulerian mesh is used by Volume-Of-Fluid (VOF) methods [START_REF] Hirt | Volume of fluid (vof) method for the dynamics of free boundaries[END_REF]. A signed distance to the interface is calculated for Level-Set methods [START_REF] Osher | Fronts propagating with curvature-dependant speed: Algorithms based on hamilton-jacobi formulations[END_REF]. Finally, the Lagrangian markers located on the interface are moved according to a transport equation with Front-Tracking method [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF].

Particle resolved numerical approach

The resolution method is based on the onefluid VOF/Front-Tracking approach implemented in TrioCFD software. Initially developed for gasliquid flows, Hamidi et al. extended the technique to fluid-solid flows [START_REF] Hamidi | Assessment of a coupled vof-front tracking/dem method for simulating fluid-particles flows[END_REF]. In this framework, only one set of equations is solved. The fluid properties vary at the interface according to a phase indicator function.

The fluid is considered incompressible in each phase. Navier-Stokes equations are given by:

∇ • u = 0 (1) ρ ∂u ∂t + (u • ∇) u = -∇p + ρg + ∇ • µ ∇u + ∇ T u ( 2
)
where u is the velocity, p the pressure, t the time, g the gravity vector, ρ the density, µ the viscosity. Velocity and pressure are solved on a staggered grid (see 2). The interface is located with the Front-Tracking method [START_REF] Unverdi | A front-tracking method for viscous, incompressible, multi-fluid flows[END_REF][START_REF] Tryggvason | A fronttracking method for the computations of multiphase flow[END_REF]. Particles are represented by a mobile mesh overlying on the fixed Eulerian mesh. This mesh is represented by the socalled Lagrangian markers, uniformly distributed over the surface. The method has the advantage of accurately describing particles of complex geometry. The interface is advected by a transport equation:

∂χ ∂t + u • ∇χ = 0 ( 3
)
where χ is the phase indicator function. Numerically, the procedure is the following. First, a trilinear interpolation (bilinear in two dimensions) of the Eulerian velocity on each Lagrangian marker is realized (see 2). Then, an average velocity is computed at the gravity center of Lagrangian facets. A Lagrangian facet is defined by three markers in three dimensions (two markers in two dimensions).

The overall average velocity of the particle p is given by:

u p = 1 S p facets S i u i (4)
where S p is the particle surface, S i and u i the surface and velocity of the facet i. Finally, the markers are advected with the velocity u p :

x t+1 i = x t i + u p ∆t (5)
where x i represent the position of the Lagrangian marker i and ∆t the time step. In this way, the particle does not deform over time and the mass conservation is satisfied. The local vorticity can initiate a rotation of the particle. However, for a spherical particle, the projection of its shape onto the Eulerian mesh will remain unchanged for any rotation about its gravity center. As Lagrangian markers are not material points but only describe the particle surface, a pure rotation of the Lagrangian markers is not required. For non-spherical particles, the current method would require some developments. Unlike gas-liquid flows, no remeshing is performed during the calculation. Then, χ is integrated on Eulerian cells to compute the fluid volume fraction. This integrated quantity divided by the mesh volume, noted as I, equals 1 in fully fluid cells and 0 in fully solid cells. In diphasic cells, the value of I lies between 0 and 1. Density and viscosity vary at the interface according to the following equations:

ρ = I ρ f + (1 -I) ρ s (6) µ = µ f µ s (1 -I) ρ f + I ρ s (7)
The subscripts f and s denote, respectively, the fluid and the solid. The precision and mesh convergence of the arithmetic model for the density and the harmonic model for the viscosity is demonstrated in [START_REF] Ritz | A numerical continuous model for the hydrodynamics of fluid particle systems[END_REF] and later validated in [START_REF] Toutant | Modélisation physique des interactions entre interfaces et turbulences[END_REF] for TrioCFD Front-Tracking method.

+ • • • • • • • • • • • • u s i • • • • • • • • • U s i V s i
Figure 2: Interpolation of Eulerian velocity on Lagrangian markers. u s i represents the velocity of the i-eme Lagrangian marker. • : pressure nodes, ▲ ▶ : velocity nodes, •: Lagrangian marker.

Discretization schemes and solvers

The mass and momentum conservation equations are solved with a first order Euler Explicit scheme with an implicit treatment for diffusion. A secondorder centered spatial discretization is employed for the convection and diffusion terms. The Poisson equation is solved with a prediction-correction algorithm. For this purpose, a conjugated gradient method with a symmetric successive overrelaxation (ssor) preconditioner is used. Finally, the Lagrangian markers are advected with a firstorder temporal scheme.

Hydrodynamic force computation

Although the number of PRS studies is very large, few of them specifically describe the hydrodynamic force computation method. One can cite [START_REF] Zastawany | Derivation of drag and lift force and torque coefficients for non-spherical particles in flows[END_REF][START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF][START_REF] Santarelli | Budget analysis of the turbulent kinetic energy for bubbly flow in a vertical channel[END_REF][START_REF] Chadil | Novel method to compute drag force and heat transfers for motions around spheres[END_REF]. In all these methods, pressure and viscous contributions are computed independently of each other. A local force is first computed at the gravity center of each Lagrangian facet and weighted by its surface (see 3). The accuracy of the computed hydrodynamic force then depends on the ability to compute the local velocity gradients and pressure at the gravity center of each Lagrangian facet. The sum of all the local components results in the global force exerted by the fluid on the particle. Zastawany [START_REF] Zastawany | Derivation of drag and lift force and torque coefficients for non-spherical particles in flows[END_REF] used a second-order extrapolation to compute the pressure on the surface and a second order Taylor series expansion to compute the velocity tensor at the interface. In this framework, diphasic cells are used to compute the velocity gradients which limits the accuracy of the viscous force estimation acting on the particle. Chadil et al. [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF] improved this method by using only fluid velocity to compute the viscous term. A third order Taylor interpolation and Lagrangian extrapolation of the same order are employed. Later, the method was enhanced by the use of Aslam extensions, as explained in [START_REF] Chadil | Novel method to compute drag force and heat transfers for motions around spheres[END_REF]. In all these studies, all the components of the viscous tensor were computed. In the method presented below, only radial derivatives of the velocity are required to compute the friction force, owing to the non-deformability of the solid interface. The calculation cost of the method is thus reduced and the precision improved, since the zero theoretical terms involved in Eq. 11 are not considered by means of the simplification detailed in Eq. 13.

The hydrodynamic force is the sum of the friction and pressure forces. It is written as follows: [START_REF] Van Der Hoef | Multiscale modeling of gas-fluidized beds[END_REF] with S p , the surface of the particle. This equation can be discretized on the surface of the particle as follows:

F h = Sp -pI 3 • n + µ f (∇u + ∇ T u) • n dS
F h = N k -p k S k I 3 • n k + µ f S k (∇u + ∇ T u) k • n k (9)
with N : the number of Lagrangian facets, p k I 3 and µ f (∇u + ∇ T u): the pressure and stress tensors at the gravity center of the Lagrangian facets.

The pressure and stress tensors are required at the center of gravity of the Lagrangian facets to compute the hydrodynamic force. However, the meshes crossed by the interface are, by definition, diphasic. The fluid velocity and pressure are therefore poorly defined at this location. It is necessary to move from the interface, along its normal, in order to extrapolate u and p from the closest single-phase fluid meshes. The method developed is inspired by that of [START_REF] Santarelli | Budget analysis of the turbulent kinetic energy for bubbly flow in a vertical channel[END_REF].

For the pressure, for each Lagrangian facet j, the computation steps are as follows:

1. Trilinear interpolation of the pressure at points P 1 and P 2 , located respectively, at a distance δ j and 2δ j from the interface (see 3). Let ∆x i , ∆y i and ∆z i , be the grid size of the Eulerian element i in which the facet j is located. n x , n y , n z , are the Cartesian components of the normal of the facet. δ n is a parameter for the user to define. It represents the distance to P 1 scaled by the mesh grid size for a Lagrangian facet which normal would coincide with one of the Cartesian axes. Then, δ j is computed as follows:

δ j = δ n n x ∆x i + n y ∆y i + n z ∆z i ( 10
)
In this way, in all geometric configurations, the element P 1 will be strictly fluid. For a uniform mesh of grid size ∆, the norm of δ j verifies: 1 ≤ δ j /∆ ≤ √ 3. This definition guarantees that the Eulerian mesh, in which the interpolation point P 1 is located, is strictly fluid (see point P 1 B in 3). Moreover, it adapts to a non-uniform mesh with, for example, near-wall refinement. 2. Linear extrapolation of the pressure at the gravity center of the Lagrangian facet j.

For the friction, the steps are as follows:

1. Interpolation of the velocity at points P 1 and P 2 in the same way as for the pressure, but with consideration of the location of the velocities. As shown in 2, the field is located at the gravity center of the faces of Eulerian elements. 2. Change of reference frame to express the velocity at points P 1 and P 2 in spherical coordinates (r, θ, Φ). 3. Computation of the stress at the gravity center of the facets:

df f = µ f (∇u + ∇ T u) r=Rp • n (11)
where R p is the radius of the particle. In spherical coordinates, the components of the tensor u that intervenes in the computation of (∇u + ∇ T u) r=Rp • n are:

    ∂Ur ∂r 1 r ∂Ur ∂θ -U θ r 1 r sin θ ∂Ur ∂Φ - U ϕ r ∂U θ ∂r - - ∂U Φ ∂r - -     r=Rp (12)
However, in the case of a solid particle, the non-deformation of the interface implies that the tangential components of the velocity are zero at the interface. This translates numerically into a constant distance between the Lagrangian markers. Thus, only radial derivatives need to be computed:

(∇u + ∇ T u) r=Rp • n =    2U r,r U θ,r -U θ r U ϕ,r - U ϕ r    r=Rp ( 13
) With a second order forward discretization scheme, the partial derivative of u along the normal is written as:

U α,r i | Rp = -U P2 α + 4U P1 α -3U i α 2δ ( 14
)
with α ∈ {r, θ, ϕ} and U i α the α component of the velocity of the facet i. 4. Change of reference frame to express the friction force in Cartesian coordinates.

The interest of this method lies in the fact that tangential derivatives of the velocity are not involved in the computation of the friction at the interface. Indeed, all derivatives involved in the projection along the interface normal of the stress tensor are only with respect to r. Thus, the only knowledge of the velocities at P 1 and P 2 allows solving the friction at the interface. Moreover, as the stress tensor is not interpolated at P 1 and P 2 , the associated computational cost is greatly reduced as well as the joint area -the number of virtual meshes that a processor must know during parallel operations. The default thickness of the joint area is 2, whereas the proposed method requires an extension to 4. However, it should be noted that it must be increased to 7 for the complete calculation of the stress tensor at P 1 and P 2 . This illustrates the memory cost reduction of the proposed method.

Stokes sedimentation

The objective of this section is to validate the force computation method and its implementation. The case of Stokes sedimentation is considered, for which the solution is known. The reference frame of the particle is considered. A reduced computation domain is defined where the theoretical solution is applied as a boundary condition. Fine meshes can then be used without incurring high computational costs. The hydrodynamic force exerted by the fluid on the particle is computed by the method described in the Sec. 2.3 and confronted with the theoretical solution.

• • • • • • • • • • • • • • • • • • × × × × × × P A 1 P A 2 P B

Simulations setup

A fictitious domain method is adopted in the current work, the flow inside the solid particle must be considered. In a numerical viewpoint, the solid particle is a highly viscous drop. The Stokes flow is solved in each phase. The closure is achieved through the Delahaye stress jump equation [START_REF] Delahaye | Jump conditions and entropy sources in two-phase systems. local instant formulation[END_REF].

In the reference frame of the particle, and in spherical coordinates, the velocity, and pressure fields can be derived analytically as shown in Tab. 3. In this set of equations, Φ µ = µp µ f is the viscosity ratio between the particle and the fluid, U ∞ the terminal velocity of the particle and R p the particle radius. A dynamic equilibrium between the weight of the particle, buoyancy and the hydrodynamic resultant force leads to the following expression of the settling velocity:

U ∞ = 2 3 gR 2 p µ f 1 + Φ µ 2 + 3Φ µ (ρ f -ρ p ) (15)

Flow solver assessment

The 3D computational domain is schematized in Fig. 8. Velocity and pressure given in Tab. 3 are Fluid flow

U r -U ∞ (1 - 2+3Φµ 1+Φµ Rp 2r + Φµ 1+Φµ R 3 p 2r 3 ) cos θ U θ U ∞ (1 - 2+3Φµ 1+Φµ Rp 4r - Φµ 1+Φµ R 3 p 4r 3 ) sin θ U φ 0 P -µ f U ∞ 2+3Φµ 1+Φµ Rp 2r 2 cos θ Inside the particle U r U∞ 2 1 1+Φµ 1 -r 2 R 2 p cos θ U θ -U∞ 2 1 1+Φµ 1 -2 r 2 R 2 p sin θ U φ 0 P µ p U ∞ 5 1+Φµ r R 2 p cos θ
Table 1: Theoretical velocity and pressure fields.

initially imposed on all the domain. A pressure condition is imposed on the top boundary, whereas velocity is imposed on the others. The exact solution of Stokes flow is thus applied on the boundaries, enabling a reduced computational domain to be used. The fluid is injected from the bottom of the domain at a prescribed velocity (see Tab. 3). The hydrodynamic force exerted by the fluid on the particle balance gravity. In other words, if the friction is well taken into account numerically, the particle must remain at rest in the calculation domain. Uz /U∞ (%) Figure 5: Mesh convergence of the relative settling velocity error, normalized by U ∞ , the theoretical settling velocity given by Eq. 18.

given in Tab. 4. The Reynolds and Stokes numbers are respectively 7.8 10 -2 and 8.7 10 -2 . These are defined as follows:

R e = 2R p ρ f U ∞ µ f ( 16
)
S t = 1 9 Re ρ p ρ f ( 17
)
ρ p ρ f µ p µ f D p g [kg.m -3 ] [kg.m -3 ] [P a.s] [P a.s] [m] [m.s -2 ]
10 4 10 3 3 3 10 -3 5.2 10 -5 10

Table 2: Physical parameters of Stokes configuration.

The relative settling velocity error is represented in Fig. 9. The fields obtained in the finest simulation ( D ∆x = 80) are represented in Fig. 10 and Fig. 11. As the referential frame of the simulation is the one of the particle, the theoretical velocity of the particle is zero. Above a grid resolution of 20 meshes per diameter, the relative settling velocity error is below 0.3 %.

Stokes hydrodynamic force

The objective of this section is to validate the force computation method and its implementation. The case of Stokes sedimentation is considered, for which the solution is known. The reference frame of the particle is considered. A reduced computation domain is defined where the theoretical solution is applied as a boundary condition. Fine meshes can then be used without incurring high computational costs. The hydrodynamic force exerted by the fluid on the particle is computed by the method described in 2.3 and confronted with the theoretical solution.

Simulations setup

A fictitious domain method is adopted in the current work, the flow inside the solid particle must be considered. In a numerical viewpoint, the solid particle is a highly viscous drop. The Stokes flow is solved in each phase. The closure is achieved through the Delahaye stress jump equation [START_REF] Delahaye | Jump conditions and entropy sources in two-phase systems. local instant formulation[END_REF].

In the reference frame of the particle, and in spherical coordinates, the velocity, and pressure fields can be derived analytically as shown in 3. In this set of equations, Φ µ = µp µ f is the viscosity ratio between the particle and the fluid, U ∞ the terminal velocity of the particle and R p the particle radius. A dynamic equilibrium between the weight of the particle, buoyancy and the hydrodynamic resultant force leads to the following expression of the settling velocity:

U ∞ = 2 3 gR 2 p µ f 1 + Φ µ 2 + 3Φ µ (ρ f -ρ p ) (18)

Flow solver assessment

The 3D computational domain is schematized in 8. Velocity and pressure given in 3 are initially imposed on all the domain. A pressure condition is imposed on the top boundary, whereas velocity is

Fluid flow U r -U ∞ (1 - 2+3Φµ 1+Φµ Rp 2r + Φµ 1+Φµ R 3 p 2r 3 ) cos θ U θ U ∞ (1 - 2+3Φµ 1+Φµ Rp 4r - Φµ 1+Φµ R 3 p 4r 3 ) sin θ U φ 0 P -µ f U ∞ 2+3Φµ 1+Φµ Rp 2r 2 cos θ Inside the particle U r U∞ 2 1 1+Φµ 1 -r 2 R 2 p cos θ U θ -U∞ 2 1 1+Φµ 1 -2 r 2 R 2 p sin θ U φ 0 P µ p U ∞ 5 1+Φµ r R 2 p cos θ
Table 3: Theoretical velocity and pressure fields. imposed on the others. The exact solution of Stokes flow is thus applied on the boundaries, enabling a reduced computational domain to be used. The fluid is injected from the bottom of the domain at a prescribed velocity (see 3). The hydrodynamic force exerted by the fluid on the particle balance gravity. In other words, if the friction is well taken into account numerically, the particle must remain at rest in the calculation domain. The physical parameters of the configuration are given in 4. The Reynolds and Stokes numbers are respectively 7.8 10 -2 and 8.7 10 -2 . These are defined as follows:

ρ p ρ f µ p µ f D p g [kg.m -3 ] [kg.m -3 ] [P a.s] [P a.s] [m] [m.s -2 ]
R e = 2R p ρ f U ∞ µ f ( 19
)
S t = 1 9 Re ρ p ρ f ( 20
)
The relative settling velocity error is represented in 9. The fields obtained in the finest simulation ( D ∆x = 80) are represented in ??. As the referential frame of the simulation is the one of the particle, the theoretical velocity of the particle is zero. Above a grid resolution of 20 meshes per diameter, the relative settling velocity error is below 0.3%.

Stokes hydrodynamic force

The validation of the forces computation method is achieved in two ways: (i) on a local scale: the pressure and friction forces calculated with the extrapolation method, on each Lagrangian facet, are confronted to their respective theoretical value. (ii) on a global scale: the calculated integral sum of the pressure and friction forces are compared to their theoretical value. In the following, the upper scripts th , dis and int stand respectively for theoretical, discretization and interpolation. Uz /U∞ (%) Figure 9: Mesh convergence of the relative settling velocity error, normalized by U ∞ , the theoretical settling velocity given by Eq. 18.

For each Lagrangian facet i, the theoretical local friction and pressure forces are written, in spherical coordinate, as follows:

dF f th ,i = - µ f U ∞ R p     2 cos θi (1+Φ i µ ) - 3Φ i µ sin θi 2(1+Φ i µ ) 0     dS i ( 21
)
dF p th ,i = - µ f U ∞ 2R p 1 + 3Φ i µ 1 + Φ i µ cos θ i I 3 • e i r dS i ( 22
)
The integral over the particle surface gives the contribution of the friction and pressure terms:

F th f = -4πµ f U ∞ R p u p ||u p || (23)
F th p = -2πµ f U ∞ R p u p ||u p || (24)
Therefore, the overall hydrodynamic force is given as:

F th h = -6πµ f U ∞ R p u p ||u p || (25)
Moreover, considering Φ µ ≫ 1, the norm of the hydrodynamic constraint is uniform over the particle surface and given by the following equation for each Lagrangian facet i:

dF th h,i dS = 3µ f U ∞ 2R p ( 26
)
This first validation step does not require PR-DNS to be performed, as only the methodology presented in 2.3 is evaluated. Instead, the theoretical Figure 13: Influence of the distance δ used in the interpolation method, to identify the points P 1 and P 2 (see Fig. 3), on the calculated force relative error. velocity and pressure fields are discretized on the Eulerian mesh and considered as an input for the interpolation of the pressure and friction forces from the fluid phase to the interface. A grid resolution of 80 meshes per diameter is adopted to reduce the discretization error. δ n is set equals to the mesh grid size as in [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF]. The surface pressure force is represented in 12. The left image presents the theoretical field obtained with Eq. 22 on the particle surface. In the middle, the equivalent surface pressure force -computed from the discretized theoretical pressure field, with the method described in 2.3 -is shown. Finally, on the right the relative error between the two fields is drawn. The maximum local relative error is less than 0.5%. With the same methodology, but utilizing the discretized theoretical velocity field as an input, a maximum relative error of 2.3% is calculated for the surface friction force. On a global scale, three assessments of the same force can be identified. Firstly, the theoretical hydrodynamic force given by the Eq. 25. Then, the sum of all the local contributions exerted on the Lagrangian facets:

|(F h -F th h )/F th h | (%)
+ • • • • • • • • • • • • - → U • • • • • • • • • V k
F th-dis h = Nb facets i dF h th ,i (27)
where dF h th ,i = dF f th ,i + dF p th ,i . In this case, dS i is the surface of the Lagrangian facet i. Finally, the force interpolated from the theoretical velocity and pressure fields, with the method described in 2.3, is noted as F th-int h . Therefore, the following relative errors are considered:

ϵ th_th-dis = ||F th-dis h || -||F th h || ||F th h || (28)
ϵ th_th-int = ||F th-int h || -||F th h || ||F th h || (29)
ϵ int_dis = ||F th-int h || -||F th-dis h || ||F th-dis h || (30)
ϵ th_th-dis characterizes the error in the hydrodynamic force caused by the discretization of the particle surface into Lagrangian facets. For a grid resolution of 80 meshes per diameter, the error over the surface, committed with the Lagrangian discretization, is equal to 0.03% leading to ϵ th_th-dis = 0, 06%. ϵ int_dis represents the error of the method, which refers to the errors committed during the interpolation and extrapolation processes of Eulerian variables on each Lagrangian facet. One should note that both F th-dis h and F th-int h are computed over the same discretized surface. Here, ϵ int_dis = 0.35%. However, the computed force should be compared to F th h which is the variable of interest. Therefore, the coefficient ϵ th_th-int is adopted to characterize the method, as it also considers the discretization error. In this configuration, it is equal to ϵ th_th-int = 0.29%. The implementation of the method is therefore validated.

A parametric study ensured that the distance δ n adopted according to [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF] was relevant to the developed method. 13 depicts the linear evolution of the error as a function of the distance chosen to interpolate the variables in P 1 . An absolute minimum is reached for a value equal to the mesh grid size which is in accordance with [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF].

A mesh convergence study of PRS is presented in 14. The computed hydrodynamic force is confronted with the theoretical value. PRS from 5 to 80 meshes per diameter were conducted. The convergence of the method is of order 2. For a grid resolution of 10 meshes per diameter, a relative error of 10% is computed. Indeed, the Eulerian grid resolution is too coarse to properly compute velocity gradients over the particles' surface. This highlights the need to develop corrections for PRS of fluid-particle flows using a grid resolution of a dozen meshes per particle's diameter.

Building a mesh-dependent correction for PR-SCS

A mesh-dependent correction is developed to account for the unresolved part of the friction and the pressure. The methodology is very similar to Breugems's work for the lubrication force [START_REF] Breugem | A combined soft-sphere collision / immersed boundary method for resolved simulations of particulate flows[END_REF]. A mesh-specific correction was developed for 2 resolutions to correct the lubrication force. In the present study, a generalized correction, therefore depending on the Eulerian grid resolution, is proposed for the hydrodynamic force.

Implementation of subgrid correction

The correlation between the force obtained in PRS, without any correction, and that obtained with PR-SCS is written as follows:

F PR-SCS = F PRS + F c ( 31
)
where F P RS is the hydrodynamic force captured with the PRS and F c is the correction force accounting for the unresolved part of the fluid-particle interactions. The correction is designed to correct the error committed on the velocity gradients and the pressure computation because the Eulerian mesh is not fine enough. The correlation assessment is described in plan::assess. It is developed to verify: lim Dp/∆x→∞ F PR-SCS = F PR-DNS [START_REF] Nomura | An arbitrary lagrangianeulerian finite element method for interaction of fluid and a rigid body[END_REF] lim

Dp/∆x→∞ F c = 0 ( 33
)
The correction is discretized over the volume of the particle and applied to the solid and diphasic faces of the Eulerian elements. For the sake of simplicity, a two-dimension representation of the discretization is given in 15. For an Eulerian face k, the component i of the modeled hydrodynamic force is applied as follows:

f c i k = I k F i c Nb faces k=0 I k V k ( 34
)
where

F i c = ||F c || U i ||u|| is the component i of F c , I k is the

Establishment of the mesh-dependent function

In the case of the Stokes sedimentation, the hydrodynamic force is well known. Thus, Eq. 31 can be rewritten as:

F PR-SCS = F Stokes = F PRS + F c (35)
With F Stokes given by the Eq. 25. The objective is to write F c as a function of the grid resolution. To do so, we assume -in the fashion of what did [START_REF] Breugem | A combined soft-sphere collision / immersed boundary method for resolved simulations of particulate flows[END_REF] for the lubrication force -that the correction is of the form:

F c = F Stokes h D p ∆x ( 36
)
where h

Dp ∆x is to be determined and depends on the Eulerian grid resolution. One should point out that the correction is Reynolds dependant as the theoretical Stokes force is a function of the Reynolds number

F Stokes = 1 2 ρ f S p U 2 ∞ 24
Re . Where S p = π * R 2 p . With results from the mesh convergence study, h can be computed as:

h D p ∆x = 1 - ||F PRS || ||F Stokes || (37)
With such a formulation, it is assumed that the force computed on the under-resolved mesh follows the correct direction, and only its magnitude must be corrected to retrieve the full hydrodynamic force. We will see later that for resolutions smaller than 5 meshes per diameter, this assumption can be defeated (see 18).

From Eq. 37, assuming h(x) = C1

x C 2 , one can be determined C 1 and C 2 with the least squares method as represented in Fig. 16. The obtained function for h can be written as follows, with x standing for Dp ∆x :

h(x) = 9.79 x 1.95 (38)
The correction is then applied in Navier-Stokes equations as explained in 4.1. Results of these PR-SCS are shown in 17. The maximum velocity error is less than 0.01% for a grid resolution of 5 meshes per diameter and less than 0.003% for the others, whereas the error without correction is higher than 4% for the coarser grid.

Although the particle settling velocity is slightly affected by this low-resolution limit, the problem is more marked on the hydrodynamic force estimate. (see [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF].

In this case, the overall hydrodynamic force obtained with PR-SCS overestimates the theoretical value. The resolved part of fluid-particle interac- tion is labeled as F P RS and the corrected part as F c in agreement with Eq. 31. One can note that the thinner the mesh, the more accurate the F P RS and the smaller the correction. The relative error on the hydrodynamic force in PR-SCS is less than 1.3% for a mesh resolution above meshes per diameter. It must be compared with a 10.6% error obtained without correction. One should point out that the same precision is obtained with the grid resolutions of 10 and 40 meshes per diameter, but the CPU cost of the simulation with only 10 meshes per diameter is 430 times lower. However, it is of 7% for the resolution of 5 meshes per di- ameter. The error in the computation of the force is more than ten times higher than the one committed on the velocity. This is explained by the interpolation method of the quantities of interest. The larger the mesh size, the further the points P1 and P2 are from the interface (see 3) and the less resolved the velocity gradients and pressure at the fluid-solid interface. In the present case, PR-SCS is therefore not adapted for PRS with less than 5 meshes in diameter. It has to be noted that the coarser grid resolution ( Dp ∆x = 5) was not used to build the correction given by Eq. 38.

Application to a dense fluid-particle flow

In a dense fluid-particle flow, particle-wall and particle-particle collisions occur. In the present method, the collision process is treated as in Discrete Element Method through a soft-sphere model developed by Hamidi et al. [START_REF] Hamidi | Assessment of a coupled vof-front tracking/dem method for simulating fluid-particles flows[END_REF]. When two particles undergo collision, the lubrication force increases as the particles get closer. When the collision occurs, the incompressible fluid film between the particles is of the order of the rugosity (approximately one micrometer). The capture of the lubrication force is therefore impossible in 3-dimensions with current numerical methods. However, to account for this phenomenon, lubrication models are widely employed in the literature [START_REF] Breugem | A combined soft-sphere collision / immersed boundary method for resolved simulations of particulate flows[END_REF][START_REF] De Motta | Numerical modelling of finite-size particle collisions in a viscous fluid[END_REF]. In the present framework, during a collision, the collision and lubrication forces predominate over the hydrodynamic forces -and thus our correction. Therefore, the proposed correction might be deactivated while the lubrication is considered.

Assessment of the Subgrid Corrected Simulations

Within the Stokes flow regime

The correlation is assessed through the study of four configurations presented in 5. These configurations are selected to verify Stokes regime hypothesis as Reynolds numbers are less than 10 -2 . Two density ratios of 10 and 3000 are chosen to represent, respectively, liquid-solid and gas-solid flows. The hydrodynamic force relative error is represented in [START_REF] Tsuji | Direct particle sim-ulation of two-dimensional fluidized bed[END_REF]. One can see that in these configurations, the hydrodynamic force relative error is independent of both the Reynolds number and the Stokes number. Furthermore, the PR-SCS force agrees almost perfectly with the theoretical Stokes force for all grid resolutions and all density ratios. Note that for small Re numbers the physics is independent of the density ratio, but the simulations are not.

case Re ρ p /ρ f St a 1.0 × 10 -3 10 1.1 × 10 -3 b 1.0 × 10 -3 3000 3.3 × 10 -1 c 1.0 × 10 -2 10 1.1 × 10 -2 d 1.0 × 10 -2 3000 3.3
Table 5: Stokes configurations for the correlation assessment.

Beyond the Stokes flow regime

To evaluate the limit of the proposed approach and especially the validity domain of the correction force proposed, simulations were performed for higher Re numbers and various density ratios (ranging from 2 to 3000). To evaluate the correction at higher Reynolds numbers, a hydrodynamic coefficient is computed and confronted with the Abraham correlation [START_REF] Abraham | Functional dependance of drag coefficient of a sphere on reynolds number[END_REF]. These coefficients are defined as follows:

C h = 2F h ρ f πR 2 p (U y -U ∞ ) 2 (39)
C dAbraham = 24 9.06 2

where F h is the hydrodynamic force in the settling direction, U y and U ∞ are, respectively, the residual velocity (numerical error) and theoretical settling velocity.

The simulation setup for this study is the same as described in 3. The imposed velocity on the boundaries is defined with the assumption that Stokes flow remains valid for these Re numbers. As the later assumption is defeated for Re > 0.1, one can expect the particle settling velocity to be nonzero. However, as it is the velocity difference between the fluid and the particle that is considered in Eq. 39, this enables the hydrodynamic coefficient to be evaluated. Results are presented in 20 for PRS and PR-SCS. The grid resolution for these simulations is 10 meshes per diameter. One can see that the correction improves the accuracy by 10% for Re numbers ranging from 10 -3 to 10 -1 , regardless of the density ratio. In this range, Abraham's correlation is valid and consistent with the Stokes theory. One can observe that the density ratio has no influence on the hydrodynamic coefficient in this creeping flow regime, as all the points are superimposed. When the Reynolds number increases, both PRS and PR-SCS deviate from the Abraham coefficient, following both the Stokes trend. In this range, the 10% improvement between the two methods seems to remain. It can be seen that, even outside its range of validity, the PR-SCS continues to closely follow the evolution predicted by the Stokes theory. This is attributed to the fact that the correction was initially designed for creeping flows. Furthermore, for Re ∼ 1, the density ratio affects the hydrodynamic coefficient, with a five percent dispersion of the different simulations. This dispersion increases for higher Re. The influence of the density ratio on the computation of the hydrodynamic force might be related to the difficulties of the Poisson equation solver in resolving strong density gradients at the interface. Finally, the correction developed in this paper is validated for Reynolds numbers up to 10 -1 and density ratios between 2 and 3000. The grid resolution of PRS must be of 10 meshes per diameter or higher to correctly capture the direction of the hydrodynamic force. To adapt the correction for higher Reynolds numbers, the structure of the correction (see Eq. 36) must be adapted to be consistent with Abraham's correlation. Furthermore, dependence with respect to the density ratio must be introduced.

In more complex flows such as fluidized beds, no theoretical solution exists to describe the velocity and pressure fields in the vicinity of the fluid-solid interface. It may also be difficult to access to purely fluid cells at point P 1 and P 2 because of the proximity of particles. To expand the correlation to higher Reynolds in fluidized beds, one key point has to be considered: how to derive a correction for higher Reynolds number, where Stokes formulation no longer holds. The targeted Reynolds numbers lie between 1 and 10 where the Abraham correlation remains valid. Thus, the Abraham drag coefficient (see Eq. 40) can be used to expand the validity of the correlation developed. At higher Reynolds number, in cases where the wake/history force cannot be neglected in front of other hydrodynamic forces. In such configurations, the proposed scheme is expected to work if the grid resolution is fine enough to resolve the wake. Indeed, the method only intends to resolve velocity and pressure gradi-ents at the fluid-solid interface. Future work will address this issue through the study of the sedimentation case of Cate [START_REF] Cate | Particle imaging velocimetry experiments and latticeboltzmann simulations on a single sphere settling under gravity[END_REF] and the Drafting-Kissing-Tumbling case described by Fortes [START_REF] Fortes | Nonlinear mechanics of fluidization of beds of spherical particles[END_REF].

Concluding remarks

The objective of this study was to clarify fluid-particle interactions for PRS with unresolved boundary layers. For this purpose, a new method to compute the hydrodynamic forces exerted by the fluid on the interface of a solid particle is detailed in this paper. Inspired of [START_REF] Chadil | Accurate estimate of drag forces using particle-resolved direct numerical simulations[END_REF][START_REF] Santarelli | Budget analysis of the turbulent kinetic energy for bubbly flow in a vertical channel[END_REF], the Lagrangian description of the interface to interpolate velocity and pressure on the Eulerian mesh is used. However, it does not require computing the full gradient velocity tensor. The originality of the present method consists in using the solid behavior of a particle that implies a non-deformation of the interface, resulting in a zero velocity gradient in the tangential components. This method thus saves time and computation memory.

A correlation was developed to correct the sub-resolution of hydrodynamic forces for coarse meshes. Applied in the volume of the particle, it corrects the action of the fluid to the particle for all resolutions thinner than 10 meshes per diameter. One may point out that only the amplitude of the force is corrected and not the direction, which is assumed to be correctly captured with these grid resolutions. Thereby, the velocity of the particle is corrected. The simulation with only 10 meshes per diameter is computed with the same precision on the particle transport velocity than the grid resolution of 40 meshes per diameter but with a CPU cost 430 times lower. A validation of the correction was realized for two Reynolds numbers and two density ratios within the Stokes flow regime. It exhibited a similar behavior for these four cases for four different grid resolutions. An assessment of the correlation beyond the Stokes flow regime enlightened that for a Reynolds number of unity, the density ratio influences the simulations. Moreover, as it was originally developed for Stokes flow regime, the error committed on the hydrodynamic force increases with the Reynolds number. This emphasizes the importance of considering both the density ratio and the Reynolds number in the correlation.

Figure 1 :

 1 Figure 1: Diagram of resolution scales for fluid-particle flows

Figure 3 :

 3 Figure 3: Diagram of the identification of neighboring Eulerian meshes for the hydrodynamic force computation. • : pressure nodes, •: Lagrangian marker, × ×: interpolation points.

Figure 4 :

 4 Figure 4: 2D diagram of stokes configuration. The size of the domain in z-direction is 3D p .

Figure 6 :

 6 Figure 6: Calculated streamwise velocity field, normalized by the sedimentation velocity, in a plane crossing the gravity center of the particle.

Figure 7 :

 7 Figure 7: Calculated pressure field, normalized by the maximum pressure, in a plane crossing the gravity center of the particle.

10 Table 4 :

 104 Physical parameters of Stokes configuration.

Figure 8 :

 8 Figure 8: 2D diagram of stokes configuration. The size of the domain in z-direction is 3D p .

Figure 10 :

 10 Figure 10: Calculated streamwise velocity field, normalized by the sedimentation velocity, in a plane crossing the gravity center of the particle.

Figure 11 :

 11 Figure 11: Calculated pressure field, normalized by the maximum pressure, in a plane crossing the gravity center of the particle.

Figure 12 :

 12 Figure 12: Local surface pressure scaled by the total theoretical surface pressure. Theoretical (left), computed from interpolation on a theoretical discretized pressure field (center), relative error (right).

Figure 14 :

 14 Figure 14: Mesh convergence of the overall hydrodynamic force relative error.

Figure 15 :

 15 Figure 15: Schematic representation of the staggered grid with the location of the variables: •: pressure nodes, ▲ ▶: velocity nodes, •: Lagrangian marker, : control volume V k of a cell face k.

  phase indicator of the face k and V k is the control volume of the face k.

Figure 16 :Figure 17 :

 1617 Figure 16: Construction of a mesh-dependent correction ×: computed with Eq. 37.-Interpolated function obtained with a least squares method.

Figure 18 :

 18 Figure 18: Overall hydrodynamic force, normalized by its theoretical value, as a function of the mesh resolution.

= 10 - 3 ,

 103 St = 1.1 × 10 -3 Re = 10 -3 , St = 0.33 Re = 10 -2 , St = 1.1 × 10 -2 Re = 10 -2 , St = 3.3

Figure 19 :

 19 Figure 19: Normalized hydrodynamic force as a function of the mesh resolution.

  = 2 ρp/ρf = 5 ρp/ρf = 10 ρp/ρf = 50 ρp/ρf = 100 ρp/ρf = 300 ρp/ρf = 500 ρp/ρf = 1000 ρp/ρf = 3000 Stokes PRS PR -SCS

Figure 20 :

 20 Figure 20: Robustness map of the proposed correction: representation of the hydrodynamic coefficient C h , normalized by Abraham's drag coefficient C dAbraham , with respect to the Re number. The black line shows the Stokes evolution in its range of validity, while the dashed line extend it for higher Re numbers. The simulation points are materialized by markers: •: PR-SCS, ×: PRS. The different density ratio ρ p /ρ f simulated are indicated by the color scale.

  •

Acknowledgments

This work was granted access to the HPC resources of CINES under the allocation 2022-A0122B11441 made by GENCI. The authors gratefully acknowledge the CEA for the development of the TRUST platform. The technical support of TRUST/TrioCFD team was also greatly appreciated. The authors acknowledge Salim Hamidi for his developments on the fluid-solid particles computation method.

Data Availability Statement

The data presented in this study are available on reasonable request from the corresponding author.

Credit authorship contribution statement


We present a new methodology for computing error bars in THz-TDS experiments and distribute open source software for this purpose. Pursuing our efforts in noise analysis and data processing of THz-TDS experiments, we followed a standard approach to derive the error bars from the simplest quantities coming from the experiments (transmission, absorbance) and a Bayesian approach when using an optimization fitting methodology to retrieve the plotted magnitude (refractive indices).

I. INTRODUCTION

IME domain spectroscopy (THz-TDS) is now in its third decade and grew from the bulky and very expensive setups of beginning using the rare femtosecond lasers of this time to very compact commercial and financially accessible apparatus. These progresses followed the ones of the femto-seconds lasers now available from optical fiber technology and as big as an external drive. The logical consequence of these developments is a spreading of the technics in many fields always further away from the initial community. Still today, the methods is not recognized as an analytical standard. We believe, with others [START_REF] Tayvah | Nelly: A user-friendly and open-source implementation of tree-based complex refractive index analysis for terahertz spectroscopy[END_REF], that one of the reason that THz-TDS did not take an even more important role in analytical field is the difficulty to derive uncertainties and error bars from the desired measured quantities.

In this conference, we will present how our efforts to fit the time trace of THz-TDS, and to analyze the noise of the systems [START_REF] Muller | Beyond the "dynamic range" approach in noise evaluation for terahertz time domain spectrometers[END_REF][START_REF] Lavancier | Heuristic approach to take up the challenge of terahertz time-domain spectroscopy for biology[END_REF] brought us further away and enabled us to derive error bars in transmission and absorbance measurements first, and even while retrieving the complex refractive index of a sample.

II. RESULTS

To test our approach we performed experiments on the socalled golden sample, a pellet of pure alpha lactose monohydrate. This sample have been widely studied and is an excellent candidate for our objective since it includes frequency region of relatively low absorption while other contains strong absorption peaks impeding any signal to cross the pellet in their frequency region.

To extract uncertainties and errors bars from quantities of relevance for our sample we have first to record the uncertainties on our signal itself. As we already [START_REF] Denakpo | Noise analysis, noise reduction and covariance estimation for time domain spectroscopy[END_REF] explained elsewhere, and following the methodology proposed by [START_REF] Mohtashemi | Maximum-likelihood parameter estimation in terahertz time-domain spectroscopy[END_REF][START_REF] Steven Dodge | A maximum-likelihood analysis framework for terahertz time-domain spectroscopy[END_REF] we first get rid of several artifacts such as drift in temperature of the optical fiber leading to a delay in the data and estimate the noise by measuring the standard deviation of our signal during the experiments including 1000 accumulations. This methodology allows the estimation of the noise first in the standard deviation in the time domain but also in the frequency domain.

Then using the frequency domain value for the uncertainties of the signal one can easily derive the one for a transmission measurements and with a little bit more of effort and assuming that the knowledge of the thickness of the sample does not add uncertainties one can derive the uncertainties on the absorbance as depicted on fig 1 .   Fig. 1. In black curve, the usual losses (being the minus of the log of the transmittance) plotted versus frequency, in blue the variance derived from the one on the signal for the sample and the reference. A nice point with this plot is that the lowering of the losses between 2.5 and 3 THz are confirmed to be true and not an artifact.

Here one can see that the uncertainty on the loss is very small in the region of strong signal going to the receiver. Hence the hypothesis of neglecting the uncertainties due to the measurements of the thickness of the sample maybe broken in this region. In the more absorbing region, one can see an increase of the uncertainties while keeping reasonable values, due to the log in the loss formula.

To go a step further we decided to extract the complex refractive index from our sample. This problem is fully invertible when dealing with non-absorbing or relatively low absorbing sample but one needs to perform an optimization as soon as the imaginary part of the refractive index become nonnegligible. Therefore getting uncertainties going through the optimization algorithm is one-step more complex. Here, we will use a Bayesian approach relating the uncertainties on the signal to the one on the retrieved n and k values by normalizing the error function in the fit routine by the noise. Hence, when the algorithm converges, the first derivative (gradient) is null and the variations are locally depicted by the second derivative (the Hessian). Using the normalization by the noise, the uncertainties on the retrieved value simply Noureddin Osseiran, Jeyan Bichon, Aditya Raj, Théo Hannotte, Sophie Eliet, Romain Peretti become the diagonal terms of the inverse Hessian. This is the methodology we used to plot the imaginary part of the complex index of refraction on fig 2. The variance goes from 1e-4 around 1 THz to 1e-2 around 3THz.

III. SUMMARY

In this conference, we will take the time to introduce how estimate the noise in THz-TDS experiments and propagate the uncertainties it introduces to the transmission, the transmittance, the refractive index and we hope that we will have reach the points of deriving the uncertainties of fitted magnitudes such as the central frequency or the width of a spectral lines.

We envision that these efforts will help to spread the technique to more analytical fields and become a new standard for myriads of applications.

Fig. 2 .

 2 Fig. 2. Imaginary part of the complex refractive index in black, in blue the variance of it.

Univ. Lille, CNRS, Centrale Lille, Univ. Polytechnique Hauts-de-France, UMR 8520 -IEMN, F-59000 Lille, FranceFrom Noise analysis to Error Bars on Refractive IndexIn THz-TDS