Hanane Bouras

Yacine Haroun
email: yacine.haroun@ifpen.fr

Francisco Fortunato Bodziony

Régis Philippe

Pascal Fongarland

Frédéric Augier

Use of CFD for Pressure drop, liquid saturation and wetting predictions in trickle bed reactors for different catalyst particle shapes

Keywords: Trickle bed reactors, pressure drop, liquid hold up, wetting efficiency, catalyst shape, CFD

The characterization of hydrodynamics in Trickle-Bed-Reactors is a complex task due to the opacity of the medium. In particular, the determination of pressure drop, liquid hold-up, wetting of the catalyst surface and catalyst shape effect on these parameters is very important for optimal catalyst use and reactor operation. Measurements under industrial conditions are limited to indirect estimations, and direct measurement methods are limited to near-ambient conditions. In this context, the objective of the present article is to use Computational-Fluid-Dynamics to investigate pressure drop, liquid saturation and wetting efficiency in Trickle-Bed-Reactors and to improve existing correlations, with a special focus on the catalyst shape effect and wetting prediction.

The Volume-Of-Fluid approach was used to simulate two-phase flow through particle loadings of spherical, trilobe and quadrilobe-shaped particles. The numerical model was

validated against literature correlations in terms of pressure drop, liquid holdup and wetting efficiency. The CFD model was then employed to explore two effects that does not reach out a consensus in existing literature, i.e effects of particle shape and gas-phase velocity on wetting efficiency. As a result, it was shown that CFD provides good predictions of pressure drop and liquid saturation for different catalyst particle shapes, the achieved deviations between CFD results and correlation estimations are below 20%. A new wetting efficiency correlation is also proposed. This new correlation is able to predict wetting efficiency with a precision of 6.99% for a wide range of liquid velocities (from 0.2 to 0.8cm/s) and gas velocities (from 5 to 20cm/s) and three particle shapes.

Introduction

Trickle-bed reactors (TBRs) are gas-liquid-solid continuous reactors, where gas and liquid flow co-currently downward through a fixed bed of solid catalyst particles. Parallel to fluids' flow, the reagents travel through the phases to react at the surface or inside the porous catalyst particles. Trickle-bed reactors owe their popularity to the unique advantages they offer for large volume processing in the chemical and petroleum industries. Indeed, TBRs are characterised with a good trade-off between (i) scalability, (ii) gas-liquid-solid interfacial contact (iii) mass, heat and momentum transports. Due to their relatively easy conception and operability, TBRs are the best technology for many industrial processes to perform hydrogenations (e.g. hydrocracking, hydrodesulfurization, hydrometallization, etc.) and are used as well in biochemical processes (e.g. Fischer-Tropsch [START_REF] Nishizawa | Use of Trickle Bed Reactor for Fischer-Tropsch Reaction over Co-Mn/Oxidized Diamond Catalyst[END_REF])) and wastewater treatment applications.

Several experimental works were conducted in the past four decades to characterise the hydrodynamic behaviour of trickle-bed reactors, leading to well-established hydrodynamic correlations for a wide range of operating conditions. Regarding two-phase pressure drop, available correlations can be classified into three main categories. The first category consists of correlations based on [START_REF] Lockhart | Proposed Correlation of Data for Isothermal Two-Phase, Two-Component Flow in Piped[END_REF] who correlated Δ𝑃 𝐺𝐿 to the following dimensionless numbers: (i) 𝜒 2 = Δ𝑃 𝐿 /Δ𝑃 𝐺 , (ii) Φ 𝐺 2 = Δ𝑃 𝐺𝐿 /Δ𝑃 𝐺 and (iii) Φ 𝐿 2 = Δ𝑃 𝐺𝐿 /Δ𝑃 𝐿 . This model was the foundation to many other experimental works [START_REF] Larkins | Two-Phase Concurrent Flow in Packed Beds[END_REF], [START_REF] Sato | Pressure Loss and Liquid Holdup in Packed Bed Reactor with Cocurrent Gas-Liquid Downflow[END_REF] , [START_REF] Midoux | Flow Pattern, Pressure Loss and Liquid Holdup Data in Gas-Liquid Downflow Packed Beds with Foaming and Nonfoaming Hydrocarbons[END_REF], [START_REF] Rao | Hydrodynamics of Two-Phase Cocurrent Downflow through Packed Beds[END_REF], [START_REF] Tosun | A Study of Cocurrent Downflow of Nonfoaming Gas-Liquid Systems in a Packed Bed. 2. Pressure Drop: Seach for a Correlation[END_REF], [START_REF] Ratman | A Model for Pressure Drop in Gas-Liquid Cocurrent Downflow through Packed Beds[END_REF]).

Nevertheless, the aforementioned dimensionless numbers (𝜒, Φ 𝐺 , Φ 𝐿 ) require the estimation of single-phase pressure drops (Δ𝑃 𝐺 𝑎𝑛𝑑 Δ𝑃 𝐿 ) in two-phase flow conditions, which proved to be challenging. Therefore, to overcome this challenge, a second correlation group based on bed characteristics and operating conditions emerged. Indeed, Turpin and Huntington (1967) suggested an independent variable Z (𝑍 = 𝑅𝑒 𝐺 1.167 /𝑅𝑒 𝐿 0.767 ) discussed and used by other authors [START_REF] Specchia | Pressure Drop and Liquid Holdup for Two-Phase Concurrent Flow in Packed Beds[END_REF], [START_REF] Rao | Hydrodynamics of Two-Phase Cocurrent Downflow through Packed Beds[END_REF], [START_REF] Sai | Pressure Drop in Gas-Liquid Downflow through Packed Beds[END_REF]). In their expressions, these correlations share some common points such as using the void fraction to the power of 3 to account for bed porosity effect on pressure drop. Finally, due to the empiric nature of the aforementioned correlation categories, a phenomenological model was presented by [START_REF] Holub | A phenomenological model for pressure drop, liquid holdup, and flow regime transition in gas-liquid trickle flow[END_REF] for the low-interaction regime, assuming that the flow inside the reactor is analogous to a flow in a 1D slit. [START_REF] Al-Dahhan | Catalyst Wetting Efficiency in Trickle-Bed Reactors at High Pressure[END_REF] later expanded on this work, by performing experiments at high pressures and widening the application range of the correlation, capable of predicting both pressure drop and liquid holdup. In an effort to investigate operating pressure effect on two-phase pressure drop, [START_REF] Larachi | Experimental Study of a Trickle-Bed Reactor Operating at High Pressure: Two-Phase Pressure Drop and Liquid Saturation[END_REF] proposed a correlation to estimate in the trickle flow regime. The authors were among the first to study the effect of high operating pressure on pressure drop. More recently, [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF] proposed a new mechanistic model to predict the pressure drop and liquid saturation inside the TBR. This correlation accounts for the tortuosity of the liquid film and provides good predictions of the hydrodynamic behaviour inside TBRs. An effort was also made by [START_REF] Bansal | Flow Regime Transition in a Trickle Bed Reactor[END_REF] to account for bed geometry at low and high-interacting regimes, for Newtonian and non-Newtonian fluids.

