domain: chim sde

text:
Hanane Bouras

Yacine Haroun
email: yacine.haroun@ifpen.fr

Régis Philippe

Frédéric Augier

Pascal Fongarland

CFD modeling of mass transfer in Gas-Liquid-Solid catalytic reactors

Keywords: Gas-Liquid-Solid reactors, mass transfer, CFD, volume of fluid method, film model, catalytic reaction

This work investigates Gas-Liquid-Solid mass transfer coupled to heterogeneous catalytic reaction using Computational Fluid Dynamics (CFD). The numerical model is based on the Volumeof-Fluid (VOF) approach, coupled with a convection-diffusion equation for mass transfer resolution.

First, the numerical method was validated on a case of falling liquid film over a semi-infinite planar surface. Then, a micro-structured reactor with α-methylstyrene hydrogenation is studied. A good agreement is found between experimental data of (Tourvieille et al. 2013) and simulation results.

Afterwards, a vertical spherical beads string is investigated. Convective transport by transversal velocities is identified as an important contributor to the overall Gas-Liquid-Solid mass transfer.

While the film model is applicable in pure diffusion regimes, the resistances-in-series model is not relevant and over-estimates the real mass transfer by nearly 30% when mass transfer occurs in liquid film flow without bulk.

The present work shows how CFD can be an effective tool for predicting hydrodynamic and catalyst geometry effect on mass transfer in Gas-Liquid-Solid reactors.

Introduction

Gas-liquid-solid reactors with a fixed bed are widely used in industry, they are encountered in several applications such as petroleum refining, chemicals production, pharmaceuticals, waste-water treatment, and many other fields. In most cases, the gas flows next to the liquid flowing around and wetting totally or partially a solid phase (catalyst), these reactors operate in co-current down-flow or up-flow configurations, as well as counter-current flow configuration. Regarding the widely encountered case where the limiting reagent is initially contained in the gas phase, due to the solute concentration gradients between phases, two interfacial gas-liquid and liquid-solid mass transfers resistance could take place before allowing a reaction on the catalyst particles. These external transfers can be merged as an overall gas-liquid-solid transfer coefficient. Reliable design of such processes depends on the ability to accurately predict the hydrodynamics, the interfacial mass transfer rates, and their couplings. However, the interfacial mass transfers are affected by hydrodynamic parameter variations such as film thickness and velocity field distribution. In addition, for the sake of simplicity, gas-liquid-solid mass transfer coefficients in such reactors are usually estimated using (i) the resistances-in-series model or (ii) the film model. On one hand, the resistances-in-series model was derived assuming a bulk concentration in the liquid phase. However, in the specific case of co-current down-flow occurring in trickle bed or in falling film reactors, as the liquid film thickness happens to be thin, the existence of a bulk concentration in the liquid phase is far from being obvious and the applicability of this widely encountered approach needs to be reviewed. On the other hand, the film model developed by [START_REF] Lewis | Principles of Gas Absorption[END_REF] defines mass transfer coefficients as the ratio of molecular diffusion coefficient to concentration boundary layer thickness. However, the difficulty lies in the right estimation of mass transfer boundary layer thickness in complex flows.