Besides, another crucial feature in trickle-bed reactors is wetting efficiency as it determines the extent of catalyst utilization. Catalyst wetting is observed at two different scales: (i) external wetting at the particle scale and (ii) internal wetting at the pore scale. The former is the fractions of external catalyst surface covered by a liquid film, and the latter is the fraction of pores' surface covered by liquid. In trickle-bed reactors, complete internal wetting is usually assumed due to strong capillarity effects [START_REF] Colombo | Solid-liquid contacting effectiveness in trickle bed reactors[END_REF]. Thus, research works focused on external wetting efficiency, henceforward referred to as wetting efficiency. Measurement techniques of wetting efficiency are either based on direct or indirect methods. The most widely used direct method is the dynamic tracer technique, which consists of injecting a dye tracer with the liquid phase into a trickle-bed-reactor, the wetted surface is then pigmented and measured to obtain the wetting efficiency. Although this approach is effective, it is usually applied in lab-scale reactors at low temperatures and pressures, operating at reduced superficial liquid velocities in order to keep the same liquid hourly space velocity as encountered in industrial-scale TBRs [START_REF] Colombo | Solid-liquid contacting effectiveness in trickle bed reactors[END_REF][START_REF] Burghardt | Experimental Studies of Liquid-Solid Wetting Efficiency in Trickle-Bed Cocurrent Reactors[END_REF][START_REF] Al-Dahhan | Catalyst Wetting Efficiency in Trickle-Bed Reactors at High Pressure[END_REF][START_REF] Van Houwelingen | Particle wetting distribution in tricklebed reactors[END_REF]Baussaron et al. 2007b;Baussaron et al. 2007a;Julcour-Lebigue et al. 2009a). Concerning indirect methods, the chemical reaction method is the most used technique to measure implicitly the wetting efficiency. Unlike the direct methods, this method is not limited to low temperature and pressure conditions, it is applicable under industrial operating conditions [START_REF] Schwartz | A New Tracer Method for Determination of Liquid-Solid Contacting Efficiency in Trickle-Bed Reactors[END_REF][START_REF] Sicardi | Catalyst Areas Wetted by Flowing and Semistagnant Liquid in Trickle-Bed Reactors[END_REF][START_REF] Al-Dahhan | Catalyst Wetting Efficiency in Trickle-Bed Reactors at High Pressure[END_REF][START_REF] Burghardt | Hydrodynamics and Mass Transfer in a Three-Phase Fixed-Bed Reactor with Cocurrent Gas-Liquid Downflow[END_REF]. However, a reactor model is needed to obtain wetting efficiencies through Residence Time Distribution (RTD), leading to a higher uncertainty.

Despite the fact that commercial reactors have an external wetting efficiency close to unity (Julcour-Lebigue et al. 2009a), many scale-up and hydrodynamic models are based on laboratory or pilot scale. [START_REF] Colombo | Solid-liquid contacting effectiveness in trickle bed reactors[END_REF] carried experimental tests in order to evaluate both external (at the pellet surface) and internal (inside the pellet's pores) effective wetting, the latter being found close to total due to capillary forces, whatever the wetting at the catalyst surface. [START_REF] Mills | Evaluation of Liquid-Solid Contacting in Trickle-Bed Reactors by Tracer Methods[END_REF] reported that the increase of superficial liquid velocity increased wetting efficiency, achieving unity at high velocities. It was also the first time a bounded function (guaranteeing 0<f<1) was proposed, using several dimensionless numbers.

More correlations were developed, many among them give dispersed results when analysing the same operating conditions. [START_REF] Van Houwelingen | Particle wetting distribution in tricklebed reactors[END_REF] used the colorimetric method to study the effect that pre-wetting conditions had on particle wetting, which may be one of the reasons why the discrepancies in wetting efficiency values happen. In addition, many of the techniques employed didn't receive appropriate evaluations regarding their precision and thus dispersed results appear. (Julcour-Lebigue et al. 2009b) modelled wetting efficiency accounting for several relevant parameters, such as liquid viscosity, bed porosity, particle diameter and wall effects. However, the effect of gas velocity has not been investigated by the authors. Results showed that the effect of particle shape on wetting efficiency is almost negligible for the investigated shapes. Besides, an increase in bed porosity decreases external wetting efficiency, this is mainly due to reduction of particle-particle contact points.

In short, large discrepancies arise between some of the available correlations when applied in the same range of conditions, mainly due to the experimental conditions employed and to the use of complex techniques.

Further investigation is needed to predict parameters of interest and cover the TBR hydrodynamic grey areas, especially for large scale reactors under high pressure and temperature operating conditions. The purpose of the present work is to use a Computational Fluid Dynamics (CFD) approach to achieve this goal.

Research into modelling multiphase flows using CFD has gained pace in recent years, especially since the proof of promising predictions in complex conditions [START_REF] Gunjal | Hydrodynamics of Trickle-Bed Reactors: Experiments and CFD Modeling[END_REF][START_REF] Kuzeljevic | Computational Modeling of Trickle Bed Reactors[END_REF][START_REF] Manoharan | Structure-Resolved CFD Simulations of Different Catalytic Structures in a Packed Bed[END_REF][START_REF] Jejurkar | Maldistribution Effects in an Industrial-Scale Trickle Bed Reactor[END_REF] are some CFD typical studies' examples. There are two main ways to approach multiphase flow in trickle-bed reactors as presented in the review of [START_REF] Wang | Modelling and simulation of trickle-bed reactors using computational fluid dynamics: A state-of-the-art review[END_REF]: (i) Euler-Euler method and (ii) interface reconstruction method. The former method considers that the particle loading is represented by an effective porous medium. This approach results in lower calculation times and simulation numerical cost, at the extent of poor local description of particular non-ideal packing effects for instance. [START_REF] Solomenko | Liquid spreading in trickle-bed reactors: Experiments and numerical simulations using Eulerian-Eulerian twofluid approach[END_REF] used the Euler-Euler method to simulate the liquid spreading in TBRs, and achieved a good agreement with volumeaveraged experimental data in terms of liquid dispersion, even when the reactor is fed by a single-point liquid injection. This latter progress proved the relevance of CFD in predicting liquid distribution inside trickle-bed reactors. Later, [START_REF] Augier | Characterization and modelling of a maldistributed Trickle Bed Reactor[END_REF] coupled the hydrodynamic model of [START_REF] Solomenko | Liquid spreading in trickle-bed reactors: Experiments and numerical simulations using Eulerian-Eulerian twofluid approach[END_REF] to a chemical reaction in order to study the effect of liquid maldistribution on reaction performances. Although this approach is effective in predicting phase distributions inside trickle-bed reactors, it has some disadvantages. The Euler-Euler approach does not enable to simulate flow at the particle scale, meaning that empirical or semi-empirical closure laws are needed to account for fluid/solid and fluid/fluid interactions. In response to this, the second approach -based on interface reconstruction -is increasingly used to simulate multiphase flows. Several numerical methods exist but the socalled Volume-Of-Fluid (VOF) approach is the most presented and used, as it is the most popular due to its simplicity and its possible use in complex geometries.

The VOF method allows tracking of explicitly defined fluid-fluid and fluid-solid interfaces, since it simulates flow over the actual physical geometry. Even though VOF method requires fine meshes, CFD studies have started using VOF method since it describes explicitly interfaces between the solid-liquid and liquid-gas phases, giving a more accurate description of the hydrodynamic phenomena. [START_REF] Augier | Numerical approach to predict wetting and catalyst efficiencies inside trickle bed reactors[END_REF] used the Volume-Of-Fluid method to investigate wetting and catalyst efficiencies in TBRs, the numerical domain consisted of three stacked spherical particles. A qualitative comparison between CFD and experimental flow snapshots shows a good prediction of the liquid distribution on the spheres. In addition, the predicted wetting efficiencies show similar tendencies as for experimental values. [START_REF] Haroun | Mass transfer and liquid hold-up determination in structured packing by CFD[END_REF] used the VOF approach to improve the understanding of gas-liquid flow behaviour on structured packings, and its effect on mass transfer performances. The same approach was then used by [START_REF] Horgue | Experimental and numerical study of two-phase flows in arrays of cylinders[END_REF] to predict the multiphase distribution on arrays of cylinders.

A comparison between CFD phase distribution and experimental snapshots shows an accurate prediction of the flow topology. More recently, [START_REF] Deng | Computational Investigation of Liquid Holdup and Wetting Efficiency Inside Trickle Bed Reactors with Different Catalyst Particle Shapes[END_REF] used the VOF approach to perform predictive simulations within packed loadings. The authors studied the effect of four different particle shapes as well as liquid properties effect on hydrodynamic parameters. The evolution of liquid holdup and wetting efficiency reported by the authors is in good agreement with previous observations in the literature. More recently [START_REF] Bouras | CFD modeling of mass transfer in Gas-Liquid-Solid catalytic reactors[END_REF] used VOF approach to simulate gas-liquid flow in structured multiphase reactors. A good agreement has been reported between experimental data and CFD simulations.

However, the VOF method has both advantages and limitations, it requires fine-enough meshes to avoid numerical diffusion at the fluid/fluid and fluid/solid interfaces, meaning that the required computational time might be very high. In addition, handling capillary flows with VOF method is intricate as spurious velocities and fluxes emerge at the interfaces for high surface tension values. Despite all these limitations, the literature shows that VOF approach is reliable and efficient.

Nowadays, with the development of computational resources, reasonable computing times can be reached with more complex geometries. The VOF method can now be applied on larger volumes of trickle-bed reactor geometries, involving hundreds of solid particles, and generated using DEM simulation software as suggested by [START_REF] Boccardo | Validation of a novel open-source work-flow for the simulation of packed-bed reactors[END_REF]. This article aims to answer two main interrogations: first, whether the VOF approach would predict pressure drop and trickling flow in representative volume of TBRs with a good accuracy and second, whether it can complete information obtained with direct experimental method to develop more precise wetting efficiency correlations.

To answer these questions, the first step is to build a realistic model of a 3D TBR, where hydrodynamic phenomena are simulated using CFD and compared with existing correlations.

In this study, about 40 points were simulated, for different fluid systems, operating conditions and particle shapes. The simulations results have been compared to well validated models from the literature. Based on the simulation data and experimental data of Julcour-Lebigue et al. (2009a), a new correlation to predict an overall external wetting efficiency is proposed.

Numerical model

This work uses Computational Fluid Dynamics in order to predict two-phase flow within trickle-bed reactors. A Volume-Of-fluid approach [START_REF] Hirt | Volume Of Fluid (VOF) Method for the Dynamics of Free Boundaries[END_REF]) is selected to accurately describe fluid-fluid and fluid-solid interfaces. The simulations are conducted using a commercial software ANSYS Fluent 19.2. In order to model the flow of immiscible fluids, the continuity equation (Eq. ( 1)) and the momentum equation (Eq. ( 2)) are solved simultaneously.

It is assumed that the fluids are Newtonian and incompressible, without phase change or heat transfer. In addition, surface tension is assumed to be uniform at the fluids' interface.

𝜕𝜌 𝜕𝑡 + ∇. (ρu ⃗ ) = 0 (1) 𝜕 𝜕𝑡 (𝜌𝑢 ⃗ ) + ∇. (𝑢 ⃗ 𝑢 ⃗ ) = -∇𝑝 + ∇(𝜇(∇𝑢 ⃗ + ∇𝑢 ⃗ 𝑇 )) + 𝜌𝑔 + 𝑇 𝜎 ⃗⃗⃗⃗ (2)
where 𝑢 ⃗ is the velocity, p is the pressure, 𝜇 and 𝜌 are respectively the mixture viscosity and density, 𝑔 is the gravitational acceleration and 𝑇 𝜎 ⃗⃗⃗⃗ is the surface tension force. The mixture density and viscosity are calculated using the following equations (Equation ( 3)) and ( 4))

𝜌 = ∑ 𝜌 𝑖 𝛼 𝑖 𝑖 (3) 𝜇 = ∑ 𝜇 𝑖 𝛼 𝑖 𝑖 (4)
where 𝛼 𝑖 is the volume fraction of the i-th phase. Interfaces' tracking is enabled through resolution of the continuity equation of the volume fraction transport for one (or more)

phases. For the 𝑞 𝑡ℎ phase, volume fraction transport is expressed as:

1 𝜌 𝑞 [ 𝜕 𝜕𝑡 (𝜌 𝑞 𝛼 𝑞 ) + ∇. (𝜌 𝑞 𝛼 𝑞 𝑢 ⃗ ) = 𝑆 𝛼 𝑞 ] (5)
where 𝑆 𝛼 𝑞 is the external mass source entering phase . This equation is not solved for the primary phase, its volume fraction is computed using the constraint in Equation ( 6):

𝛼 𝐿 + 𝛼 𝐺 = 1 (6) Volume fraction 𝛼 𝐿 is unity when only the liquid phase is present in the computational cell and zero when only the gas phase is present in that cell. The in-between values refer to the gas-liquid interface.

To account for capillarity effect, the Continuum Surface Force method (CSF) is employed. This non-conservative method is implemented through a source term 𝑇 𝜎 ⃗⃗⃗⃗ in the momentum equation. Proposed by [START_REF] Brackbill | A Continuum Method for Modeling Surface Tension[END_REF], the method assumes that surface tension is a continuous volume force. To further describe capillarity effects, a contact angle is specified at the catalyst's wall. Combined with surface normal, the contact angle allows to specify the interface's curvature, necessary to adjust the surface tension force 𝑇 𝜎 ⃗⃗⃗⃗ . Given the complexity of contact angle measurements in porous particles, a low contact angle value is usually employed [START_REF] Augier | Numerical approach to predict wetting and catalyst efficiencies inside trickle bed reactors[END_REF]). This numerical choice is supported by experimental observations, since porous particles behave like a sponge and are commonly assumed to be fully wetted internally. For the sake of simplicity, a static contact angle is applied at the solid surface as it has been proved to give satisfying results by [START_REF] Augier | Numerical approach to predict wetting and catalyst efficiencies inside trickle bed reactors[END_REF]. In this study a contact angle of 1° is specified at the catalyst particle surface.

This surface normal one cell away from the wall is given by Equation ( 7).

𝑛 ⃗ = 𝑛 𝑤 ⃗⃗⃗⃗⃗ 𝑐𝑜𝑠𝜃 𝑤 + 𝑡 𝑤 ⃗⃗⃗⃗ 𝑠𝑖𝑛𝜃 𝑤 (7)
where 𝜃 𝑤 is the contact angle at the wall, and 𝑛 𝑤 ⃗⃗⃗⃗⃗ and 𝑡 𝑤 ⃗⃗⃗⃗ are the unit vectors normal and tangential to the wall respectively.