In recent years, Computational Fluid Dynamics (CFD) has proved to be a powerful tool for small scales and local phenomena analysis which are difficult to investigate experimentally. CFD enable access to velocity, pressure and concentration fields in complex systems and geometries through the resolution of continuity, momentum and concentration transport equations respectively. Several numerical works investigated the multiphase flow along solid surfaces, including mass transfer restrained to gas-liquid interfacial transfer. [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film[END_REF]) studied interfacial reactive mass transfer in a 2D stable falling liquid film, using the Volume-Of-Fluid approach (VOF) to solve the film hydrodynamics, and a modified concentration transport equation including the thermodynamic equilibrium at the gas-liquid interface. The simulated Sherwood numbers were compared to the analytical solution suggested by [START_REF] Higbie | The rate of absorption of a pure gas into still liquid during short periods of exposure[END_REF]) and a good agreement was achieved. [START_REF] Xu | Mass transfer across the falling film: Simulations and experiments[END_REF] used the VOF method to study the gas-liquid interfacial mass transfer in a 2D wavy falling film, in order to improve the understanding of the interface velocity fluctuations effect on gas-liquid mass transfer, the results were compared to an empirical experimental expression. (Haroun et al. 2010a;[START_REF] Haroun | Mass transfer and liquid hold-up determination in structured packing by CFD[END_REF][START_REF] Haroun | Use of Computational Fluid Dynamics for Absorption Packed Column Design[END_REF] applied the numerical model developed in [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film[END_REF]) on a 2D structured packing element, to study the effect of geometry on liquid film flow and mass transfer, a good agreement is found with [START_REF] Higbie | The rate of absorption of a pure gas into still liquid during short periods of exposure[END_REF] solution. [START_REF] Marschall | Numerical simulation of species transfer across fluid interfaces in free-surface flows using OpenFOAM[END_REF] investigated the mass transfer for 3D free-surface flows, namely a stagnant liquid film and gas bubble, using the VOF method and a modified concentration transport equation including Henry's law, the comparison with analytical and experimental results shows a good agreement. Later on, particles. Because of coarse meshing resolution, the authors found better predictions using Euler-Euler as opposed to VOF model. Later on, [START_REF] Jejurkar | Maldistribution Effects in an Industrial-Scale Trickle Bed Reactor[END_REF]) investigated the gas-liquid-solid mass transfer in a trickle-bed reactor, the two-phase flow is modelled using an Euler-Euler model where the gas-liquid and liquid-solid interfaces are not tracked explicitly. In summary, there are no contributions on numerical investigation of gas-liquid-solid mass transfer coupled with heterogeneous catalytic reaction using the VOF approach.

The aim of the present work is to use computational Fluid Dynamics to improve the understanding of mass transfer in gas-liquid-solid contactors, when the transfer is accompanied by catalytic heterogeneous reactions on the solid phase. The multiphase hydrodynamic resolution is based on the Volume-Of-Fluid (VOF) numerical approach, while mass transfer resolution is investigated by a convection-diffusion transport equation.

First, the CFD modeling is validated on a 2D laminar falling film with mass transfer along a vertical semi-infinite catalyst plane by comparing the results to available models from the literature. The model is then applied on a 3D falling film micro-structured reactor and pellet string reactor to couple hydrodynamics, mass transfers and heterogeneous reaction. The modelling results are compared to the experimental work of [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF], for α-methylstyrene catalytic hydrogenation in a micro-structured reactor. For all the aforementioned cases, a particular attention is attributed to the overall mass transfer calculation. The gas-liquid-solid overall mass transfer coefficient obtained by CFD is compared to (i) the gas-liquid and liquid-solid resistances-in-series modeling and (ii) the onefilm mass transfer modeling. The results are further discussed and an adapted modeling methods suggested.

Simulations of two phase flow are carried out with the volume of fluid method (VOF) with the commercial CFD software ANSYS Fluent 19. In the following section, the governing equations are described. Section 3 presents the model validation on a 2D falling liquid film over a vertical plane.

The results for a falling film micro-structured reactor and pellet string reactor are presented in Section 4.

Computational model

hydrodynamic model

The present work aims to couple multiphase flow with gas-liquid-solid mass transfer and reaction at the solid surface. To do so, the multiphase flow is solved using the volume of fluid approach (VOF), which is a numerical method suitable to describe immiscible fluid flows by resolving the transport equation ( 1 ) of the volume fraction in the computational domain.

⃗⃗

( 1 )

Where represents the liquid volume fraction ranging from 0 to 1. A value of 1 stands for liquid filled cells, whereas a value of 0 stands for gas filled cells. The interface is represented by the values in between. The VOF method tracks and describes explicitly the interface between the flowing fluids using interface reconstruction or finite volume discretization schemes.

The VOF approach might face limitations for complex flow configurations, such as high aspect ratio geometries, where spurious fluxes appear at the interface which becomes unstable. For this reason, the Coupled Level-Set (LS) and VOF (CLSVOF) method is used [START_REF] Sussman | A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows[END_REF].

To begin with, the Level-Set method is a numerical scheme which deals with fluid-interface motion, usually used for an accurate resolution of the interface topology. The LS method solves the evolution of a level-set indicator function as the following:

⃗⃗ ( 2 )
Where is positive on one side of the interface, negative on the other side, and equal to zero at the interface. The LS formulation captures the interface with a higher accuracy than VOF but suffers from mass conservation issues. However, these issues are overcome by combining the LS and VOF methods (CLSVOF) [START_REF] Sussman | A Coupled Level Set and Volume-of-Fluid Method for Computing 3D and Axisymmetric Incompressible Two-Phase Flows[END_REF]. In the scope of this work, CLSVOF is used when convergence issues or interface stability issues are encountered.

In addition, the conservation of mass (eq. ( 3)) and momentum equations (eq. ( 4)) are also solved for the gas-liquid incompressible and isothermal flow.

( ⃗⃗)

( 3 )

( ⃗⃗ ⃗⃗ ⃗⃗) ⃗⃗ ⃗ ⃗ ( 4 )
The force ⃗ considered in equation ( 4) accounts for external body forces other than gravity, pressure and viscosity. In this work, this force includes the effects of surface tension on the gasliquid and liquid-solid interfaces, it is modelled by the continuum surface force (CSF) model proposed by (Brackbill et al. 1992) as the following :

( ) ( 5 )
Where is the interface curvature, is the density of phase i, the volume-averaged density, and is the surface tension. The VOF model requires fine meshes for the gas-liquid interface to be sharp, coarse meshes may induce errors in the interface curvature values, which may lead to aberrant values of velocity, also known as spurious fluxes [START_REF] Meier | A novel technique for including surface tension in PLIC-VOF methods[END_REF].

To guarantee resolution accuracy of the volume fraction transport equation, time step size should fulfil the CFL condition (Courant-Friedrichs-Levy). This condition depends on the mesh and velocity profile. As the velocity profile changes at each iteration, the variable time-step option was used to ensure a maximum CFL value of 1.

The numerical simulations are performed using ANSYS Fluent v19.2, which is a commercial code.

The gas-liquid interface was advected using the Compressive Interface Capturing Scheme for Arbitrary Meshes (CICSAM), the mass and momentum conservation equations were discretized using the second order upwind numerical scheme, and the Coupled algorithm was used for the pressurevelocity coupling.

mass transfer model

For multiphase flows, ANSYS Fluent allows solving an additional scalar transport equation in the mixture, which is a User Defined Scalar (UDS) equation. To obtain the mass transfer coefficients between the different phases, the limiting reagent concentration was computed inside the domain using the latter approach as follows:

( )

( 6 )
Where u is the velocity, C is the solute concentration inside the domain, is the diffusion coefficient in (m²/s) and S is a user defined source term.

This work investigates the gas-liquid-solid mass transfer in the presence of a heterogeneous catalytic reaction at the surface of the solid catalyst using two methods. The first one consists of simulating a gas absorption into the liquid, coupled to its consumption by a heterogeneous reaction at the solid surface. The second one is based on the resistances-in-series model, where the gas-liquid and liquidsolid mass transfer coefficients are estimated separately, to be used to estimate the overall mass transfer coefficient.

Implementation of the heterogeneous catalytic reaction

In many cases, the heterogeneous reaction takes place in the porous catalyst layer volume, through a coupled reaction/diffusion process. It can be modelled by the addition of a source term through this catalyst volume, meaning that the diffusion inside the catalyst needs to be taken into account, as well as an accurate description of the porous catalytic layer. Implementing a catalyst discretisation in the CFD calculation would impact significantly computation time. Therefore, in this work, the heterogeneous catalytic reaction is accounted for by a surface reaction modulated with a surface effectiveness factor to take into account possible internal diffusion limitation. The following flux condition at the surface of the catalyst is derived.

(

) 7
Where is the intrinsic reaction rate, is the Pd mass fraction in the catalyst, and are respectively the catalyst volume and density, and is the surface efficiency factor which is estimated classically through the Thiele modulus . Since the flux condition (equation ( 7)) is applied at the catalyst's external surface, the concentration within the catalyst pores is not accessible. Smaller scale simulations are required to access such information. The -methylstyrene catalytic hydrogenation reaction was considered with a Pd/Al 2 O 3 catalyst. Due to its very fast intrinsic kinetics, this reaction is very often operated in full external mass transfer regime and is a good candidate to probe chemically overall external mass transfer efficiency of different GLS reactor configurations [START_REF] Meille | Kinetics of α-Methylstyrene Hydrogenation on Pd/Al2O3[END_REF][START_REF] Meille | Effect of Water on α-Methylstyrene Hydrogenation on Pd/Al2O3[END_REF]. The reaction rate simplifies to an n th order reaction in the selected experimental conditions:

( 8 )
Where and k is the reaction rate constant [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF]. As the reaction rate depends on hydrogen concentration, the concentration transport equation (eq. ( 6)) is solved for hydrogen only. The reaction is accounted for using the flux method explained above (eq. ( 7)).

Even though the flux condition is applied on the surface of the catalyst, the diffusion throughout the catalyst volume is accounted for by the efficiency factor, which is defined as the following for an nth order catalytic reaction and is estimated at each mesh cell:

√ ( ) ( 9 )
Where is the catalyst thickness, k is the reaction rate constant and is the effective diffusivity. This latter parameter is estimated as since the tortuosity ranges from 2 to 3, and the internal porosity is nearly 0.6 [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF]) and there is no Knudsen diffusion term.