The simulations are conducted in transient conditions. Since trickle-bed reactors usually operate at low liquid Reynolds numbers, laminar flow conditions are considered with a pressure-based solver. The PISO scheme is employed to ensure pressure-velocity coupling, while PRESTO! Is selected for pressure interpolation. The momentum equations are discretised using Second Order Upwind method. The resolution of two-phase flow throughout trickle-beds is achieved through the Volume-Of-Fluid approach, using the geo-reconstruct scheme. Finally, to guarantee numerical stability and convergence, a variable time step configuration is employed in order to maintain CFL condition under 1.

In order to simulate hydrodynamic performance within the considered loadings, boundary conditions shown in Figure 1 are considered. Regarding the multiphase flow, the gas and liquid enter the domain from the top, where their corresponding inlet velocities are specified. To describe the flow near the solid particles, a no-slip boundary condition is applied at the catalyst's surface. Since the numerical domains were extracted from larger loadings, the lateral boundaries are set to symmetry-type conditions. Finally, at the bottom, a regular pressure outlet condition is applied.

Packed bed loadings were generated in the open-source software Blender, as explained in the work of [START_REF] Boccardo | Validation of a novel open-source work-flow for the simulation of packed-bed reactors[END_REF]. Considering a particle shape, the packing generation is grounded on the Bullet Physics Library (BPL) available in Blender. This advanced physics simulation library solves the Newton's second law for a number of particles N, and is able to detect collisions between solid particles as well as their final position. According to [START_REF] Boccardo | Validation of a novel open-source work-flow for the simulation of packed-bed reactors[END_REF], Blender is capable of handling complex particle shapes. Therefore, three particle shapes are considered in this work, namely spheres, trilobes and quadrilobes.

Spherical particles are loaded in a mono-dispersed configuration, meaning that all the spheres have equal diameters. In contrast, poly-dispersion is included in trilobes and quadrilobes loadings, with a mean particle length of 4.8mm and standard deviation of 1.8mm. The packings were generated within cylindrical containers of equal diameter (D=3.5cm) and height (H=3cm), their characteristics are summarized in Table 1. Spheres and trilobes loadings have similar void fractions, while the quadrilobes loading is characterised by a higher void fraction, as can be seen from Table 1. Since loadings were generated following the exact same procedure, high quadrilobe bed porosity is attributed to the particle's geometrical features.

In contrast with spheres and trilobes, quadrilobes exhibit larger concavities.

In order to decrease computational cost, cuboid Representative Elementary Volumes (REV) are extracted from larger loadings in cylindrical containers. In order to obtain a surface-fitted mesh, the REVs are meshed using SnappyHexMesh following the recommendations of [START_REF] Boccardo | Validation of a novel open-source work-flow for the simulation of packed-bed reactors[END_REF], as shown in Figure 2.

The meshes were generated following a set of constraint on mesh quality metrics to avoid skewed and/or low orthogonal quality mesh cells. In addition, the mesh was refined near the catalyst surfaces to resolve sufficiently the contact points between particle shapes.

Consequently, the contact points were represented by bridges to ensure all previous criteria are respected. Since VOF is sensitive to mesh quality, SnappyHexMesh was used to generate the different meshes. The artificial bridges between particles represent less than 1% of the total void volume.

Three different fluid systems were considered in this study: (i) heptane-nitrogen at standard conditions, (ii) a gas-liquid hydrocarbon mixture with properties commonly encountered in the middle distillate hydrocracking processes, and (iii) isohexadecane-nitrogen at standard conditions. The properties of the aforementioned systems are presented in Table 2.

It is important to note that simulations were carried out in transient regime to prevent numerical instabilities or poor stabilization of transient phenomena. However, the results were analyzed after reaching steady state. The latter state is achieved after 3 second of flow time within the domain.

Numerical resolution

Mesh convergence study

As explained previously, the VOF approach requires a good quality mesh to provide accurate results. Therefore, a mesh independence study was performed on a representative volume of spheres in order to find the optimum mesh size. Four mesh densities are tested: (i) 1.5 Million, (ii) 3.4 Million, (iii) 4.3 Million and (iv) 11.2 Million cells. In order to capture the geometrical details, particularly at the solid-solid contact points, these meshes present a similar refinement level near the catalyst wall. The mesh independence study is conducted for constant liquid superficial velocity 𝑉 𝑆,𝐿 and gas superficial velocity 𝑉 𝑆,𝐺 conditions of: 𝑉 𝑆,𝐿 = 0.5𝑐𝑚/𝑠 and 𝑉 𝑆,𝐺 = 5𝑐𝑚/𝑠.

Figure 3 shows the level of refinement in each mesh and Table 3 summarises their corresponding resolutions and wetting efficiencies. Following the increase of mesh resolution, the predicted wetting efficiency values increase until reaching a plateau. In addition, almost tripling mesh density (from 4.3M to 11.2M) increases wetting efficiency by only 1.37%. It is interesting to note that simulation time increases from two days (4.3M) to fourteen days (11.2M) at constant CPU usage (128 cores).

Figure 3 illustrates contour plots of liquid volume fraction within the numerical domain for the considered meshes. The plots show clearly that the gas-liquid and liquid-solid interfaces are more distinguished for higher mesh resolutions. This means that numerical diffusion decreases with mesh resolution. Indeed, the coarser mesh (1.7 Million) exhibits considerable numerical diffusion, leading to a non-realistic representation of the gas-liquid interface. In fact, as the grid gets sharper, numerical thicknesses decrease, resulting in sharp fluid-fluid and fluid-solid interfaces.

In addition, mesh sharpness ensures a good definition of the boundary layer near the particles. As a result, the real particle loading geometry is better fitted, particularly at the contact points between particles. The selected gas and liquid velocities lead to the most severe conditions that are very sensitive to mesh resolution. The lowest liquid velocity is the lowest liquid film thickness is. This case requires fine mesh resolution in the thin liquid film to describe properly hydrodynamic properties.

The 4.3 Million cell mesh gives similar results to those obtained with the 11.2 Million cell mesh, as shown in Table 3. In addition, both meshes allow a sharp description of the gasliquid inteface. However, when compared to the 4.3 Million cell mesh, the 11.2 Million cell mesh takes 6 to 7 times more time to reach convergence. Therefore, the 4.3 Million cell mesh is employed in the hereinafter simulations since it provides a trade off between results' precision and resonable calculation times.

Inlet Effect on gas-liquid flow distribution and stabilization within the calculation domain

Liquid inlet effect on hydrodynamics is studied in order to identify the flow establishment zone within calculation domains. Indeed, since the liquid feed consists of 5 point sources, liquid distribution might require a certain bed depth to stabilise. Also called "calming bed depth", this stabilisation distance is achieved when liquid reaches a steady distribution. In order to check the distance needed for flow stabilisation, wetting efficiency is monitored in different regions of the numerical domains for HDT flow in several 𝑉 𝑆,𝐿 and 𝑉 𝑆,𝐺 conditions.

The numerical domain is divided into four sub-domains: (i) from 0 to 25% of the length (L), (ii) from 25%L to 50%L, (iii) from 50%L to 75%L and (iv) from 75%L to 100%L, given that inlet is located at z=0 and outlet at z=L. Then, the wetted solid surface is calculated within each subdomain as the integral of liquid volume fraction with respect to the solid surface. This study was conducted for all particle shapes at different operating conditions. However, results are presented only for spherical shape since equivalent behaviours are found for all particle shapes.

Figure 4 presents wetting efficiency evolution obtained at different bed depths for spherical particles' loading. In the investigated 𝑉 𝑆,𝐿 and 𝑉 𝑆,𝐺 conditions, it can be observed that wetting efficiency reaches a plateau for bed depths beyond 50%L, indicating that the developed flow is reached. Indeed, wetting efficiency increases on average by 6% between (i) and (ii) subdomains, then by 1.49% from (ii) to (iii) sub-domains, and finally by 0.03% from (iii) to (iv) subdomains. Moreover, wetting efficiency follows a similar evolution regardless of gas and liquid superficial velocities. Due to complex phase interactions occurring near the inlet, and in order to achieve accurate and consistent results analysis, post-processing results presented hereinafter are computed within the developed flow domain only.

Results and discussions

Pressure drop predictions

A preliminary study was performed in order to validate pressure drop across the reactors. First, liquid downward flow was simulated in single-phase conditions for various Reynolds numbers at steady-state. These simulations were performed for loadings of spheres, trilobes and quadrilobes. Regarding boundary conditions, liquid superficial velocities ranging from 0.02mm/s to 30.75mm/s are investigated. The latter superficial velocities correspond to 0.05 < 𝑅𝑒 𝐿 ′ < 100, where 𝑅𝑒 𝐿 ′ is a modified Reynolds number defined as follows:

𝑅𝑒 𝐿 ′ = 𝑉 𝑆,𝐿 𝑑 𝑒 ′ 𝜌 𝐿 𝜇 𝐿 (1 -𝜀 𝐵 ) (8)
𝑑 𝑒 ′ is a modified equivalent diameter of the particles defined as the following:

𝑑 𝑒 ′ = Ψ. 𝑑 𝑒 = Ψ. 6𝑉 𝑃 𝑆 𝑃 (9)
Since Ergun law [START_REF] Ergun | Fluid flow through packed columns[END_REF]] was originally developed for spherical particles, a sphericity factor is required to apply Ergun law [START_REF] Ergun | Fluid flow through packed columns[END_REF]] on non-spherical particles [Trahan et al. 2014]. Indeed, the sphericity factor Ψ given by equation ( 10) is a measurement of how close a particle's shape is to a sphere.

Ψ = 𝜋 1 3( 6𝑉 𝑃 ) 2 3 𝑆 𝑃 ( 10
)
where 𝑉 𝑃 and 𝑆 𝑃 are respectively the particle's volume and surface area. This factor is described as the ratio of (i) the surface area of a sphere with the same volume as the given particle to (ii) the surface area of the non-spherical particle.

Pressure drop estimations in porous media are achieved through Ergun law [START_REF] Ergun | Fluid flow through packed columns[END_REF]]. It is possible as well to estimate pressure drop from CFD simulations. To do so, the first step is to define two horizontal plane cuts within the developed-flow domain, at different height locations (𝑧 1 and 𝑧 2 ). Total pressure is then area-averaged within each plane (𝑃 1 and 𝑃 2 ), and pressure drop is obtained as Δ𝑃 = 𝑃 1 -𝑃 2 𝑧 1 -𝑧 2

. A comparison between CFD results and Ergun law [START_REF] Ergun | Fluid flow through packed columns[END_REF]] estimations is presented in Figure 5. The parity plot highlights the accuracy of CFD in predicting single-phase pressure drop, since 95% of results are within ±20% and the majority is close to the identity line.

After checking the ability of CFD to handle pressure drop in single phase flow conditions, the two-phase flow behaviour through three different loadings was predicted using the Volume-Of-Fluid approach. Table 4 summarises the investigated gas and liquid superficial velocities, as well as particle shapes. Thus, for each particle shape, 9 different (𝑉 𝑆,𝐿 , 𝑉 𝑆,𝐺 ) couples were simulated for the hydrocarbon-hydrogen system. 2007) is presented in Figure 6.

The results show a very good agreement between CFD simulation and the model. Indeed, 80% of the plotted points lie within the 20% deviation envelope, with an average relative error of 12.78. Pressure drop obtained with quadrilobes is low in comparison to spheres and trilobes as shown in Figure 6. This behaviour is expected due to the high void fraction of the quadrilobe particles' loading. In addition, it is observed that the most significant deviation between correlation and simulation results is obtained at the highest 𝑉 𝑆,𝐺 values. The deviation could be attributed to a transition from laminar to turbulent flow, which is not accounted for in these CFD simulations.

Gas-liquid flow behaviour in trickle bed reactors

Liquid saturation

Similarly to two-phase pressure drop, liquid saturation is computed from the predicted liquid volume flowing through the particule loadings in established flow region. Figure 7 shows typical multiphase flow inside the investigated representative volumes; after reaching developed flow conditions. In the results presented hereinafter, the liquid saturation is calculated by the ratio of CFD liquid volume to bed void volume.

For both investigated fluid systems, i.e hydrotreatment and heptane-nitrogen systems, predicted liquid saturation values are compared to the values calculated by the correlation of [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF], as shown in Figure 8. Even though CFD results slightly over-estimate liquid saturation compared to [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF] model's predicted value, the simulation results

show an interesting agreement with an average relative deviation of 34%. It is interesting to note that the highest deviations are obtained for quadrilobes. It is also important to point out that the model of [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF] has not been validated on this type of shape, which can explain the observed higher deviations. Moreover, it would be interesting to confront these CFD predictions to experimental data, not yet available.

From these results, it is observed that the higher the gas superficial velocity is, the lower the liquid saturation becomes. Such effect is a result of increased shear force at the gas-liquid interface, subsequently decreasing liquid film thickness through the particle loading. In addition, at constant gas superficial velocity, increasing liquid velocity results in increasing liquid saturation since liquid films become thicker as expected. Furthermore, quadrilobe shaped particles have particularly deep concavities leading to liquid accumulation (cf. Figure 7). Therefore, higher liquid saturation values are obtained for quadrilobes. Besides, in comparison to spheres and trilobes, the void fraction was previously reported to be higher for the quadrilobe loading, promoting lower interstitial gas velocities within the particle bed, and thus higher liquid saturations.

Liquid film thickness and flow regime investigation

It is well known that the liquid film thickness has significant effect on mass transfer and thus is of high importance in trickle-bed reactors. Since CFD allows access to the local flow parameters, the CFD predicted average liquid film thicknesses are compared to the ones predicted by Nusselt's falling film theory estimations [START_REF] Nusselt | Die Oberflachenkondensation des Wasserdampfes[END_REF] assuming total catalyst wetting, as shown in Figure 9. Initially developed for liquid flow over inifinite vertical planes, Nusselt's theory describes the evolution of liquid film thickness with liquid superficial velocity.

This theory is assumed to be a reasonable reference candidate to describe liquid flow over complex surfaces, such as particle loadings. Commonly, in steady-state conditions, the average liquid film thickness is defined as the ratio of liquid volume to the wetted catalyst surface. To avoid inlet effects, liquid film thickness is computed within the developed flow domain.