Since the reaction rate and flux condition (eq. ( 7)) depend on the hydrogen concentration, only one concentration convection-diffusion equation is solved for concentration. The overall mass transfer coefficient can be post-processed after the solute concentration computation, the total hydrogen consumption flux is used to calculate the global mass transfer coefficient and the AMS conversion as follows directly linked to the 1:1 stoichiometry between the 2 reactants:

( ̅̅̅ ) ( 10 ) ( 11 )
is the gas-liquid interfacial area, is the liquid volume inside the reactor, is the AMS concentration fixed to 1 mol/m3 at the inlet, is the liquid volumetric flow rate. The gas-liquid interfacial area has been chosen arbitrarily.

is the thermodynamic equilibrium concentration of H 2 in the liquid phase (mol/m 3 ) given by the following expression [START_REF] Herskowitz | Solubility of hydrogen in.alpha.-methylstyrene[END_REF] :

( ( ) ) ( ) ( 12 )
Where T and P are the operating temperature and pressure.

̅̅̅ used in equation ( 10) is the mean hydrogen concentration at the catalyst's surface. Because of various possible definitions and for verification purpose, it was determined using two methods.

The first one consists in calculating the average concentration at the catalyst surface at the postprocessing stage of the CFD simulation. The second method is based on the hydrogen consumption flux at the catalyst surface. At the permanent regime, the hydrogen consumption flux is equivalent to the reactive flux at the catalyst surface, therefore the following model can be used to estimate ̅̅̅ :

̅̅̅̅̅̅ ( ̅̅̅ ) √ ( ) ̅̅̅ ̅̅̅ ( 13 )
Where ̅̅̅̅̅̅ is the averaged hydrogen consumption flux at the catalyst surface. Equation ( 13) is a non-linear equation which can be solved using an optimization algorithm such as a Levenberg-Marquardt one. In this work, the two approaches lead to very similar results, as the hydrogen surface concentration estimated by both methods are close, therefore only the first method is finally used in the following sections.

Resistances-in-series model analysis

Gas Liquid Solid reactor models usually use the resistances-in-series model to estimate the overall mass transfer coefficient, termed to distinguish it from presented in equation ( 10), considering that the overall external mass transfer inside such reactors is governed by two drivers in series, gas-liquid and liquid solid mass transfer resistances. First, a gaseous component diffuses from the gas bulk to the liquid film through the gas-liquid interface, the solute concentration stabilizes inside the liquid volume and reaches a liquid bulk concentration, then it is consumed by the heterogeneous reaction at the catalyst surface.

The overall mass transfer coefficient is determined using the resistances-in-series model, defined as follows:

( 14

)
This model is developed assuming the presence of a liquid bulk inside the liquid film. The gas phase consists of pure H 2 , and as the hydrogen is poorly soluble in the liquid phase, the gas-liquid mass transfer coefficient represents the liquid-side mass transfer.

The mass transfer coefficients and are determined respectively from simulations of:

 Physical gas absorption into the liquid through the gas/liquid interface  Physical solid dissolution into the liquid through the liquid/solid interface First of all, the physical gas absorption mass transfer coefficient is obtained after solving equation ( 6), considering a zero-concentration flux at the catalyst surface. The solute diffuses throughout the gas-liquid interface and accumulates inside the liquid film. At steady state, equation ( 6) takes the following form:

( ) ( 15 )
The convective flux difference ( ) between the inlet and the outlet of the domain is calculated to estimate the diffusion flux. The mass transfer coefficient is then defined from the diffusion flux as the following:

( ) ( 16 )
Where is the concentration at the interface, equal to the thermodynamic equilibrium concentration. Whereas is the mean solute concentration in the liquid film.

In order to determine the mass transfer coefficient , the solid dissolution is simulated considering the diffusion properties of Hydrogen from the solid to the liquid. The concentration convection-diffusion equation (eq. ( 6)) is solved considering a zero-concentration flux at the gasliquid interface. is calculated using the following equation:

( ) ( 17 )
Where is the concentration at the surface of the catalyst, and is the mean solute concentration in the liquid film.

The coefficients and are then used to estimate the overall mass transfer coefficients using the resistances-in-series model presented in equation ( 14).

Validation case