Results plotted in Figure 9 highlight the effect of gas and liquid superficial velocities on liquid film thickness. Indeed, the results show that increasing liquid superficial velocities leads to thicker films, while the increase of gas flow rate results in the opposite effect. These conclusions are consistent with previous observations on liquid saturation. As the gas velocity increases, the gas-liquid shear interactions are enhanced, resulting in liquid entrainment and liquid film thickness reduction. Besides, for indentical superificial velocities, liquid film thicknesses are on average higher for quadrilobes than the rest of the shapes. As was discussed previously, quadrilobe's geometrical characteristics (concavities) as well as lower interstitial gas velocities (high bed porosity) tend to promote an increase in film thickness.

Furthermore, results show that CFD and Nusselt's model [START_REF] Nusselt | Die Oberflachenkondensation des Wasserdampfes[END_REF]) predict a predominant effect of liquid superficial velocity on liquid film thickness. Thus, the model can be used to obtain an order of magnitude of the mean liquid film thickness, since mean realtive deviations of 14.7%, 18.1% and 20.6% are achieved respectively for spheres, trilobes and quadrilobes. Nevertheless, the effect of 𝑉 𝑆,𝐺 is not captured by the model.

Wetting efficiency predictions

Wetting efficiency is defined as the ratio of predicted wetted solid area to total solid area. For the sake of consistency, wetting efficiencies are computed in the developed flow domain, beyond the calming depth. The results are compared to the correlation of Julcour-Lebigue et al. (2009a), given by the following equation: -1.55 ] (11) where 𝜀 𝐵 is the bed porosity, 𝐹𝑟 𝐿 and 𝑀𝑜 𝐿 are respectively the liquid Froude and Morton dimensionless numbers, which accurately describe both liquid velocity and physical properties effects.

𝑓 = 1 -exp[-1.986𝐹𝑟 𝐿 0.139 𝑀𝑜 𝐿 0.0195 𝜀 𝐵
Figure 10 shows a comparison between wetting efficiencies obtained by CFD and the correlation of Julcour-Lebigue et al. (2009a). Several observations arise from the analysis of this parity plot. For spheres and trilobes, predicted wetting efficiency values are in a good agreement with the correlation of Julcour-Lebigue et al. ( 2009a) at low gas superficial velocities. However, a significant deviation is observed for quadrilobe particles as the offset between CFD and correlation values is much higher than 20%. Although Julcour-Lebigue et al.

(2009a) developed a wetting efficiency correlation for different particle shapes, their investigation did not include quadrilobe shaped particles. In addition, the authors performed experiments in the absence of gas flow, thus equation (11) does not account for gas superficial velocity 𝑉 𝑆,𝐺 effect. For these reasons, high deviations between CFD and correlation are reported in Figure 10.

It is interesting to note that CFD predictions show a negative impact of gas velocity 𝑉 𝑆,𝐺 on wetting efficiency. This is attributed to liquid disruption resulting in low liquid adherence to the solid surface. Wetting efficiency decreases with gas superficial velocities has also been observed in other CFD simulations, more recently on the work of [START_REF] Deng | Computational Investigation of Liquid Holdup and Wetting Efficiency Inside Trickle Bed Reactors with Different Catalyst Particle Shapes[END_REF]. The authors argued that these decreases are mainly due to higher gas-liquid interactions. Similarly [START_REF] Herskowitz | Global rates of reaction in trickle-bed reactors: effects of gas and liquid flow rates[END_REF] observed a global negative impact of 𝑉 𝑆,𝐺 on the global rate of a catalytic reaction. The proposed explanation was the increase of dry catalyst area with gas velocity, thus reducing the apparent effectiveness of the catalyst. This has been confirmed by [START_REF] Burghardt | Experimental Studies of Liquid-Solid Wetting Efficiency in Trickle-Bed Cocurrent Reactors[END_REF], who employed a mathematical description of the dynamic tracer method in order to review literature results. The authors reached the same observations. In this context, a new correlation is developed on the basis of CFD simulations performed in this work and former experimental data of , Julcour-Lebigue et al. (2009a) to include the effect of catalyst shape and gas superficial velocity

Development of a new wetting efficiency correlation

In their work, Julcour-Lebigue et al. (2009a) correlated wetting efficiency to three parameters, namely bed porosity 𝜀 𝐵 , liquid Morton number 𝑀𝑜 𝐿 and Froude number 𝐹𝑟 𝐿 . In this section, an improvement of the aforementioned correlation is proposed in order to include gas velocity and particle shape effects. In addition to predictions of hydrotreatment system and heptane-nitrogen, isohexadecane-nitrogen system was investigated by CFD in the conditions presented in Table 4. In an effort to broaden physico-chemical properties, isohexadecanenitrogen system was selected for its high liquid viscosity and surface tension.

Analysis of simulation results unveiled important physical scales to account for in the new correlation: (i) gas superficial velocity 𝑉 𝑆,𝐺 , (ii) gas-phase properties (𝜌 𝐺 , 𝜇 𝐺 ) and (iii) particle characteristics (𝑑 𝑃 , 𝑎 𝑆 ). Therefore, the proposed general formulation is given by:

𝑓 = 1 -exp [-𝐴 × 𝐹𝑟 𝐿 0.139 𝑀𝑜 𝐿 0.0195 𝜀 𝐵 -1.55 𝑆 𝐶 0 (1 -𝐶 1 ∏ 𝑁 𝑖 𝐶 𝑖 𝑛 𝑖=2 )] (12)
where 𝐶 0 , 𝐶 𝑖 , … , 𝐶 𝑛 are fitted constants, 𝑁 𝑖 are different gas-phase dimensionless groups, 𝐴 is a multiplication factor and 𝑆 is a shape factor.

Gas-phase dimensionless groups 𝑁 𝑖 are chosen from the following dimensionless numbers:

𝑅𝑒 𝐺 , 𝐹𝑟 𝐺 , 𝐺𝑎 𝐺 . Mathematical optimisation is conducted to determine the best group combination and fit the new correlation. The Sum of Squared Errors (SSE) between correlation and CFD/experimental values is minimised using a Generalised Reduced Gradient (GRG) algorithm. Table 5 summarises the obtained standard deviations with different group combinations. The minimum standard deviation was reached in Case 2 and Cases 4-7.

Therefore, for the sake of simplicity, the gas Froude number is selected to include gas phase properties effect on wetting efficiency. It is important to point out that for liquid flow at 𝑉 𝑆,𝐺 = 0, 𝐹𝑟 𝐺 is equal to zero and the correlation resumes to the model of Julcour-Lebigue et al. (2009a). Even though gas phase properties effect on the gas Froude number was not investigated in the scope of this work, the developed CFD model could be employed in a future work to establish this effect.

First, particle shape effect is described using a constant shape factor 𝜙, where 𝜙 depends solely on the particle shape and is equal to unity for spherical particles. The following correlation is then obtained:

𝑓 = 1 -exp[-1.986Fr L 0.139 Mo L 0.0195 𝜖 𝐵 -1.55 𝜙(1 -0.437𝐹𝑟 𝐺 0.367 )] (13)
Table 6 shows the various constant shape factors for each shape. As one can see, the equation is similar to the original one developed by Julcour-Lebigue et al. (2009). The correlation presented in equation ( 13) has not been evaluated for particle shapes other than spheres, trilobes and quadrilobes. Therefore, it cannot be used to estimate the wetting efficiency for any particle shape.

In order to expand the correlation of Julcour-Lebigue et al. (2009a) to different particle shapes and sizes, a different shape factor S is introduced as shown in equation ( 14). First, the shape factor 𝑎 𝑆 𝑑 𝑃 /𝜀 𝐵 2 used by [START_REF] Mills | Evaluation of Liquid-Solid Contacting in Trickle-Bed Reactors by Tracer Methods[END_REF] was tested, but a good wetting efficiency estimation could not be reached. Nevertheless, a good fit was achieved using the shape factor proposed by [START_REF] Ellman | A New, Improved Liquid Hold-up Correlation for Trickle-Bed Reactors[END_REF] :

𝑆 = 𝑎 𝑆 𝑑 ℎ 1 -𝜀 𝐵 (14)
This factor accounts for particle's shape and size through the hydraulic diameter 𝑑 ℎ . Whereas bed characteristics are included through specific area of the bed 𝑎 𝑆 and its porosity 𝜀 𝐵 .

Considering the aforementioned shape factor, the wetting efficiency correlation is expressed as the following:

𝑓 𝑐𝑜𝑟𝑟 = 1 -exp [-0.649𝐹𝑟 𝐿 0.139 𝑀𝑜 𝐿 0.0195 𝜀 𝐵 -1.55 ( 𝑑 ℎ 𝑎 𝑆 1 -𝜀 𝐵 )

1.147

(1 -0.436𝐹𝑟 𝐺 0.371 )] (15) [START_REF] Ellman | A New, Improved Liquid Hold-up Correlation for Trickle-Bed Reactors[END_REF] observed that the shape factor proposed in equation ( 14) could be further simplified since:

𝑎 𝑆 𝑑 ℎ 1 -𝜀 𝐵 ∝ 𝑐𝑜𝑛𝑠𝑡 × 𝜀 𝐵 3 (1 -𝜀 𝐵 ) 2 (16)
Due to the proportionality observed in equation ( 16), the correlation expression given by equation ( 15) can be further simplified. Another correlation expression, accounting for both bed void fraction and particle shape, is presented in the following equation:

𝑓 = 1 -exp [-4.065𝐹𝑟 𝐿 0.
139 𝑀𝑜 𝐿 0.0195 𝜀 -1.55 ( 𝜀 𝐵

3

(1 -𝜀 𝐵 ) 2 ) 0.376

(1 -0.434𝐹𝑟 𝐺 0.376 )] (17)

The standard deviations for each correlation with different particle shape factors are shown in Table 7. It can be seen that very good agreement is achieved by either expression, since deviation is under 10%.

The same standard deviation value is reached using equation ( 15) and equation ( 17), owing to the proportionality observed between shape factors used in both correlation expressions.

Wetting efficiency is estimated using equation ( 15), and compared to experimental and CFD results for over 110 data points. Figure 11 Without any doubt, Table 7 and Figure 11 confirm the good agreement between simulated and experimental values, highlighting the relevance of adding gas velocity and particle shape effects to the wetting efficiency correlation. Validity ranges of the wetting correlation are summarised in Table 8, excluding water and ethanol.

Conclusion

In this work, Computational Fluid Dynamics were employed to investigate the influence of catalyst shape, fluids propreties and gas-liquid flow operating conditions on pressure drop, liquid saturation, liquid flow behaviour, wetting efficiency in trickle-bed reactors.

A numerical model grounded on the Volume-Of-Fluid approach was used. The multiphase CFD results were compared to validated and well-established correlations from literature for pressure drop, liquid saturation, liquid film thickness and wetting efficiency. A good agreement was achieved between CFD results and existing correlations ( [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF]) and (Julcour-Lebigue et al. 2009a)). Therefore, these results validated the numerical model and showed that CFD can be a powerful tool to predict complex flow patterns in trickle bed reactors. In addition, CFD allows access to local phenomena and spatial distributions of transport parameters within the loadings, such as liquid film thicknesses.