2D simulation of mass transfer in a liquid falling film

This section presents the simulations performed to develop and validate the numerical model. It consists of a 2D falling liquid film over a vertical plate, the aim is to guarantee a good agreement between the predictions and adequate analytical solutions, and to verify the validity of the resistances-in-series and film models in this numerical set-up.

Set-up & boundary conditions

The computational domain consists of a vertical plane plate of 1mmx250mm, water and hydrogen enter the domain at the top, and flow co-currently downward over the plate. The hexahedral mesh has a high cell density near the plane plate in order to improve the interface sharpness and concentration gradients in the liquid film

The boundary conditions are presented in Figure 1. The gas and liquid enter the domain at and respectively. The concentration is fixed to the equilibrium concentration at the gas inlet and gas-adjacent wall, and to at the liquid inlet. At the remaining boundaries, Neumann boundary conditions ( ) are specified.

Mesh convergence study

In their work, [START_REF] Haroun | Volume of fluid method for interfacial reactive mass transfer: Application to stable liquid film[END_REF] recommended to have at least 5 cells in the concentration boundary layer for an accurate concentration gradient description and gas-liquid mass transfer calculation. In order to determine the adequate mesh resolution to obtain accurate numerical description of hydrodynamics and overall gas-liquid-solid mass transfer, a mesh sensitivity study is performed. All the tested meshes are structured, for which the mesh resolutions are varied as shown in Figure 2. The criterion on which the final mesh is selected is mesh convergence calculation of the overall Sherwood number.

Figure 2 shows the evolution of gas-liquid-solid Sherwood number with mesh resolution. It can be observed that the overall Sherwood number does not change for meshes above 4.7 10 5 cells, corresponding to 14 cells within the liquid film thickness, indicating the mesh convergence of the simulation. It is interesting to note that the mesh used in this work are hexahedral with higher cell resolution near the catalyst surface in order to improve the interface sharpness and the concentration gradients calculation.

It is to be noted that for the other study cases : 3D Falling Film Micro-Reactor and 3D falling film on spherical particles were meshed with mesh resolution much higher than 14 mesh cells within the liquid film. Moreover, it is interesting to note that in the Falling Film Micro-Reactor, a very sharp mesh with a mesh resolution higher than the requirement of 14 cells is used in order to capture the film's meniscus shape.

Analytical solution

The liquid flow over such a semi-infinite vertical plane plate was studied by [START_REF] Nusselt | Die Oberflächenkondensation des Wasserdampfes[END_REF], who suggested an analytical expression for the interface velocity ( ) and liquid film thickness ( ) depending on the flow and physico-chemical properties. The analytical solutions are expressed as follows:

( 18 ) ( ) ( 19 )
Where is the liquid film thickness. This analytical solution was used to validate the numerical hydrodynamic results.

Concerning mass transfer validation, [START_REF] Higbie | The rate of absorption of a pure gas into still liquid during short periods of exposure[END_REF]) developed a model to obtain the local and global gas-liquid mass transfer coefficients in such films under laminar flow depending on the contact time, the coefficients are expressed as follows:

√ √ ( ) ( 20 )
The contact time is estimated based on the interface curvilinear length ( ) from the inlet to a position x.

Figure 3 shows the gas volume fraction distribution. It can be seen that the liquid film flows continuously over the plane plate, this film seems to have a constant thickness. The liquid film is slightly curved at the establishment zone near the inlet because the inflow velocity value is low.

The interface velocity and the liquid film thickness were compared to the analytical solution. As one can see from the Table 2, the numerical values are in good agreement with the analytical values with very low relative deviation for both evaluated parameters. When the flow is established, the simulated liquid film maintains a constant thickness throughout the domain, the liquid film thickness reported in Table 2 corresponds to the converged film thickness.

Mass transfer validation

Figure 4 shows a comparison between the analytical and numerical local mass transfer coefficients, for the 2D falling film over a vertical solid plate. The numerical mass transfer coefficient is obtained using equation ( 16), where the bulk concentration is taken at the furthest location from the gas-liquid interface, at the solid surface. As one can see, there is a high difference between the simulated and analytical values at the inlet of the domain, due to the difference between concentration boundary conditions at adjacent gas and liquid inlets, causing a high concentration gradient therein. However, overall, the numerical model reproduces the analytical solution with a good agreement in the rest of the domain.

The gas-liquid mass transfer coefficient is well predicted by CFD simulations as seen above.

Similarly, considering the dissolution of chemical species from a solid surface into the liquid film, and following the method described in the section 2.4, the liquid-solid mass transfer coefficient can be computed. The diffusive flux is determined by Fick's law, considering the diffusion coefficient of hydrogen.

The H 2 concentration profile in the liquid film when H 2 is consumed at the solid surface is represented in Figure 5. As one can see, the obtained profile is linear within the liquid film, which means that mass transfer takes place in the absence of a liquid bulk and due to the specific kinetics (zeroth order in substrate) and gas conditions (pure and isobaric H2 flow) this profile remains the same all along the falling length. The corresponding overall mass transfer coefficient obtained through equation( 10), for this reactive simulation is of 5.99.10 -5 m.s -1 .It can be compared to , its physical counterpart obtained following the resistance in series strategy (eq. ( 14)). In , the gas-reported in Table 3. This means that the diffusion layer thickness is nearly equal to the actual overall hydrodynamic liquid film thickness.

Study cases

3D simulation of a Falling Film Micro-Reactor (FFMR)

The developed numerical model has been applied to simulate one channel of a falling film microstructured reactor in order to validate its ability to predict hydrodynamic and mass transfer parameters in a more complex set-up, where the liquid film thickness takes different values depending on the location within the channel. This feature is obtained due to the structuration at the sub-millimetre-scale of the channel, thus involving a flow driven by capillary and viscous forces.

The purpose is to improve the understanding of mass transfer behaviour in similar film morphologies, to identify the most suitable mechanism to describe GLS mass transfer, and to compare the numerical results with available experimental results published by [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF]) on the hydrogenation of -methylstyrene with a coated Pd/Al 2 O 3 catalyst under mass transfer controlled regime but also on liquid film profiles established by fluorescent confocal microscopy under various flow conditions. These results were used to further validate the model both for hydrodynamic and mass transfer predictions.

The FFMR (Large version) is a falling film containing several parallel micro-structured channels, where the gas and liquid flow counter-currently, transferring the reagents from gas to the liquid and from the liquid to the solid, in the case of the presence of a coated heterogeneous catalyst. Microstructuration applied to falling film contact mode brings very stable G-L interfaces and intensification through liquid film thinning due to capillary forces. More details about the reactor can be found in [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF]) and in (Vankayala et al.).

Only one longitudinal half of a single channel has been simulated to reduce computation time, using a symmetry boundary condition in the middle of the channel.

As in the experimental work where two different plates were used, two numerical domains were considered: a first one without catalyst coating for the hydrodynamic validation, and a second one taking into account the catalyst layer thickness of 40µm.

As we assume that all the channels have an equivalent behaviour, the liquid is evenly distributed at the inlet, all the liquid channels have the same liquid and gas flow rates [START_REF] Tourvieille | Innovating microstructured gas-liquid-solid reactors: a contribution to the understanding of hydrodynamics and mass transfers[END_REF].

Moreover, liquid enters the domain at , and the concentration of hydrogen in the gas phase is constant for it consists of pure H 2 , and the interface H 2 concentration is equal to the equilibrium concentration given by equation ( 12 ).

To be fully representative of the experimental work, and because starting procedure plays a significant role on wetting, the channel is initially overflowed with the liquid at the inlet liquid velocity. The contact angle was not characterized experimentally, thus different contact angle values were tested in the model, and it is found that the liquid film profile is well represented with contact angle values of 20°.

Liquid flow over a flat plane can take several forms. At large liquid flow rates, a continuous liquid sheet is formed on the solid surface, which might breakup at certain locations. Falling film microreactors are known for preventing film breakup, they facilitate the stability of the liquid film and the corresponding gas-liquid interface and gas-liquid-solid 3-phase contact lines. Due to the combination of capillary forces and small channel widths, liquid is pulled up along the sides of the channels and covers up a significant portion of the channel width. The most desirable flow configuration is when the liquid flows in a meniscus shape as shown in Figure 6, that is to say for small Weber numbers, as the interface shape is governed by the interfacial energy.

Downward co-current two-phase flows were simulated, and converged liquid film profiles were compared to the corresponding experimental results. In order to avoid compatibility issues between the liquid phase and solid material of the FFMR, the hydrodynamic experiments were performed using ethanol/hydrogen instead of MCH+AMS/hydrogen. For consistency purposes, hydrodynamic simulations are performed using ethanol/hydrogen, and reactive mass transfer simulations are carried out with MCH+AMS/hydrogen. The properties of ethanol/hydrogen are given in Table 1.

Figure 6 shows this comparison for different inlet liquid flow rates. The numerical results are in good agreement with the experimental values indicating a very good prediction of the hydrodynamics by the CFD simulations.

The gas-liquid interface curves upwards moving from the channel's centre to the three-phase contact lines (always located at the edges of the channel due to the starting wetting procedure), and the interface is not disrupted across the numerical domain. The more the liquid flow rate increases, the flatter the interface becomes. Such behaviour was found in many experimental studies [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF][START_REF] Yeong | Characterisation of liquid film in a microstructured falling film reactor using laser scanning confocal microscopy[END_REF]. Besides, the liquid film is contained inside the channel regardless of tested liquid inlet flow rates in the investigated flow rate range. The results shown in Figure 7 were obtained with a contact angle value of 20°.

Logically, as shown in Figure 7, the corresponding specific gas-liquid area values are in good agreement with the experimental results, with a relative deviation to the experimental results only up to 10% at the lowest liquid flow rate. The variation of the specific gas-liquid area illustrates well the effect of liquid flow rate, as this latter increases, the interface flattens, the mean liquid thickness increases too leading to the decrease of the specific area.

Mass transfer

In this section, the two-phase flow of hydrogen and methylcyclohexane/ -methylstyrene mixture is investigated, and two channel dimensions were simulated for further validation. The coated channels characteristics are summarized in Table 4.

As one can see from Figure 8, diffuses from the gas phase to the liquid phase, and is consumed by the heterogeneous catalytic reaction at the solid wall. Due to the liquid film thickness variation inside the channel, the liquid is saturated with near the three phase contact line as the liquid film thickness is lower at this location. Whereas, at the channel's centre, the liquid film thickness reached its maximum value leading to a lower concentration at the catalyst surface. In addition, the mass transfer fully developed regime is reached at 5mm from the channel's inlet, thus the concentration profile is the same at different plane cuts of the channel as shown in Figure 8.

Experimental liquid volumes were estimated using a developed correlation in [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF], which can explain the 20% difference between numerical and experimental liquid volumes reported in Table 5.

As one can see from Table 5, the experimental and CFD predicted consumption fluxes are in the same orders of magnitude and follow the same evolution with liquid flow rate. The relative deviation between CFD and experimental results are 17%, 2% and 0.8% at liquid flow rates of 3ml/min, 5ml/min and 7ml/min respectively. The same findings apply on the -methylstrene conversions. However, the overall mass transfer coefficients show significant deviations between experimental values and numerical values, this is mainly due to experimental liquid volumes estimations. According to [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF], evaporation was noticed at low pressure in the experiment and a significant liquid volume was found at the gas outlet. Thus, the effective liquid volume inside the reactor in reaction conditions was lower than the one predicted through the correlation and the overall mass transfer coefficients were overestimated. However, evaporation was negligible at 5 bar, which means that data measured at 5 bar are the most reliable ones, especially concerning the H 2 consumption flux [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF].