The CFD model is then employed to investigate effects of (i) gas and liquid velocity, (ii) gas and liquid properties and (iii) catalyst shape on flow behaviour and wetting efficiency. It is found that increasing gas velocity leads to the decrease of wetting efficiency, with stronger effects noticed at higher gas densities. This was already observed in previous studies (such as the study developed by Al-Dahhan and Dudukovic (1995)). Moreover, particle geometry effect on hydrodynamic parameters was found greater for irregular shapes (quadrilobes). Overall results analysis revealed the significant influence of gas superficial velocity and particle shape on hydrodynamic parameters, particularly wetting efficiency. Therefore, an expansion of the wetting efficiency correlation of Julcour-Lebigue et al. ( 2009a) is proposed. In order to account for gas velocity and particle shape effects, two additional dimensionless numbers are added into the initial correlation. The gas velocity effect is included through gas Froude number 𝐹𝑟 𝐺 , while particle shape effect is accounted for through two shape factor expressions, depending on available characteristics.

The new correlation, developed in this study, can predict wetting efficiency with good accuracy (less than 10% deviation) for CFD experiments and already published experimental results of Julcour-Lebigue et al. (2009a). The correlation is valid for a wide range of physicochemical conditions, at high and low pressures, involving highly wetting liquid-solid systems as hydrocarbon and alumina. Further investigations may be required to extend the models to much less wetting liquid systems like water and ethanol.