The deviation in even for the best experiment remains large because of an important difference in the liquid volume determination as already mentioned. For the experimental results, it is approached through a correlation and for this work it is precisely simulated using the hydrodynamic model. That's why in the following, comparisons with experimental data will be carried out on the H 2 consumption flux (which is equivalent to a comparison in ). within 10% around the experimental ones, this allows us to conclude that CFD was able to capture and predict the effect of channel dimension on this coupled situation involving hydrodynamics, mass transfer and heterogeneous reaction. One interesting fact, both numerical and experimental mass transfer coefficients reach higher values with bigger channels which is counterintuitive in microstructuration. Indeed, micro-structuration is needed to obtain a curved and thinned liquid interface in comparison to unstructured conventional falling. films. However, because of the presence of the free interface, too much structuration can be detrimental in this case, leading to thicker films with the smaller channels at the same mean flow velocity or flowrate. This holds true as long as the 3phase contact lines remain stable and at the edges of the channels.

Figure 10 shows a comparison between experimental and predicted values for -methylstyrene conversion at P=5bar. As one can see, numerical predictions are in good agreement with experimental values thus validating the developed CFD model, the liquid flow rate effect on conversion is also well represented.

The resistances-in-series model has been tested on this case study as well. The results are summarized in Table 6. They show on average 32% deviation between the computed overall mass transfer coefficients, , and the resistances-in-series model, . These results show that the resistances-in-series model is not adapted to describe the overall mass transfer coefficient, the liquid film thicknesses obtained in FFMRs are low, leading to the absence of a bulk concentration within the liquid film.

Similarly to the semi-infinite planar case, the applicability of the film model was tested. An interesting difference lies in the definition of the correct characteristic length because of the encountered variable thickness in the cross section of it. Thus, different liquid film thicknesses were used to define the film model mass transfer coefficient, namely the liquid thickness at the centre of the channel , the arithmetic mean thickness , the geometric liquid thickness defined as the ratio of the liquid volume to the wetted surface, and finally the harmonic mean thickness .

The film model underestimates the mass transfer coefficient using , and strongly overestimates it using ; the lowest mean relative deviation reported in Table 7 for corresponds to the arithmetic mean liquid film thickness , the mean relative deviation is above 30% for the remaining liquid film thickness definitions. These differences are caused by the liquid film morphology, where the film thickness is very low near the three-phase contact line and reaches a maximum at the centre of the channel, causing important local differences in H 2 consumption at these locations. Nonetheless, to reconciliate the sophisticated CFD simulations with the simple film model, a correction factor A can be introduced to the film model to predict well the mass transfer:

The results are presented also in Table 7 and the best fit was obtained using . However, is considered to be the most interesting thickness value, for its easier acquisition both experimentally and numerically.

3D simulation of a reactive falling film on string of spherical catalyst particles

The previously studied cases consisted of falling films over a regular plane or micro-channel. In order to highlight the effect of the solid catalyst shape and to go towards more complex geometries encountered in fixed beds, simulations of reactive two-phase flow over stacked 1mm diameter spheres were carried out. The purpose of this case study is to investigate the effect of liquid film thickness distribution within the domain on GLS mass transfer involving a heterogeneous catalytic reaction, and to check the validity of resistances-in-series and film models on a such complex configuration. Contrary to the two previous cases, here the chosen geometry is likely to induce convection in other directions than the main flow direction and its impact on overall external transport from the gas phase to the solid catalyst surface will be investigated.

The three dimensional simulation domain consists of 12 spheres of 1mm diameter stacked inside a cylinder, where the co-current hydrogen/water flow takes place by gravity from inlet to outlet. The mesh was generated using SnappyHexMesh, with high refinement near the spheres to increase the number of cells in the liquid film.

For the first hydrodynamic investigation, the gas and liquid inlet Reynolds numbers are 25 and 30 respectively, and the flow takes place in isothermal (T=283K) and isobaric conditions. The hydrogen concentration at the liquid inlet is equal to zero, and at the gas inlet, it is homogeneous in the gas phase a it consists of pure hydrogen, and reaches the equilibrium concentration given by equation ( 12) at the gas-liquid interface.

The liquid forms a continuous film over the catalyst particle string at steady state as shown in is thin and vice-versa, due to mass conservation inside the computational domain. The total number of spheres was chosen in order to ensure the fully developed regime is reached.

Mass transfer analysis

The convection diffusion equation of concentration (eq. ( 6)) was solved at different pressure conditions, that is to say at different thermodynamic equilibrium concentrations . In order to investigate the effect of particle shape on mass transfer, purely diffusive conditions are also considered. The pure diffusion is an artificial calculation since the convective term is artificially considered equal to zero in the concentration transport equation, but the liquid film thickness profile remains the same as for convection-diffusion conditions.

As one can see from Figure 13, when convection is neglected in transport analysis, the diffusion boundary layer is equivalent to the whole liquid film thickness (Figure 13-a). Whereas when the convection effects are taken into account, the diffusion boundary layer is closer to the particles (Figure 13-b), the liquid film is loaded with hydrogen since the concentration is maximum in a large proportion of the liquid film, the concentration then decreases near the surface of the spheres where it is consumed by the heterogeneous catalytic reaction.

In order to bring out the effect of convection, Figure 14-a compares the concentration profiles in the liquid film at the great circle of a particle. As one can notice, the concentration profile is linear when convection is neglected. In contrast, the convection boosts the solute transfer since the concentration is nearly constant far from the spheres, and the diffusion is predominant near the spheres as the concentration profile becomes linear. convection. The flux is significantly increased by convection in comparison to pure diffusion, it is on average 6 times higher. This result is substantially different from the two previous study cases for the radial velocity profile contributes to radial transport enhancement in this case. In addition, one can also notice that the flux profile shape is smooth for pure diffusion and reaches peaks at the great circle of particles. In contrast, when convection is considered, the flux profile shape is asymmetric and reaches a peak at the top of each particle, following the shape of the radial velocity as shown in The same evolution is noticed for the distribution of the solute surface concentration on the spheres surface. Figure 16 shows the difference in surface concentration of the solute with and without convection on the 10 th sphere. For pure diffusion, the surface concentration is symmetric, whereas the radial velocity boosts the concentration at the first hemisphere of the particles when convective transport is considered.

The local overall external mass transfer obtained for convection-diffusion is represented in is calculated at each mesh cell of the gas-liquid interface, and defined as the ratio of the local hydrogen consumption flux and the local concentration difference ( ̅̅̅̅̅̅ ). follows the same trend as the surface hydrogen consumption flux and surface concentration, Figure 17 reports only the data at P=1bar because the profiles are exactly the same for P=3bar and P=5bar. After going through the flow stabilization zone at the inlet, for the 2-4 first spheres, the mass transfer coefficient profile becomes similar from one sphere to another and is maximized at the first hemisphere of each particle, with the liquid acceleration, then decreases at the second hemisphere with the liquid deceleration.

The average external overall mass transfer coefficients are obtained by integrating the local overall external mass transfer coefficient presented in Figure 17 throughout one sphere in the fully developed mass transfer regime, that is to say on the 10 th sphere for instance. As done previously, the apparent diffusion layer thickness is obtained from the average using the film model. The apparent diffusion layer thicknesses reported in Table 8. The results show that the radial convection enhances mass transfer. The apparent diffusion layer is low when convection is accounted for, this means that the diffusion takes place in a small liquid film thickness near the surface, which was explicitly noticed in Figure 14-a.

As shown in Table 8, The mass transfer coefficients determined by resistances-in-series model are on average 37% higher than the obtained overall mass transfer coefficient, once again because the gas-liquid and liquid-solid mass transfer layers overlap, since there is no bulk concentration in the liquid film.

Concerning the film model, neither the geometric thickness nor the average thickness gives good results, as the apparent diffusion layer thickness is very small compared to the liquid hydrodynamic thickness.

Further investigation of the convection effect on mass transfer has been made, considering several liquid Reynolds number conditions to investigate the liquid velocity and viscosity effects, the diffusion coefficient remains fixed at 4.88.10 -9 m²/s. As one can see from In order to gather all these data in a single correlation, a gas-liquid-solid Sherwood number corresponding to the overall external mass transfer has been estimated. Figure 19 shows that this Sherwood number increases with the Reynolds number, and seems to reach an asymptotic value of 6 in the investigated conditions since the flow regime is laminar. When the convection is neglected, the mass transfer coefficient was equal to , leading to a Sherwood number equal to 1. Thus, a tentative correlation has been derived for this overall Sherwood number, taking the following form:

( ) ( ) ( 25 )
Where is the Sherwood number for pure diffusion and is the asymptotic value of the Sherwood number equal to 6 in this case. The constant A was fitted on the obtained numerical data, and the tentative correlation is given by:

( ) ( 26 )
As one can see from Figure 19, Sherwood numbers estimated by the correlation presented in equation ( 26 ) and by CFD are close, the maximum relative error is 11% and the average relative error is 4.74%.

In summary, This case illustrates the effect of a more complex geometry corresponding to pellet string reactor on hydrodynamic and mass transfer performances. It is found that the convection enhances drastically the mass transfer rate, in contrast with the two previous study cases where diffusion was the exclusive radial transport phenomenon. However, since the string of few particles are stacked in a perfectly aligned pattern, the two-phase flow behavior is different from the one encountered in trickle-bed reactors, where the catalyst particles are stacked randomly and are more packed, leading to more complex phase distributions inside the reactor. Nevertheless, this work

shows that CFD simulation allows investigating realistic complex reactive systems such as trickle-bed reactors.

Conclusion

A numerical model has been developed in order to improve the understanding of reactive falling liquid films flow over different catalytic surfaces. This model coupled a robust hydrodynamic VOF description with a consistent gas-liquid-solid mass transfer and a surface heterogeneous catalytic reaction. Three case studies have been simulated, namely a 2D semi-infinite falling liquid film over a vertical planar surface, a 3D falling film micro-reactor and a 3D falling liquid film over a string of spherical particles.

First, the CFD modeling of two-phase flow in bidirectional domain and three dimensional microchannel reactor have been investigated. The simulation results were compared respectively to Nusselt's model [START_REF] Nusselt | Die Oberflächenkondensation des Wasserdampfes[END_REF] and experimental data of [START_REF] Tourvieille | Mass transfer characterisation of a microstructured falling film at pilot scale[END_REF]). The simulation results show a very good agreement in both cases.

Regarding mass transfer, the overall external gas-liquid-solid mass transfer behavior encountered under heterogeneous catalytic reaction conditions is significantly different in the three study cases. Indeed, For the 2D semi-infinite plate, the gas-liquid mass transfer coefficients is found to be in good agreement with the analytical solution of [START_REF] Higbie | The rate of absorption of a pure gas into still liquid during short periods of exposure[END_REF]. In addition, the gas-liquid-solid mass transfer is well described by the film model, when a pure radial diffusion regime is reached at steady state. Moreover, this work shows that the resistances-in-series model is not suitable to describe mass transfer for this case, since it overestimates the mass transfer coefficient by 27%.

Regarding the 3D falling film micro-reactor (FFMR), a good agreement is obtained between experimental and simulated overall mass transfer coefficients. It is found that the gas-liquid-solid mass transfer is predominantly diffusive, a correction factor is proposed and introduced to the film model to account for the non-uniform liquid film thickness in the channel. For FFMR system, it is also found that the resistances-in-series model overestimates the overall mass transfer coefficient by 32% for similar aforementioned reasons as the 2D falling liquid film.

Finally, using the validated CFD model, mass transfer in trickling flow conditions over a string of 12 spherical catalyst particles with a heterogeneous catalytic reaction at the catalyst surface has been studied, using the validated CFD model. It is found that the gas-liquid-solid mass transfer regime is mainly dominated by the convection induced radially by this geometry. The convection enhances mass transfer, leading to thin mass transfer boundary layers. In addition, it is found that the film model does not give representative results for this case, and a new Sherwood number correlation is developed to correct the model. This work shows also that the overall GLS external mass transfer coefficients estimated by the resistances-in-series model are on average overestimated by 37% again because of bulk concentration absence.