In addition, this work shows that CFD VOF approach is an effective tool for the prediction of pressure drop, wetting efficiency and trickling flow in representative volume of TBRs, with a good accuracy and allow to develop more precise wetting efficiency correlations. A such work methodology presented in this work could be used for the innovation and optimisation of catalyst shape particles.

  [START_REF] Boyer | Hydrodynamics of trickle bed reactors at high pressure: Two-phase flow model for pressure drop and liquid holdup, formulation and experimental validation[END_REF] investigated pressure drop and liquid saturation in two-phase flow in the trickling regime and proposed predictive models for these two physical parameters. The authors validated these models over several experimental conditions and for different catalyst shapes. A two-phase pressure drop parity plot between CFD results and the model ofBoyer et al. (

  presents a parity diagram of correlation estimations against CFD and experimental results. It is possible to observe from Figure 11 that all multiphase CFD results and almost all experimental results (𝑉 𝑆,𝐺 = 0) are well fitted with the new developed correlation, even for more irregular particle shapes. It should be noted that higher deviations are observed for liquids -such as water and ethanol-. A corrective term was developed by Julcour-Lebigue et al. (2009a) which reasonably accounts for the contact angle effect through a critical liquid phase Froude number.

FiguresFigure 1 :

 1 Figures

Figure 2 :

 2 Figure 2: The numerical domains (upper row) and slices of their respective meshes (a) spheres (1.0x1.0x1.5 cm 3 , 𝒅 𝒆 =2.00mm, 𝜺 𝑩 =37.2%) (b) trilobes (1.0x1.0x2.1 cm 3 , 𝒅 𝒆 =1.80mm, 𝜺 𝑩 = 𝟑𝟕. 𝟕%) and (c) quadrilobes (1.0x1.0x1.5 cm 3 , 𝒅 𝒆 =1.95mm, 𝜺 𝑩 =62.8%).

Figure 3 :

 3 Figure 3: Minimum and maximum mesh densities of the spherical loading: (a) 1.7Million, (b) 11.2Million cell meshes and Liquid fraction contour plot for (c) 1.7Million (d) 11.2Million cell meshes.

Figure 4 :

 4 Figure 4: Wetting efficiency evolution within spherical particles' packing for HDT flow in different gas and liquid superficial velocities conditions.

Figure 5 :

 5 Figure 5: Parity diagram of single-phase pressure drop comparison between Ergun and CFD. Sphericity factors are respectively 𝚿 = 𝟏, 𝚿 = 𝟎. 𝟖𝟏 and 𝚿 = 𝟎. 𝟒𝟐 for spheres, trilobes and quadrilobes. 20% deviation envelope.

Figure 6 :

 6 Figure 6: Parity diagram of two-phase pressure drop comparison between (Boyer et al. 2007) and CFD for 𝑽 𝑺,𝑳 = 𝟎. 𝟐 𝒄𝒎/𝒔 (small markers), 𝑽 𝑺,𝑳 = 𝟎. 𝟓 𝒄𝒎/𝒔 (medium markers), 𝑽 𝑺,𝑳 = 𝟎. 𝟖 𝒄𝒎/𝒔 (large markers) and 𝑽 𝑺,𝑮 = 𝟓, 𝟏𝟎, 𝟐𝟎 𝒄𝒎/𝒔. Sphericity factors are respectively 𝚿 = 𝟏, 𝚿 = 𝟎. 𝟖𝟏 and 𝚿 = 𝟎. 𝟒𝟐 for spheres, trilobes and quadrilobes. 20% deviation envelope.

Figure 7 :

 7 Figure 7: HDT flow through the different numerical domains: (a) Spheres, (b) Trilobes and (c) Quadrilobes at 𝑽 𝑺,𝑳 = 𝟎. 𝟓𝒄𝒎/𝒔 and 𝑽 𝑺,𝑮 = 𝟏𝟎 𝒄𝒎/𝒔. The liquid is represented by the blue color.

Figure 8 :

 8 Figure 8: Parity diagram of liquid saturation comparison between (Boyer et al. 2007) and CFD for 𝑽 𝑺,𝑳 = 𝟎. 𝟐 𝒄𝒎/𝒔 (small markers), 𝑽 𝑺,𝑳 = 𝟎. 𝟓 𝒄𝒎/𝒔 (medium markers), 𝑽 𝑺,𝑳 = 𝟎. 𝟖 𝒄𝒎/𝒔 (large markers) and 𝑽 𝑺,𝑮 = 𝟓, 𝟏𝟎, 𝟐𝟎 𝒄𝒎/𝒔. 20% deviation envelope.

Figure 11 :

 11 Figure 11: Wetting efficiency parity diagram. Comparison between the correlation derived in this work (Equation (15)), CFD results and former experimental results of Julcour-Lebigue et al. (2009a). Dotted lines represent the  10% error envelope.

Table 4 -Particle shapes, bed porosity, liquid and gas superficial velocities tested to validate the model.

 4

	741
	Particle shape	Bed void fraction 𝜀 𝐵	𝑉 𝑆,𝐿 [𝑐𝑚/𝑠]	𝑉 𝑆,𝐺 [𝑐𝑚/𝑠]
			0.2	5, 10,
	Sphere	0.372	0.5	5, 10,
			0.8	5, 10,
			0.2	5, 10,
	Trilobe	0.377	0.5	5, 10,
			0.8	5, 10,
			0.2	5, 10,
	Quadrilobe	0.628	0.5	5, 10,
			0.8	5, 10,
	742

Table 5 -Standard deviations obtained with the different dimensionless numbers group combinations.

 5

	743
	Case			Exponents 𝐶 𝑖 for i>1	Standard
		C1				deviation 𝜎
	Number		𝑅𝑒 𝐺	𝐹𝑟 𝐺	𝐺𝑎 𝐺	[%]
	1	-0.092	0.283	0	0	7.40%
	2	-0.437	0	0.367	0	4.98%
	3	0.074	0	0	0.141	8.72%
	4	-0.459	-0.006	0.362	0	4.98%
	5	-0.451	0	0.359	-0.003	4.98%
	6	-0.451	0.718	0	-0.362	4.98%
	7	-0.451	0.118	0.300	0.062	4.98%
	744

Table 6 -Shape factor values for each particle shape.

 6

	745
	Shape	Spheres	Trilobes	Quadrilobes
	𝝓	1.00	0.91	2.41
	746

Table 7 -Standard deviation obtained by the model of Julcour-Lebigue et al. (2009a) and the new developed correlation for three shape factor expressions. Julcour-Lebigue et al. (2009a) correlation Shape factor for the new correlation

 7

	Eq. (13)	Eq. (15)	Eq. (17)

Table 8 -Range of physical properties of gases and liquids, geometric properties of packings and operating conditions. Physical properties of gases and liquids used

 8 Liquid density, 𝜌 𝐿 680 𝑘𝑔/𝑚 3 ≤ 𝜌 𝐿 ≤ 830 𝑘𝑔/𝑚 3 Gas density, 𝜌 𝐺 1.249 𝑘𝑔/𝑚 3 ≤ 𝜌 𝐺 ≤ 15 𝑘𝑔/𝑚 3 Liquid viscosity, 𝜇 𝐿 0.11 × 10 -3 𝑃𝑎. 𝑠 ≤ 𝜇 𝐿 ≤ 3.75 × 10 -3 𝑃𝑎. 𝑠 Gas viscosity, 𝜇 𝐺 1.3 × 10 -5 𝑃𝑎. 𝑠 ≤ 𝜇 𝐺 ≤ 1.7 × 10 -5 𝑃𝑎. 𝑠 Surface tension, 𝜎 𝐿 3.44 × 10 -3 𝑁/𝑚 ≤ 𝜎 𝐿 ≤ 28.1 × 10 -3 𝑁/𝑚 Operating conditions Liquid mass flow flux, 𝐿 1.34 𝑘𝑔/𝑚 2 𝑠 ≤ 𝐿 ≤ 60.1𝑘𝑔/𝑚 2 𝑠 Gas mass flow flux, 𝐺 0 𝑘𝑔/𝑚 2 𝑠 ≤ 𝐺 ≤ 3 𝑘𝑔/𝑚 2 𝑠 Liquid superficial velocity, 𝑉 𝑆,𝐿 0.2 𝑐𝑚/𝑠 ≤ 𝑉 𝑆,𝐿 ≤ 0.8 𝑐𝑚/𝑠 Gas superficial velocity, 𝑉 𝑆,𝐺 5 𝑐𝑚/𝑠 ≤ 𝑉 𝑆,𝐺 ≤ 20 𝑐𝑚/𝑠

	Geometric properties of packings used
	Equivalent diameter, 𝑑 𝑝	1.44𝑚𝑚 ≤ 𝑑 𝑃 ≤ 7𝑚𝑚
	Bed porosity, ε B	0.367 ≤ 𝜀 𝐵 ≤ 0.611
	Particle shape	Spheres, trilobes and quadrilobes
	Packing specific surface area, 𝑎 𝑝	516𝑚 -1 ≤ 𝑎 𝑃 ≤ 2604𝑚 -1

Tables Table 1: Characteristics of particles and packings.

Particle shape

Particle