To conclude, this work shows that the mass transfer in two-phase laminar flows with reaction occurring at the solid surface are dependent on the solid geometry and behave differently regarding the overall external mass transfer. Not all the cases are equivalent to a simple falling liquid film, thus each case needs to be studied separately. In addition, this work proves that CFD can be an a powerful tool, not only to predict complex flow patterns, but also to simulate physically relevant mass transfer processes coupled with heterogeneous reactions, as well as to improve the understanding of multi-physics phenomena in gas-liquid-solid fixed bed reactors. Nowadays , thanks to the development of high performance calculations resources, this approach might be extended to fixed-bed reactors. However, the predictions would be limited to a few hundreds of particles instead of the entire reactor. Even though the simulated scales are far from the reactor scale, the presented numerical model can be applied to improve understanding of local mass transfer mechanisms.

Publication bibliography

(

  [START_REF] Sebastia-Saez | 3D modeling of hydrodynamics and physical mass transfer characteristics of liquid film flows in structured packing elements[END_REF]) evaluated the two-phase film flow and gas-liquid mass transfer along a structured packing element, the results show a reasonably good agreement with theoretical and experimental data, such as the liquid film thickness and Sherwood number.(Chii-Dong et al. 2011) conducted 2D and 3D simulations of the two-phase flow inside a falling film micro-reactor, the numerical model is based on the VOF approach, this preliminary study highlights the complexity to predict flow behavior in micrometric high aspect ratio geometries, since convergence difficulties were faced particularly for the interface stability. Despite the numerous scientific studies on twophase free-surface flow simulation and interfacial mass transfer, only few numerical studies simulated the gas-liquid-solid mass transfer considering a heterogeneous catalytic reaction. The majority of these works are based on the Eulerian framework, where the gas-liquid and liquid-solid interfaces are not tracked explicitly. For instance,Lopes et al. (2007) andLopes and Quinta-Ferreira (2007) proposed a CFD model coupling hydrodynamics to catalytic wet air oxidation (CWAO) of vanillic and phenolic acids respectively. The authors validated their CFD model against hydrodynamic experimental data, before conducting predictive simulation of the reactive performance. In the same context, Lopes and Quinta-Ferreira (2010) used both the Euler-Euler and VOF method to predict reactive performance of CWAO. The numerical domain consisted of 200 non-overlapping spherical

Figure 9

 9 Figure 9 compares experimental and numerical , as well as the evolution with liquid

Figure 11 -

 11 Figure 11-a, the thickness varies throughout the domain. As one can see from Figure 12-a, the liquid

Figure 14 -

 14 Figure 14-b shows a comparison between the hydrogen consumption flux with and without

  Figure 15.

  Figure 17.

  Figure 18-a, the geometric thickness increases with Reynolds number as well as viscosity. The overall mass transfer coefficient is represented in Figure 18-b since the liquid film thickness increases with the Reynolds number, the mean mass transfer layer becomes thicker and the mass transfer coefficient decreases. However, convection enhances mass transfer since the convective mass transfer coefficient remains higher than the purely diffusive mass transfer coefficient.

  Brackbill; J. U.; Kothe; D. B.; Zemach; C. (1992): A Continuum Method for Modeling Surface Tension. In Journal of Computational Physics 100, pp. 335-354. Chii-Dong, H.; Chang, H.; Chen, H. J.; Chang, C. L.; Li, H. H.; Chang, Y. Y. (2011): CFD simulation of the two-phase flow for a falling film microreactor. In International Journal of Heat and Mass Transfer 54 (15-16), pp. 3740-3748. DOI: 10.1016/j.ijheatmasstransfer.2011.03.015.

FiguresFigure 1 :Figure 3 :Figure 5 :Figure 6 :Figure 8 :

 13568 Figures

Figure 9 :Figure

 9 Figure 9 : Comparison between experimental data and numerical results at T=283K and P=5bar. For the 600*200µm² channel, the liquid flow rate varies from 3ml/min to 7ml/min and experimental results are drawn from (Tourvieille et al. 2013).For the 1200*400µm² channel, the liquid flow rate varies from 5ml/min to 13ml/min and experimental data are drawn from (Tourvieille 2014). Error bars represent ±20%

Figure 14 :

 14 Figure 14 : (a) Concentration profile comparison at the equatorial plane of the 10th sphere between pure diffusion conditions and convective-diffusive conditions at P=1bar, T=283K, and =30 and (b) axial evolution of hydrogen consumption flux for pure diffusion and convection-diffusion regimes at P=1bar, T=283K, and =30

liquid and liquid-solid mass transfer resistances were obtained from the physical absorption and dissolution cases respectively, is equal to 7.95.10 -5 m.s -1 . The results show a 27% deviation, as the resistances-in-series model overestimates the mass transfer coefficient. This deviation is due to the absence of a bulk concentration in the liquid film, as the film thickness is low and the flow is laminar and unidirectional, meaning that the radial mixing is poor.

was also compared to the analytical overall coefficient obtained from the equations below for a similar situation but under pure diffusion only (no fluid convection), at steady state:

Where is the heterogeneous catalytic flux condition at the catalyst surface defined in equation ( 7), and is the hydrogen equilibrium concentration given by equation ( 12). The corresponding analytical concentration profile is given by the following expression:

The corresponding analytical overall mass transfer coefficient is derived from equation ( 23), and is expressed as the following:

This and the overall mass transfer coefficient obtained by CFD are equal, indicating that in this configuration, the vertical convective down-flow does not impact the radial transport which is ensured only by pure diffusion in the liquid film thickness.

Similarly, the given by the film model gave a good agreement with the computed overall mass transfer coefficient , the relative deviation being equal to 1.3% as