Qiang Chen

Gildas Guillemot

Charles-André Gandin

Michel Bellet
email: michel.bellet@mines-paristech.fr

Numerical modelling of the impact of energy distribution and Marangoni surface tension on track shape in selective laser melting of ceramic material

Keywords: selective laser melting, laser interaction, surface tension, Marangoni effect, track profile, ceramics

Introduction

In the last decade, Additive Manufacturing (AM) has been largely developed. This innovative technology reflects the current evolution of manufacturing methods in industries. It is sometimes presented as the symbol of the next industrial revolution [START_REF] Berman | 3-D printing: The new industrial revolution[END_REF]. Indeed this new technology reduces the time gap between design and manufacturing in order to deliver parts with few finishing work or for direct end-use applications. As parts can be directly printed based on a computer-aided design model, AM processes remove the intermediate costly steps usually required in industries to manufacture pieces. In addition, few constraints are imposed on the geometry of parts, giving the possibility to produce complex shapes that would have required an assembly of several parts with traditional manufacturing technologies.

Different AM processes have been developed in recent years with applications covering a wide range of materials, from polymers to metals. Among them, Selective Laser Melting (SLM) is attractive for the manufacturing of parts requiring high level of mechanical properties, especially in the fields of aerospace and medical orthopedic for metallic alloys hard to shape with conventional cutting, forming or casting technologies. Demonstrations include aluminum [START_REF] Olakanmi | A review on selective laser sintering/melting (SLS/SLM) of aluminium alloy powders: Processing, microstructure, and properties[END_REF], stainless steel [START_REF] Akita | Defect-dominated fatigue behavior in type 630 stainless steel fabricated by selective laser melting[END_REF], titanium [START_REF] Liu | Microstructure, defects and mechanical behavior of beta-type titanium porous structures manufactured by electron beam melting and selective laser melting[END_REF] [START_REF] Schwab | Microstructure and mechanical properties of the near-beta titanium alloy Ti-5553 processed by selective laser melting[END_REF], cobalt chromium [START_REF] Averyanova | Manufacture of Co-Cr dental crowns and bridges by selective laser Melting technology[END_REF]. With SLM, processed material is totally melted by the use of a focused laser with energy concentrated on a narrow domain. Fully dense parts can thus be achieved and near net shape is possible.

However, if the process parameters are not optimized, specific defects like cracks, porosities or poor surface quality may be encountered. All these defects have detrimental consequences on final mechanical properties for direct end-use applications [START_REF] Olakanmi | A review on selective laser sintering/melting (SLS/SLM) of aluminium alloy powders: Processing, microstructure, and properties[END_REF]. This justifies efforts made in order to improve the knowledge of physical mechanisms occurring in SLM processes at the local scale as well as the development of associated modelling tools to master and control the formation of elementary tracks. Indeed, an essential difficulty in SLM is to obtain a well attached and fully dense elementary track without cracks. The track should be continuous, stable and with regular shape. Such fine tracks may reduce the formation of porosity and improve the densification of final part [START_REF] Wei | The AlSi10Mg samples produced by selective laser melting: single track, densification, microstructure and mechanical behavior[END_REF]. This requires a deep understanding of physical phenomena taking place in SLM.

The specific application of SLM to ceramic materials is basically motivated by the fact that ceramics like alumina or zirconia have outstanding mechanical strength and excellent thermal and wear resistance [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF]. However, SLM processing applied on ceramics is known for some difficulties, which explains that industrial applications are not available yet [START_REF] Wilkes | Additive manufacturing of ZrO2-Al2O3 ceramic components by selective laser melting[END_REF]. The first difficulty lies in their low energy absorption to Nd:YAG laser at the wavelength of 1064 nm. This causes high energy loss and reduces production rate. Another difficulty is the formation of cracks in solidified tracks due to the low fracture toughness of ceramics [START_REF] Sing | Direct selective laser sintering and melting of ceramics: a review[END_REF] [START_REF] Gan | Properties of selective laser melted spodumene glassceramic[END_REF]. A practical solution of crack elimination was proposed by Hagedorn et al. [START_REF] Hagedorn | Net shaped high performance oxide ceramic parts by selective laser melting[END_REF] through preheating the whole part to be just below the melting point.

As the driving force of the fusion is the laser energy, the interaction between laser and powder is of prime importance. When a laser flux irradiates a powder bed, powder particles are heated by incident radiations. However, in addition to the partial absorption process, radiations can also be reflected and transmitted by powder particles. Multiple reflections between particles usually occur, leading to an irradiated and heated region having a larger diameter than the nominal transverse section of the laser beam and being deeper than in bulk materials [START_REF] Gusarov | Modelling of radiation transfer in metallic powders at laser treatment[END_REF]. Moreover, the penetration can be quite deep for ceramic materials such as alumina (

), zirconia ( ), or a mix of those two, which can be transparent to the Nd:YAG laser radiation. Consequently, an adequate heat source model should address all these effects in order to propose an efficient modelling of the energy input. The ray tracing method at the powder particle scale was used by Wang and Kruth [START_REF] Wang | A simulation model for direct selective laser sintering of metal powders[END_REF] who found a deeper penetration for highly reflective metallic powder particles. This model was also integrated by Khairallah et al. [START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF], who pointed out that melting is not uniform and noted more heat accumulation inside powder particles compared with the homogeneous laser deposition. They showed that this effect is induced by low heat extraction through the substrate because of narrow point contact. However, the ray tracing method needs to represent powder particles and to analyze their packing structure, making the simulation complex and time consuming.

Gusarov and Smurov [START_REF] Gusarov | Radiation transfer in metallic powder bed used in laser processing[END_REF] described a model assimilating the powder to a homogeneous absorbing and scattering medium. They showed that the radial transport of radiative energy due to scattering in the powder layer can reduce the energy density at the spot center. This phenomenon is more pronounced for narrow laser beams and materials with low absorptivity.

This model was later integrated by Hodge et al. [START_REF] Hodge | Implementation of a Thermomechanical Model for the Simulation of Selective Laser Melting[END_REF] and King et al. [START_REF] King | Overview of modelling and simulation of metal powder bed fusion process at Lawrence Livermore National Laboratory[END_REF]. Other models like the one based on the Beer-Lambert law can also be found with heterogeneous absorption being taken into account [18][19].

Besides laser/material interaction, the dynamics of melt pool has a significant influence on the shape of the solidified track. Fluid flow in the melt pool is principally driven by external forces like vaporization pressure, surface tension and Marangoni effect [START_REF] Markl | Multi-scale modeling of powder-bed-based additive manufacturing[END_REF].

Vaporization often occurs in SLM as the laser energy is concentrated and the vicinity of melt pool under laser spot can easily exceed the boiling temperature. This phenomenon reduces the maximum temperature by taking away a large amount of energy, sometimes spraying out liquid drops [START_REF] Liu | Investigation into spatter behavior during selective laser melting of AISI 316L stainless steel powder[END_REF]. Surface tension has a tendency to smooth the surface of the melt pool inducing different effects on the track shape depending on the scanning speed. At low scanning speed, the melt pool is usually continuous. However, when the velocity increases, the ratio of the melt pool length to its circumference reaches the limit of Plateau-Rayleigh instability [START_REF] Gusarov | Heat transfer modelling and stability analysis of selective laser melting[END_REF]. In order to minimize surface energy, the long melt pool is interrupted into fragments, which get spheroidized due to surface tension. This so-called balling effect was observed and analyzed by Gu and Shen [START_REF] Gu | Balling phenomena in direct laser sintering of stainless steel powder: Metallurgical mechanisms and control methods[END_REF], who proposed control methods depending on the kind of balling effect. It was reproduced in numerical simulation by Khairallah et al. [START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF],

who pointed out the possibility to control the fluctuation magnitude of melt pool, and thus the balling effect by controlling the heat input over time.

In the present context, an important feature regarding surface tension is that it depends essentially on temperature. That explains the Marangoni effect which is induced by the tangential surface tension gradient along the melt pool surface, where a high tangential temperature gradient is present. This effect is considered to have a marked impact on the melt pool geometry and the temperature distribution inside it [START_REF] Chan | A two-dimensional transient model for convection in laser melted pool[END_REF]. In the context of SLM, Yuan

and Gu [START_REF] Yuan | Molten pool behavior and its physical mechanism during selective laser melting of TiC/AlSi10Mg nanocomposites: simulation and experiments[END_REF] developed a 3D finite volume model for TiC/AlSi10Mg nanocomposites considering the Marangoni effect. A wider and shallower melt pool was observed when this effect was taken into account. However, their model does not include the formation of the track shape and consequently the influence of Marangoni effect on its evolution was not investigated. Qiu et al. [START_REF] Qiu | On the role of melt flow into the surface structure and porosity development during selective laser melting[END_REF] took into account the surface tension and Marangoni flow with recoil pressure in their model with regularly packed Ti-6Al-4V powder particles. They demonstrated the significant influence of Marangoni force and recoil pressure on the melt splashing and thus the melt pool instability. Khairallah et al. [START_REF] Khairallah | Laser powder-bed fusion additive manufacturing: physics of complex melt flow and formation mechanisms of pores, spatter, and denudation zones[END_REF] used a fine-scale numerical model with randomly packed 316L steel powder particles. They concluded that the Marangoni flow induces the circulation of liquid and hence cools the laser spot region. This flow can also result in liquid spattering away from the surface due to the low liquid metal viscosity.

The authors of the present paper proposed in a previous article [START_REF] Gusarov | Radiation transfer in metallic powder bed used in laser processing[END_REF] a Finite Element (FE) model with a level set formulation dedicated to SLM. This model included a volume heat source based on the Beer-Lambert law as well as an implementation of surface tension. As a consequence, the melting of the powder, the formation of liquid droplets connected to the powder bed as well as the dynamic of the liquid-gas interface upon feeding of the melt pool with liquid droplets could be modeled. However, prediction of both convection and heat flow were limited by the absence of Marangoni convection. This is added in the model presented hereafter and gives access to the prediction of the track height variation, an output of simulations that is rarely available in the literature. In the present paper, the phenomenon of multiple reflections is considered but simplified by a simple enlargement of the Gaussian distribution. The influence of modified source distribution on the temperature field, the melt pool geometry and the shape of solidified track are compared with a non-modified Gaussian distribution. The Marangoni effect (tangential force) due to the surface tension gradient is integrated, in addition to the surface tension in the normal direction. This model is applied to simulate the formation of a single track for pure alumina with different process and material parameters. The resulting convection flow is presented and its role on the homogenization of temperature field is demonstrated. The instability of the melt pool with and without Marangoni effect is also investigated, together with their consequences on the melt pool geometry and the final track shape.

Modelling

The numerical finite element formulation developed to solve the different conservation equations on the simulation domain has been previously detailed by the authors [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramics materials[END_REF]. The following section present the essential features of this model.

Level set method

Figure 1: Schematic representation of AM by SLM as modelled using a LS method. The purple interface is between liquid ceramic and gas, whereas the yellow interface is between liquid ceramic and non-melted powder bed.

The system consists of two domains. The material domain, , is separated from the gas domain, , as shown in Figure 1. The boundary of these two domains is immersed in the system and a level set function is defined to track it. This signed distance function is taken as positive (resp. negative) in gas (resp. material) domain. A transition zone is introduced around the boundary ( ) with a half thickness in which the Heaviside function, , evolves from 0 to 1 and keeps constant outside (1 in gas and 0 in material domain):

{ [ ( )] if { | | (1)
Contrary to domain D 2 , which is made of a single gas phase, the material domain, D 1 , is multiphase. Its average volumetric properties 〈 〉 are expressed by:

〈 〉 ∑ (2)
where is the volume fraction of phase  in the domain D 1 and is the intrinsic property of phase . The average system properties { } are computed over the two-domain system by using the Heaviside function [START_REF] Desmaison | A level set approach for the simulation of the multipass hybrid laser/GMA welding process[END_REF][28]:

{ } 〈 〉 〈 〉 (3)
This averaged approach is used to define a set of thermo-physical properties including density, enthalpy, thermal conductivity and viscosity. It aims at defining thereafter a set of equations averaged over the whole system corresponding to energy, momentum and mass conservation.

Energy conservation

The transient heat transfer equation, taking into account convection and diffusion phenomena, is expressed as follows:

{ } { } { } ̇ ̇ (4)
where average quantities are the temperature, , the density, , the specific enthalpy, , and the thermal conductivity, . The convection velocity, , is the solution of the Navier-Stokes equation (see Section 2.3). The weak form of the non-linear Eq. ( 4) is solved on the whole system by use of a preexisting FE solver [START_REF] Saad | Temperature-based energy solver coupled with tabulated thermodynamic properties -Application to the prediction of macrosegregation in multicomponent alloys[END_REF]. This solvernot presented hereis based on a non-linear temperature resolution using the average enthalpy which is provided by tabulated phase transformation paths and tabulated phase properties. This general description of phase transformation allows studying any phase change during heating and cooling in the case where the phase fractions depend on the temperature only. As regards evaporation, the phenomenon is less important with ceramics as compared to the processing of metals, due to high boiling points (2977 for alumina). Hence, the liquid/gas phase transformation is actually not taken into account in the present model. This is validated a posteriori by observing that the maximum temperature found in the following simulations is below the boiling point. The next two paragraphs detail the specific expression of the right hand side terms of Eq. (4).

Surface heat loss

In the context of the level set formulation, the heat loss by radiation at the surface of the material domain ( ) is taken into account by introducing a volume heat extraction rate:

̇ (5)
where is the Stefan-Boltzmann constant, is the environment temperature and is the emissivity of the powder bed surface [START_REF] Sih | The prediction of the emissivity and thermal conductivity of powder beds[END_REF]. Here the Continuous Surface Force (CSF) method [START_REF] Brackbill | A continuum method for modeling surface tension[END_REF] is used to transform the surface heat loss into a volume heat extraction by multiplying it by the Dirac function ( ).

Laser interaction

Because of the low absorption of ceramic with respect to the 1064 nm Nd:YAG laser radiation used for processing, the volume heat source ̇ is modelled based on the Beer-Lambert law [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramics materials[END_REF]:

̇ ( ) ( ∫ ) (6)
Here is the reflection coefficient at the powder bed surface, is the nominal laser power, is the local absorption coefficient, is the radial distance to the laser beam axis. In this model, the absorption coefficient appears in the first and last terms of the right-hand side, the latter representing the attenuation in depth. Note that due to energy dispersion in the powder layer, the interaction radius is increased compared with the nominal beam radius [START_REF] Gusarov | Modelling of radiation transfer in metallic powders at laser treatment[END_REF]. In the numerical implementation, the heat source is imposed in the region , corresponding to 98.9% of the nominal power. Normalization is achieved to deliver the nominal net global power . to transform them into volume forces. Note that differs from the previously defined function, as it is computed by another level set associated to gas/liquid (purple curve) and powder/liquid (yellow curve) interfaces. The volume force associated with surface tension can then be expressed as: [START_REF] Sing | Direct selective laser sintering and melting of ceramics: a review[END_REF] where is the surface tension coefficient, is the average curvature and is the unit vector normal to the interface. In addition to the surface tension force which is normal to the interface, the volume force associated with the Marangoni surface force is expressed by: [START_REF] Gan | Properties of selective laser melted spodumene glassceramic[END_REF] where is the surface temperature gradient.

In view of FE resolution, the weak form of Eq. ( 7) can be established. Combined with Eqs. ( 8), ( 10) and ( 11) it has the following expression:

∫ ∫ ( ) ∫ ̇ ̇ ∫ ∫ ∫ [ ( ( ) ( ))] ∫ ∫ (12)
Here, braces are omitted for simplification and w denotes any vector test function. Note that the surface tension term is implemented with a semi-implicit formulation, as detailed in [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramics materials[END_REF], providing a more stable numerical resolution.

Mass conservation

The average momentum conservation is coupled with the average mass conservation expressed by the divergence of the velocity field [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramics materials[END_REF]:

{ } 〈 〉 ( 〈 〉 〈 〉 〈 〉 ) 〈 〉 〈 〉 ( ) (13)
Note that it vanishes in the gas domain D 2 . The weak form of this equation is discretized using FE and solved together with Eq. ( 12) in a mixed velocity/pressure formulation. The obtained velocity is then used to update the gas/material interface by the transport of the level set function :

where denotes the velocity at the interface between the two domains tracked by the LS method. However, the update of the level set function does not guarantee the Eikonal property, ‖ ‖ . In fact, this resolution is only valid for the interface . Therefore, the signed distance function is recomputed after the resolution of Eq. ( 14) by a geometrical method [START_REF] Shakoor | An efficient and parallel level set reinitialization method -Application to micromechanics and microstructural evolutions[END_REF] with respect to the updated position .

Mesh adaptation

In order to capture the evolving gas/material interface represented by , dynamic mesh adaptation is required. In the present model, an anisotropic mesh is used with a local size defined by a suitable metric tensor computed at each mesh node. The latter enables to follow precisely the gas/material interface as well as the evolution of the thermophysical variables (e.g. the liquid fraction). The mesh metric tensor is computed based on an error estimation combining several calculated fields [START_REF] Coupez | Metric construction by length distribution tensor and edge based error for anisotropic adaptive meshing[END_REF], particularly density, liquid fraction, temperature and heat source distribution. This method refines the regions in which the second derivative (Hessian matrix) of the selected fields is high. For instance, by tracking the density field, the gas/material interface can be followed. In addition, the melt pool boundary can be refined by tracking the liquid fraction. The advantage of this method is that it minimizes the interpolation error of tracked variables while keeping a stable number of elements. Its efficiency will be shown in the next section.

Application to SLM numerical modelling

Simulation configuration

The dimension of the whole simulation system in present applications is as illustrated in Figure 3. The bottom part corresponds to the substrate (grey color) with a height of 1 mm. A layer of powder (red color) with a thickness of 30 µm and porosity fraction of 50% is initially present on the substrate and the rest is gas (light blue). The simulations are restricted to a single pass developed in the X direction. This configuration is assimilated to the first passes of the first layer at the beginning of process when the substrate and powder are still at room temperature. To apply the model to a representative domain during construction, more complex initial temperature should be considered in order to account for the residual heat. In addition, the current substrate should be replaced by the previously consolidated part.

The material is considered as pure alumina. A laser beam with an effective power and an effective radius µm is applied. It moves from position S with coordinates to position E with coordinates at a constant scanning speed. At the beginning, the heat source is imposed progressively from 20 % to 100 % of effective power from to . At the end, the heat source is switched off but the simulation continues until all the liquid is solidified. The initial and ambient temperatures are 20 °C in the whole system. Boundary conditions for thermal and fluid mechanical resolutions are detailed in Table 1. Principal material properties are provided in Table 2 and more details can be found in Ref [START_REF] Chen | Three-dimensional finite element thermomechanical modeling of additive manufacturing by selective laser melting for ceramics materials[END_REF]. Specific additional comments are given hereafter for the viscosity and surface tension of liquid alumina, which are both considered as temperature dependent. First, the expression of the dynamic viscosity (in Pa s) as a function of temperature (in K) is obtained by fitting the results of Langstaff et al. [START_REF] Langstaff | Aerodynamic levitator furnace for measuring thermophysical properties of refractory liquids[END_REF] with an Arrhenius law:

( ) (15)
where R is the gas constant: 8.314 J mol -1 K -1 . Viscosity values are found between 42.4 mPa s at the melting temperature T m = 2054 °C and 31.9 mPa s at 2300 °C. Secondly, the surface tension is assumed to vary linearly with temperature: [START_REF] Hodge | Implementation of a Thermomechanical Model for the Simulation of Selective Laser Melting[END_REF] where is the surface tension at melting point T m . The temperature slope of the surface tension coefficient, , is taken from the literature where different measurement methods yield values varying between and N m -1 K -1 [START_REF] Paradis | Surface tension and viscosity measurements of liquid and undercooled alumina by containerless techniques[END_REF]. This variation domain being large, it is proposed here to develop SLM simulations in order to show firstly the influence of this coefficient on fluid dynamics and track formation. The parameters of presented cases are detailed in Table 3. The results of Paradis and Ishikawa [START_REF] Paradis | Surface tension and viscosity measurements of liquid and undercooled alumina by containerless techniques[END_REF] with = 0.64 N m -1 and N m -1 K -1 are used as the reference case #1. Then the extreme literature value N m -1 K -1 is applied in case #2, while in case #3 an artificial positive value aims at studying the influence of the sign of . The effect of the scanning speed and the effective interaction radius on the temperature distribution and melt pool stability are investigated in cases #4 and #5. The time step is set to 2 µs for cases #1 to #3 and 1.5 µs for cases #4 and #5. The number of mesh elements increases from around 800 000 to 1 200 000 from the beginning to the end of simulations in order to have a good representation of the track surface. 3). Note that the gas domain is removed in all figures so as to see the material domain. (see attached video SingleTrack_Case_1)

Results and discussion

Temperature field with isocontours of the reference case #1 is shown in Figure 4 (a)

when the laser reaches position X = 1.4 mm (at time t = 4 ms), i.e. approximately at middistance of the total track length (see attached video SingleTrack_Case_1). The melt pool is elongated behind the laser spot with a length of about 240 µm. The isocontours show an evolution of the tail shape during cooling. It changes from a convex shape at high temperature to a non-convex shape at low temperature, with an enlargement into the powder bed. At high temperature, for example for T > 1100 °C in the region corresponding to the red window used to identify the Figure 4 (c), the solidified track is hotter than the powder on its lateral sides.

However, when it is cooled to 700 °C, the track becomes colder than the neighboring powder.

This evolution is caused by the insulant effect of powder due to its lower thermal conductivity ( 0.2 W m -1 K -1 ) as a porous material [START_REF] Gusarov | Radiation transfer in metallic powder bed used in laser processing[END_REF][START_REF]Dry air properties[END_REF].

Moreover, one significant phenomenon in SLM is the high temperature gradient induced by rapid heating. Generally, the temperature gradient ahead of the laser spot is much higher than that behind it, as the heat input at front cannot be rapidly conducted and transferred. The average temperature gradient in the track just behind the melt pool is about K m -1

(estimated with isocontours T = 1500 °C and T = 1900 °C), while it is about K m -1

(estimated with isocontours T = 2014 °C and T = 2300 °C) in the melt pool, as the convection tends to homogenize the temperature field. The melt pool shape in Figure 4 (b) shows a long tail and a penetration of about 50 µm in the substrate. It is largely influenced by the laser energy distribution (depending on the laser power, reflection, interaction radius and material absorption) and the scanning speed. Consequently, a desired melt pool penetration in substrate can be obtained by a good control of the process parameters. At the front of the melt pool, liquid droplets form by powder melting and spheroidization under the surface tension effect.

These droplets lead to the discontinuity of the liquid zone. They fall into the melt pool under gravity and cause instabilities or wave propagation in the melt pool. This phenomenon is captured by the use of the LS method and the efficient mesh adaptation approach previously detailed. The latter is an important issue in the present FE method and the mesh size control is shown in Figure 4 (c) and (d). By tracking two solution fields (density and liquid fraction), the melt pool with high curvature is well refined, especially for the melted powder region with a minimum mesh size around 0.5 as shown in Figure 4 (d). However, once the track is solidified, liquid fraction has no more effect on the mesh refinement and a larger mesh size is retrieved (Figure 4 (c)). In addition, elements are elongated in the X direction as low density gradient variation is encountered in this direction, leading to an anisotropic mesh. The Marangoni effect on the temperature distribution can be investigated by comparing cases #1 to #3 in Figure 5 (laser at X = 1.4 mm and t = 4 ms). With a higher negative value of in case #2 (Table 3), the melt pool (red color) is enlarged both in width and length.

This also occurs for the temperature isocontours above the liquidus temperature.

Consequently, the isocontours in the solidified track are also extended. By contrast, a positive value of (case #3) reduces the melt pool, leading to a narrower and shorter pool shape.

This effect can result in the necking of isocontours, for example at T = 2300 °C, which is also indicated in the longitudinal cut view. The consequence of the necking is the interruption of isocontours, such as the one corresponding to T = 2500 °C. By comparing cases #1 to #3, it is shown that the melt pool becomes unstable when the Marangoni coefficient varies from the negative to the positive value. However, the Marangoni effect on the depth of melt pool is not prominent for these cases #1 to #3.

As a second step in simulation strategy, the scanning speed has been increased to 500 mm s -1 . This leads to a reduction in the deposited energy and consequently to a smaller melt pool, as shown in case #4 (Figure 5). In addition, the high scanning speed results in the discontinuity of temperature isocontours. On the other hand, with a higher interaction radius as proposed in case #5 (Figure 5), powder is heated in a wider region. Correspondingly, more powder is melted at the melt pool front and larger droplets are formed, leading to a more unstable melt pool. In both cases #4 and #5, the melt pool is much shallower than that in case #1 and it can just melt the substrate surface. Temperature isocontours are the same as in Figure 5.

The temperature distribution in the melt pool is largely influenced by the convection flow. Top and transverse cut views are presented in Figure 6. Several phenomena contribute to fluid flow. Firstly, the melted powder at the front of the melt pool forms droplets under the effect of surface tension. Discontinuous fall events of these droplets in the pool under the effect of gravity create waves that propagate toward the rear of the melt pool. This explains the velocity distribution for case #1 shown in Figure 6 where the surface velocity field is essentially oriented toward the rear of the melt pool near section A-A (green arrows) after the fall event. This convection flow is also affected by the Marangoni effect, leading to spreading of the melt toward the two sides of the melt pool. The latter is more evident in case #2 when the magnitude of is increased. Because , the convection flow is induced from hotter region (low surface tension) to colder region (high surface tension). This flow drives heat outwards and extends the boundary of the melt pool, as well as the temperature isocontours. The transverse cut view shows clearly that two opposite convection cells are formed. At the surface of melt pool, centrifugal flow is driven from the center (hotter region) to lateral sides (colder region), leading to a wider melt pool. In the transversal cutting plane, the flow is upward to the melt pool surface. These convection cells result in the transport of fluid and increase heat exchange, thus changing the isocontours (T = 2300 °C and T = 2500 °C). In contrast to case #2, a positive with the same magnitude as in case #3 shows convection flow in the opposite direction. The flow is driven from colder to hotter regions, as surface tension coefficient increases with temperature. This leads to the reduction of melt pool surface and temperature isocontours, or even discontinuity (T = 2500 °C). Two convection cells in opposite directions compared to those of case #2 are observed in the transverse cut view. Fluid flows inward at the surface of the melt pool and downward along the longitudinal symmetry plane. However, the convection cells are limited to a weak penetration into the substrate, probably because of the relatively high viscosity of liquid alumina (compared to the viscosity of liquid metals). Consequently, the melt pool depth is almost unchanged, comparing cases #2 and #3 with the reference case #1.

One of the possible consequences of the melt pool dynamics is the irregular profile of the solidified tracks, with interruption of track or balling effect. The surface heights of the final tracks in cases #1 to #5 are shown in Figure 7, together with isocontours. A transition region at the beginning is observed where a swelling occurs. However, the domain of interest is the quasi-steady state which follows. The regularity of the track profile is perceptible by the color difference along the track and by the discontinuity of isocontours. The high altitude regions (hills) are highlighted by small circles and the necks of isocontours are the low altitude regions (valleys). Compared to case #1, the number of hills and valleys decreases in case #2, meaning a more stable surface height when is decreased. By contrast, the increase of leads to more hills and valleys, thus a higher track irregularity in case #3. The Marangoni effect on the track regularity can be explained by the convection flow in Figure 6.

If we consider a melt pool with two successive liquid hills and a valley between them, convection flow with negative tends to drive fluid from hills to valleys, thus reducing the height difference between them. Conversely, a positive value of can promote this difference and enhance irregularities. This means that a high magnitude negative is favorable to obtain a smoother track surface, and may be helpful to avoid balling effect. In the simulation of processing stainless steel with SLM, Khairallah et al. [START_REF] Khairallah | Mesoscopic simulation model of selective laser melting of stainless steel powder[END_REF] obtained the height variation of solidified track surface. Although the Marangoni effect was not taken into account, they guessed that negative could smooth the track surface. This remark is confirmed by the present work. Case #4 with increased scanning speed shows more hills and valleys, thus more surface instabilities. This reveals the possible tendency to the balling effect at high scanning speed, which is coherent with the results of Gu and Shen [START_REF] Gu | Balling phenomena in direct laser sintering of stainless steel powder: Metallurgical mechanisms and control methods[END_REF]. In addition, the solidified track is narrower than in case #1. A larger track is obtained in case #5 with a higher laser interaction radius, leading to an increase of melted powder. Its merge into the existing melt pool leads to instabilities and to a more marked irregular profile of the solidified track. Track regularity is illustrated by the height profiles in Figure 8, taken in the longitudinal median plane of cases #1 to #5. Two indicators for profile regularity are selected among those adopted for roughness evaluation along a linear profile of length L. The first one is the total roughness with x the spatial coordinate along the linear profile and z(x) the altitude difference with respect to the average altitude over the whole profile. The second one is the average quadratic measure for roughness: √ ∫ . Less height variation with a R q value of 0.38 µm is formed in case #2 compared to the value of 0.65 µm in case #1. As already noted in Figure 7, case #3 presents a more irregular track shape with a higher R q (0.80 µm) and a higher R t (3.7 µm). In cases #4 and #5, a higher interaction radius or a higher scanning speed increase irregularities in track height, leading to higher values of roughness indicators. Finally, the volume and the surface area of the melt pool (including droplets) are plotted in Figure 9 (a) and (b), in order to study its stability in cases #1 to #5. The transition region at the beginning of scanning is clearly shown by the peaks of melt pool volume and surface area.

The stable regime begins when the laser beam arrives at position X 0.7 mm. The melt pool volume and surface area can be compared from this position to the switch off position of the laser beam at X = 2.8 mm. For a given case, the scattering of the surface area is much higher than that of the melt pool volume, as the whole liquid volume is conserved while this is not the case for surface area during the fall of droplets. These variations may be considered as indicators of the instability of the melt pool. Comparing with the reference case (#1), the average volume is slightly increased in case #2 and decreased in case #3. However, the difference is much less significant compared with the surface area in Figure 9 (b). Due to the centrifugal surface flow along the melt pool surface in case #2, the surface area is largely increased. By contrast, the centripetal surface flow in case #3 leads to the reduction of the area of the melt pool surface. Both cases #4 and #5 have a much smaller melt pool volume than case #1, which can be indicated by their shallower melt pool as shown in Figure 5 (longitudinal cut view). The melt pool in cases #4 and #5 has almost the same surface area while its volume in case #5 is higher than that in case #4 as the melt pool is deeper in case #5 (Figure 5).

Conclusions

In the present approach, the effects of surface tension and Marangoni force are considered in a finite element model with level set method at the track scale, expressed in terms of continuum mechanics and heat transfer. A compressible behavior of powder is introduced to model powder shrinkage when melted. The important influence of the Marangoni effect on temperature distribution and melt pool shape is discussed by using different values. It is demonstrated that a negative value leads to the expansion of the melt pool, while a narrower and shorter melt pool can be obtained with a positive value. Similar evolutions are observed for temperature isocontours. These observations are discussed and explained considering evolution of tangential forces induced by Marangoni effect and their direction. However, its influence on melt pool depth is not the main effect as is commonly observed for metals, due to the high viscosity of liquid ceramic.

In addition, the Marangoni flow has an effect on the regularity of the profile of solidified track. A high negative value of tends to smooth the surface and avoid the fragmentation of the track and associated balling effect. The scanning speed and interaction radius are also investigated. Results show that the melt pool becomes smaller and discontinuous with the increase of scanning speed or laser/material interaction radius. Future work will focus on the comparison with experimental observations currently in progress and the quantitative validation of the present model. Then the construction of process map under different working conditions with different track morphologies will be carried out to better control the process. Moreover, the integration of solid mechanics should complete this model by predicting stresses in the freshly solidified track and in its vicinity, which is important regarding the definition of process windows to avoid cracking during construction.

Figure 2 :

 2 Figure 2: Energy distribution with different values of r int and  (Eq. (5)) under the same R and P L , isocontours corresponding to ̇ and ̇ , where ̇ .

  function

Figure 3 :

 3 Figure 3: Configuration for unidirectional single pass simulation. A layer of powder (red) with 30 µm thickness lies at the top of a substrate (grey) with 1 mm height. The top gas domain (light blue) is displayed with a transparent color in order to show the scanning trajectory and the solidified track on the substrate.

Figure 4 :

 4 Figure 4: (a) Temperature field (top view) and isocontours (black lines) from 300 to 2300 °C with 400 °C interval. The liquidus line is also presented (red contour) for T = 2104 °C; (b) melt pool shape with the surface position of substrate (black line); (c)-(d) mesh in the solidified track and at the melt pool front when the laser arrives at position X = 1.4 mm (t = 4 ms) for case #1 (Table3). Note that the gas domain is removed in all figures so as to see the

Figure 5 :

 5 Figure 5: Top (left column) and longitudinal cut view (right column) of temperature field for cases #1 to #5 when the laser arrives at position X = 1.4 mm. Temperature isocontours (black) range from 300 to 2700 °C with 200 °C interval and the red isocontour corresponds to the melt pool. The vertical white line at left is for comparing the melt pool length and the horizontal one at right indicates the surface of substrate.

Figure 6 :

 6 Figure 6: Velocity field (arrows colored by magnitude) at the surface of melt pool and the section for cases #1 to #3 when the laser arrives at position X = 1.4 mm for t = 4 ms.

Figure 7 :

 7 Figure 7: Height (mm) of solidified tracks with respect to the substrate surface and isocontours for cases #1 to #5. The image of tracks is scaled by a factor 2 in the Y direction for better visualization.

Figure 8 :

 8 Figure 8: Height profiles in the longitudinal median plane for cases #1 to #5 between X = 0.5mm and X = 2.5 mm (i.e. the region between the two black lines in Figure7). Indicators for track profile regularityas defined in the textare given on the right.

Figure 9 :

 9 Figure 9: Evolution of (a) melt pool volume and (b) surface area for cases #1 to #5.

Table 1 :

 1 Boundary conditions for heat transfer and fluid flow resolutions.

	Properties	Symbol	Value	Unit	Ref.
	Alumina density		3800	kg m -3	[34]
	Gas density		1.3	kg m -3	[35]
	Specific enthalpy of alumina		Fig.5 (a) in Ref [19]	J kg -1	[36]
	Specific heat of gas		1000	J kg -1 K -1	[35]
	Thermal conductivity of alumina			W m -1 K -1	[37]
	Thermal conductivity of gas		0.024	W m -1 K -1	[35]
	Gas viscosity			Pa s	[35]
	Powder absorption	〈 〉	4	mm -1	[19]
	Dense matter absorption	〈 〉	3	mm -1	[19]

Table 2 :

 2 Material properties of alumina and gas.

		Parameters		Case
	Symbol	Unit	#1 (reference)	#2	#3	#4	#5
		N m -1 K -1	-8.2	-48	+48	#1
	v L	mm s -1	300	#1		500	#1
	r int	µm	35		#1		50

Table 3 :

 3 Set

of process and physical parameters used for the simulation cases. Values in bold red color highlight the modified parameters with respect to the reference case.

Acknowledgement

This work has been conducted within the framework of the CEFALE project, part of the ACLAME program funded by the Institute CARNOT M.I.N.E.S (Paris, FR). The authors would like to thank Liliana Moniz Da Silva Sancho, Christophe Colin, Jean-Dominique Bartout and Marie-Hélène Berger from the MINES ParisTech, Centre des Matériaux (Evry, FR) for information and discussion regarding AM by SLM applied to ceramic materials.

The energy distribution with different values of and under the same and is illustrated in Figure 2. The increase of always leads to shallower isocontours, while it becomes narrower (e.g. isocontour 0.4) near the center but wider (e.g. isocontour 0.1) far away from the center. On the other hand, the decrease of  always results in narrower isocontours, while its influence on the isocontour depth is not monotonous. With smaller values for  coefficient, the region with high energy becomes shallower (e.g. isocontour 0.4)

and the one with low energy becomes deeper (e.g. isocontours 0.1 and 0.2).

Momentum conservation

The velocity is computed in the whole system by the resolution of the average momentum conservation equation:

where u is the velocity, is the stress tensor, is the gravity vector, and are the volume forces corresponding to surface tension and Marangoni forces, respectively. These contributions are detailed hereafter.

Behavior law

As the powder is described as a continuous medium, the apparent density of the material domain D 1 varies strongly during SLM. As a consequence, shrinkage takes place during consolidation from powder to compact material. Thus a compressible Newtonian constitutive law is used in the whole system including material and gas domains:

where is the pressure, is the dynamic viscosity, is the deviatoric part of , is the identity tensor and ̇ is the strain rate tensor. Here { } is a term related to the shrinkage rate deduced from the mass conservation equation as described below.

Surface forces implemented as volume forces

Both the surface tension and Marangoni force are imposed at the gas/liquid and powder/liquid interface, shown as the purple and yellow curves in Figure 1, respectively. As these interfaces are immersed in the system, the surface forces are multiplied by a Dirac


EXPERIMENTAL OPTIMIZATION OF AN INNOVATIVE BIOGAS UPGRADING PROCESS ADAPTED TO THE AGRICULTURAL CONTEXT

Keywords: upgrading, methane, absorption, membrane contactor, process intensification

The development of the large European biomethane potential is facing the challenge of the biogas upgrading technologies adaptation to the agricultural context (small and medium size anaerobic digestion units, rural environment). Membrane contactors are an opportunity to adapt the robust Water Scrubbing (WS) process and to decrease the high specific investment costs reported for available technologies. Microporous polypropylene (PP) hollow fiber membrane contactors were combined in an absorption/desorption closed-loop using water as the absorbent. A Design of Experiments methodology highlighted the gas and liquid flowrates, the inlet biogas composition and the absorption gas pressure to be the most significant parameters impacting biomethane quality with a low level of interactions. From the subsequent optimization, gas-grid quality biomethane (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 97.6%) was produced with a moderate recovery rate (𝑅 𝐶𝐻 4 = 83.7%). The limited mass transfer efficiency may be compensated by a larger membrane area. Alternatively, the addition of a methane recycling loop, similarly to the conventional WS process, significantly improved the methane recovery rate (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 95.9% / 𝑅 𝐶𝐻 4 = 93.5%). These experimental results were confirmed by on-site upgrading of raw farm biogas.

INTRODUCTION

In the last decades, biogas has been considered as an energy source to build a sustainable energy mix and substitute fossil fuels. The Anaerobic Digestion (AD) of organic matter (agricultural waste, industrial waste, biowaste, sewage sludge, energy crops…) produces biogas essentially composed of methane and carbon dioxide [START_REF] Wellinger | The biogas handbook -Science, production and applications[END_REF]. Its calorific content can be converted through several pathways: heat production, electricity production or conversion to biomethane [START_REF] Pöschl | Evaluation of energy efficiency of various biogas production and utilization pathways[END_REF]. The latter offers a substitute to fossil natural gas: it can be stored and transported in the gas-grid infrastructure and/or utilized as vehicle fuel [START_REF] Budzianowski | Biomethane storage: Evaluation of technologies, end uses, business models, and sustainability[END_REF].

Several processes are available for biogas upgrading and relies on different separation principles [START_REF] Yang | Progress and perspectives in converting biogas to transportation fuels[END_REF], [START_REF] Bauer | Biogas upgrading -Review of commercial technologies[END_REF]: gasliquid absorption (physical such as water or organic scrubbing, chemical such as amine scrubbing), adsorption (PSA), gas membrane permeation or cryogeny. Among them, Water Scrubbing (WS) offers the advantages of a chemical-free process with a limited energy consumption 0.212 𝑘𝑊ℎ. 𝑁𝑚 𝑏𝑖𝑜𝑔𝑎𝑠 -3 [START_REF] Xiao | CO2 Removal from Biogas by Water Washing System[END_REF]. Its robustness is well-adapted to the agricultural context with non-technician operators [START_REF] Rotunno | Energy and economic analysis of a water scrubbing based biogas upgrading process for biomethane injection into the gas grid or use as transportation fuel[END_REF], [START_REF] Benizri | Epuration du biogaz à la ferme : EPUROGAS, une solution énergétique et économique d'avenir. -Etude expérimentale et modélisation d'un procédé d'absorption de dioxyde de carbone avec de l'eau sous pression à une échelle industrielle[END_REF]. These technologies aim at producing a gas-grid quality biomethane (𝑦 𝐶𝐻 [START_REF] Yang | Progress and perspectives in converting biogas to transportation fuels[END_REF] 𝑜𝑢𝑡 > 97% for H gas in France)

with the highest methane recovery rate 𝑅 𝐶𝐻 4 (generally above 97%). Europe and France specifically hold a massive potential for biomethane production [START_REF] Thrän | European biomethane potentials[END_REF], [START_REF] Cowi | State of the Art on Alternative Fuels Transport Systems in the European Union[END_REF]. Within it, agricultural waste from small and medium size farms represent a significant share [START_REF] Aebiom | A Biogas Roadmap for Europe[END_REF]. Though, the development of biogas upgrading technologies is limited due to their rising specific investment costs as the size of the AD unit decreases [START_REF] Bauer | Biogas upgrading -Review of commercial technologies[END_REF]. Gas permeation and WS are the leading technologies on this specific market. WS is particularly adapted to a rural context due to its tolerance to biogas impurities (such as volatile organic compounds or hydrogen sulfide).

Hollow fiber membranes, a technology derived from the artificial lung technology, is a gas-liquid contacting device which offers several advantages [START_REF] Luis | Non-dispersive absorption for CO2 capture: from the laboratory to industry[END_REF], [START_REF] Zhao | Status and progress of membrane contactors in post-combustion carbon capture: A state-of-the-art review of new developments[END_REF]. Its modular nature offers an easy process design while its compactness makes it visually more attractive than conventional packed columns. The hydrophobic membrane confines both gas and liquid phases to their respective compartments and provides an operational flexibility (no flooding or channeling issues).

As another advantage, the use of membrane contactors in both absorption and desorption steps avoids the depressurization of the solvent generally needed for air stripping desorption occuring in packed columns; this represents a significant energy saving as the solvent recompression is considered to account for 20-30% of the process energy requirement [START_REF] Chabanon | Gasliquid separation processes based on physical solvents: opportunities for membranes[END_REF].

These advantages make it an attractive solution for gas-liquid absorption processes. Most of the studies focused on post-combustion carbon capture [START_REF] Luis | Non-dispersive absorption for CO2 capture: from the laboratory to industry[END_REF]- [START_REF] Gomez-Coma | Carbon dioxide capture by [emim][Ac] ionic liquid in a polysulfone hollow fiber membrane contactor[END_REF] but a few of them experimentally explored its potential for biogas [START_REF] Mcleod | Quantifying the loss of methane through secondary gas mass transport (or "slip") from a micro-porous membrane contactor applied to biogas upgrading[END_REF]- [START_REF] Atchariyawut | Separation of CO2 from CH4 by using gas-liquid membrane contacting process[END_REF] or natural gas upgrading [START_REF] Atchariyawut | Separation of CO2 from CH4 by using gas-liquid membrane contacting process[END_REF]- [START_REF] Rahim | Absorption of CO2 from natural gas using different amino acid salt solutions and using hollow fiber membrane contactors[END_REF]. The largest membrane areas reported were tested with PP hollow fibers (highly available for processability and cost considerations).

The mass transfer resistance is typically described by a resistance in serie model divided into the gas, membrane and liquid phases [START_REF] Fougerit | Gas-liquid absorption in industrial cross-flow membrane contactors: Experimental and numerical investigation of the influence of transmembrane pressure on partial wetting[END_REF]. In the case of pure CO2 absorption, the limiting mass transfer resistance was demonstrated to be located in the liquid phase [START_REF] Fougerit | Gas-liquid absorption in industrial cross-flow membrane contactors: Experimental and numerical investigation of the influence of transmembrane pressure on partial wetting[END_REF], [START_REF] Mansourizadeh | Experimental study of CO2 absorption/stripping via PVDF hollow fiber membrane contactor[END_REF]. Similarly, in the case of biogas upgrading, McLeod proved that the carbon dioxide absorbed flow 𝑄 𝐶𝑂 2 𝑎𝑏𝑠 in water was driven by the liquid velocity when a physical absorbent was used [START_REF] Mcleod | Quantifying the loss of methane through secondary gas mass transport (or "slip") from a micro-porous membrane contactor applied to biogas upgrading[END_REF].

Considering only absorption, 85% biomethane was produced but the methane slip was already significant (𝑅 𝐶𝐻 4 = 94.78 %). Using a 1M NaCl solution instead of pure water, 𝑅 𝐶𝐻 4 was increased to 96.00 % thanks to the salting out effect. In a process development perspective, only one recent study carried out experimental work to implement an absorption/desorption cycle applied to biogas upgrading, as presented by Teplyakov [START_REF] Teplyakov | Computer design of recycle membrane contactor systems for gas separation[END_REF]. With water as the absorbent, Kim produced a gas-grid quality biomethane (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 98%) though with a lower recovery rate (𝑅 𝐶𝐻 4 = 85 %) [START_REF] Kim | Practical designs of membrane contactors and their performances in CO2/CH4 separation[END_REF]. Coupling absorption and desorption steps, 𝑅 𝐶𝐻 4 was further decreased to 75%. Despite biogas upgrading standards were not reached, the process potential was confirmed by a stability over 50h.

The main concern regarding the industrial application of hollow fiber membrane modules to gas-liquid absorption processes is the membrane wetting phenomenon [START_REF] Kreulen | Determination of mass transfer rates in wetted and non-wetted microporous membranes[END_REF]. Pores of the hydrophobic membrane are partially filled with solvent, thus increasing the membrane mass transfer resistance. Absorbed CO2 flux was divided by 2 after 800 hours of operation with 30%wt monoethanolamine (MEA) solution [START_REF] Chabanon | Membrane Contactors for Postcombustion Carbon Dioxide Capture: A Comparative Study of Wetting Resistance on Long Time Scales[END_REF]. The overall mass transfer resistance was increased by 21-53 % in the case of pure CO2 absorption in water after some hours of operation [START_REF] Mavroudi | A study of mass transfer resistance in membrane gas-liquid contacting processes[END_REF]. The compatibility of the solvent with the membrane material is an essential selection criterion to limit this phenomenon [START_REF] Mosadegh-Sedghi | Wetting phenomenon in membrane contactors -Causes and prevention[END_REF]. Current research on membrane materials is addressing this wetting issue. The deposition of a thin dense layer or the development of liquid-filled membrane are potential solutions [START_REF] Zhao | Status and progress of membrane contactors in post-combustion carbon capture: A state-of-the-art review of new developments[END_REF], [START_REF] Gomez-Coma | Membrane modules for CO2 capture based on PVDF hollow fibers with ionic liquids immobilized[END_REF]. Another strategy consists in increasing the membrane hydrophobicity by membrane surface treatment [START_REF] Abdulhameed | Carbon dioxide capture using a superhydrophobic ceramic hollow fibre membrane for gas-liquid contacting process[END_REF], [START_REF] Lee | Effect of hydrophobic modification on carbon dioxide absorption using porous alumina (Al2O3) hollow fiber membrane contactor[END_REF].

The results from these previous studies indicate the potentialities of hollow fiber membrane modules for biogas upgrading using water as the absorbent. This paper aims at experimentally optimizing the coupled absorption/desorption process using a liquid closed-loop. Experiments were run in a partial wet mode to account for the membrane wetting limitation. The effects of operating parameters (flowrates, pressure, inlet biogas composition) have been studied through a Design of Experiments. From these initial results, the influence of the membrane area at the absorption step was tested in an open-loop configuration. At last, a new process configuration mimicking the intermediate flash tank from the WS process with a methane recycling loop was suggested and assessed. The membrane module was selected from a range of available industrial modules (Liqui-Cel® Extra-Flow) in order to easier the foreseen process upscaling step. The geometric features of the modules and the data on both membrane types used in this study are respectively presented in Table I and Table II. This module range is built with a diverting baffle at its center to enhance the mass transfer efficiency. As the liquid flows in the shellside, the flow pattern is in between the countercurrent parallel flow and the cross-flow configurations. The fiber bundle of the 2.5"x8" module is estimated at 10 200 fibers.

The X-50 fiber shows a higher porosity compared to the X-40 type (40-45% against 20-25%) and a lower thickness (80 µm against 100). These features are likely to enhance the gas-liquid mass transfer. Thus, a X-50 module is used for the absorption step. The X-40 type is selected at the desorption step for its higher mechanical resistance: it indeed faces a 13-fold higher transmembrane pressure.

Experimental set-up for biogas upgrading

A modular set-up was designed to evaluate the biogas upgrading process with different configurations. The reference configuration used in the Design of Experiments (DoE) is detailed in Figure I. Carbon dioxide from biogas is dissolved into reverse osmosis water in the absorption module. The CO2-rich liquid is then degassed by applying a vacuum in the desorption module. Vacuum is preferred to air stripping firstly to avoid air traces in the outlet biomethane but also to produce a non-diluted offgas stream. The vacuum pump (Vacuubrand MZ 2C NT VARIO) establishes a pressure level in the range 50-500 mbar a.

The inlet gas flowrate and its composition are controlled by a 2-way mixer equipped with mass flow controllers (FC 41 and FC 42, Bronkhorst In-Flow CTA, ± (1 %RD + 1 %FS)). Pure gases (CO2, CH4 -Air Liquid AlphaGaz 1 standard) were used to produce synthetic mixtures. After CO2 removal, the outlet gas flowrate is measured in another mass flow controller (FC 43, Bronkhorst Low-ΔP-Flow, ± 1 %FS). The process gas pressure is regulated with a PID controller connected to the solenoid valve of the outlet mass flow controller. The degassed gas flow, named offgas, is measured in a mass flowmeter (Bronkhorst In-Flow CTA, ± (1 %RD + 1 %FS)). It takes about 5 minutes for the gas flows and pressure to stabilize at their setpoints. As mass flowmeters were calibrated using air, they require the use of corrective factors depending on gas composition analysis which induce an additional uncertainty ±2%. The DoE investigated the available operating conditions range presented in Table III.

𝑜𝑢𝑡

and the methane recovery rate 𝑅 𝐶𝐻 4 . The process is evaluated in its steady state regime. Gas composition analyses were performed 20 minutes after the process parameters reached their setpoints. Similar stabilization times are reported in the literature [START_REF] Mcleod | Quantifying the loss of methane through secondary gas mass transport (or "slip") from a micro-porous membrane contactor applied to biogas upgrading[END_REF], [START_REF] Kim | Practical designs of membrane contactors and their performances in CO2/CH4 separation[END_REF]. A second gas composition analysis was performed 15 minutes after the first one to confirm the steady state regime had been reached.

The methane recovery rate 𝑅 𝐶𝐻 4 can be derived from inlet and outlet gas flowrates and compositions (Eq. 1).

𝑅 𝐶𝐻 4 = 𝑄 𝐶𝐻 4 𝑜𝑢𝑡 𝑄 𝐶𝐻 4 𝑖𝑛 = 𝑦 𝐶𝐻 4 𝑜𝑢𝑡 𝑄 𝑜𝑢𝑡 𝑦 𝐶𝐻 4 𝑖𝑛 𝑄 𝑖𝑛 Eq. 1
Though, the uncertainties are stronger on the gas flowrates than on the gas composition analysis. Therefore, a method based on the composition of the three gas flows was preferred. The combination of two mass balances (total gas and methane only) leads to:

𝑅 𝐶𝐻 4 = 𝑦 𝐶𝐻 4 𝑜𝑢𝑡 𝑦 𝐶𝐻 4 𝑖𝑛 (𝑦 𝐶𝐻 4 𝑖𝑛 -𝑦 𝐶𝐻 4 𝑜𝑓𝑓 ) (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 -𝑦 𝐶𝐻 4 𝑜𝑓𝑓 ) Eq. 2
The efficiency of the CO2 mass transfer is assessed for the open-loop experiments. First, the logarithmic mean of the transfer potential ∆𝑥 𝑀𝐿 is calculated.

∆𝑥 𝑀𝐿 =

(𝑥 𝑜𝑢𝑡, * -𝑥 𝑜𝑢𝑡 ) -(𝑥 𝑖𝑛, * -𝑥 𝑖𝑛 ) ln(𝑥 𝑜𝑢𝑡, * -𝑥 𝑜𝑢𝑡 ) -ln(𝑥 𝑖𝑛, * -𝑥 𝑖𝑛 ) Eq. 3

𝑥 𝑜𝑢𝑡, * and 𝑥 𝑖𝑛, * are derived from gas composition analysis and values from Henry coefficients [START_REF] Sander | Compilation of Henry's law constants, version 3.99[END_REF]. 𝑥 𝑜𝑢𝑡 is calculated by a CO2 mass balance over the absorption unit (Eq. 4). 𝑥 𝑖𝑛 is reasonably assumed to be zero in the open-loop configuration (initially gas-free inlet absorbent).

𝑥 𝑜𝑢𝑡 = 𝐺 𝑖𝑛 𝐿 𝑦 𝐶𝑂 2 𝑖𝑛 - 𝐺 𝑜𝑢𝑡 𝐿 𝑦 𝐶𝑂 2 𝑜𝑢𝑡
Eq. 4

A constant liquid overall mass transfer coefficient 𝐾 𝐿 0 is then derived from experimental data. It is related to 𝐾 𝐿 coefficient expressed in 𝑚. 𝑠 -1 by Eq. 6.

𝐾 𝐿 0 = 𝐿(𝑥 𝑜𝑢𝑡 -𝑥 𝑖𝑛 ) 𝑎 ∆𝑥 𝑀𝐿 𝑆𝑍 Eq. 5

𝐾 𝐿 = 𝐾 𝐿 0 𝑐 𝑡𝑜𝑡 𝑙 ≅ 𝐾 𝐿 0 𝑐 𝐻 2 𝑂

𝑙

Eq. 6

3 DESIGN OF EXPERIMENTS

Influence of process parameters

In a first approach, the process parameters and their expected influence on the performances have been summarized in Table IV. The flows 𝑄 𝑔 and 𝑄 𝑙 are the main parameters of the gas-liquid absorption processes, typically described by their ratio 𝑄 𝑔 𝑄 𝑙 ⁄ . The gas pressure 𝑖𝑛 necessarily impacts the process performances. At industrial scale, the biogas quality is an input parameter which depends on the substrates and on the AD process. Among the listed parameters, five are selected for further investigation for their supposed influenced on the biogas upgrading process: both gas and liquid flowrates 𝑄 𝑔 and 𝑄 𝑙 , gas absorption and vacuum pressures 𝑃 𝑔 and 𝑃 𝑣𝑎𝑐 , and inlet biogas composition 𝑦 𝐶𝐻 4 𝑖𝑛 .

Design of Experiments

The DoE methodology was adopted to compare the relative influence of the process parameters and their interactions in a range of operating conditions. Five parameters were selected and three levels were defined to investigate the quadratic influences and interactions (Table V). The temperature of the absorbent was measured but not regulated (T l = 25°𝐶 ± 3°𝐶). Biomethane quality 𝑦 𝐶𝐻 [START_REF] Yang | Progress and perspectives in converting biogas to transportation fuels[END_REF] 𝑜𝑢𝑡 was the selected process response for its low uncertainty (±0.3% 𝑣𝑜𝑙.).

A central composite design was used to reduce the number of experiments (243 tests initially). The central composite design is divided in 3 matrix:

-A factorial design with 2 levels (-1 and +1): 32 tests; -A series of tests at the center of the domain: 𝑁 0 tests; -A series of tests on the axis of the domain (-θ and +θ for each parameter): 10 tests.

The value of θ=2.378 is chosen from tabulated values to satisfy the optimality criteria: rotability, orthogonality and uniform precision. The criteria of orthogonality and uniform precision cannot be satisfied at the same time. Therefore, the design is qualified of quasi-orthogonal.

For a given number of parameters (here k = 5), the number of tests 𝑁 0 at the center of the domain is defined according to the criterion that must be satisfied (17 and 10 respectively for orthogonality and uniform precision criteria). In practice, 15 tests were performed at the center of the domain.

Statistical processing

A quadratic model with interactions was used to interpret the experimental data. Three models are successively constructed in a MATLAB program: -Model1 which is the raw model with all the coefficients considered; -Model2 which is model1 with only the same significant coefficients conserved; -Model3 which is an optimization of the coefficients with the significant parameters only.

To determine whether the effect of the parameter is significant, the ratio of its effect 𝑎 𝑃 𝑖 against the model standard error 𝑠 𝑚𝑜𝑑𝑒𝑙 is compared to a Student law coefficient with 𝑁 -𝐶 degrees of freedom (N number of experiments, C number of coefficients in the model).

| 𝑎 𝑃 𝑖 𝑠 𝑚𝑜𝑑𝑒𝑙 | > 𝑡 1-𝛼 2 ⁄ ,(𝑁 -𝐶) 𝑠𝑡𝑢𝑑𝑒𝑛𝑡

Eq. 7

The value of α determines the confidence interval. 𝛼 = 0.05 is a typical choice which describes the confidence interval at 95 %.

The same student coefficient and the mean squared error 𝜎 𝑚𝑜𝑑𝑒𝑙 determines the precision of the model 𝑃 𝑚𝑜𝑑𝑒𝑙 by:

𝑃 𝑚𝑜𝑑𝑒𝑙 = 𝑡 1-𝛼 2 ⁄ ,(𝑁 -𝐶) 𝑠𝑡𝑢𝑑𝑒𝑛𝑡
. √𝜎 𝑚𝑜𝑑𝑒𝑙 Eq. 8

Data rejection

The tests have been performed at different time periods by two different operators, with membranes subject to degradation. To perform an efficient statistical analysis, systematic mass balances (total gas, methane and carbon dioxide) have tested the result consistency. The mass balance errors were calculated according to:

𝜀 𝑡𝑜𝑡 =
𝑄 𝑜𝑢𝑡 + 𝑄 𝑜𝑓𝑓 -𝑄 𝑖𝑛 𝑄 𝑖𝑛 Eq. 9

𝜀 𝑖 = 𝑦 𝑖 𝑜𝑢𝑡 𝑄 𝑜𝑢𝑡 + 𝑦 𝑖 𝑜𝑓𝑓 𝑄 𝑜𝑓𝑓 -𝑦 𝑖 𝑖𝑛 𝑄 𝑖𝑛 𝑦 𝑖 𝑖𝑛 𝑄 𝑖𝑛 Eq. 10
Most of the experiments display a very good mass balance (𝜀 < 5%) for the three indicators (Figure II). In such a process, a 15% error can even be considered admissible. As a general trend, the CH4 and CO2 errors tend to be positive and negative respectively. At the end of the campaign, a batch of tests has a stronger CO2 error (10 -25 %). These tests are performed at the extremity of the parameter axis. Three of these tests had replicates with a better mass balance. Therefore, no experiment with a mass balance error over 15% was integrated in the statistical analysis. The observed effects of the process parameters displayed in Figure V are in accordance with the preliminary discussion in section 3.1. Only the pressure 𝑃 𝑣𝑎𝑐 applied at the desorption step is not found to be significant in the domain (100 -300 𝑚𝑏𝑎𝑟 𝑎). The inlet biogas quality is strongly significant both in the first order, the second order and by its interactions with other parameter (𝑄 𝑙 And 𝑃 𝑔 ). This finding means that the process performances is highly dependent on the inlet biogas quality 𝑦 𝐶𝐻 4

𝑖𝑛 . The flowrates have limited quadratic effects:

the gas flowrate quadratic effect 𝑄 𝑔 is almost not significant (𝑎 ¨𝑃5 2 = 0.39). This other important result indicates that both parameters can be aggregated into the 𝑄 𝑔 𝑄 𝑙 ⁄ ratio for next experimental campaigns. In order to optimize the biomethane quality, 𝑃 𝑔 must be set high while the 𝑄 𝑔 𝑄 𝑙 ⁄ ratio must be low. From the DoE findings, process parameters were varied to produce the targeted biomethane quality in a liquid closed-loop configuration. For a given gas flowrate (𝑄 𝑔 = 153 NL. ℎ -1 ), a 𝑄 𝑔 𝑄 𝑙 ⁄ ratio of 1.30 𝑁𝐿. 𝐿 -1 produced biomethane (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 97.6%) with a methane recovery rate 𝑅 𝐶𝐻 4 of 77.6%. A higher 𝑄 𝑔 𝑄 𝑙 ⁄ (1.69 𝑁𝐿. 𝐿 -1 ) was still fulfilling the quality requirement with an improved recovery rate (𝑅 𝐶𝐻 4 = 83.7%). The recovery rate are still low compared to the available technologies which allow to recover more than 97% of the inlet methane flow [START_REF] Yang | Progress and perspectives in converting biogas to transportation fuels[END_REF], [START_REF] Bauer | Biogas upgrading -Review of commercial technologies[END_REF]. Nonetheless, this experimental set-up produces on one side a methane-rich fraction which can be injected into the gas grid, and on the other side a CO2-rich fraction. Compared to other available technologies, this offgas flow is not diluted in an air flow (as for air stripping) and still contain a significant methane fraction (respectively 𝑦 𝐶𝐻 4 𝑜𝑓𝑓 = 25.3 𝑎𝑛𝑑 20.0%). This methane content is high enough to burn this offgas flow in a low calorific burner and cover the digester heat demand. This process point is thus a first proof-of-concept for a biogas upgrading process of industrial interest. In an ideal exchanger, a decrease in the water temperature improves the absorption performance. It is a way towards process intensification. Under 𝑇 𝑙 = 20°𝐶, a theoretical 𝑄 𝑔 𝑄 𝑙 ⁄ ratio of 5.29 𝑁𝐿. 𝐿 -1 is needed to obtain a 97% quality biomethane. The DoE model built from experimental data reaches 90% biomethane quality only for ratios below 1.10 𝑁𝐿. 𝐿 -1 . It is well below the performances of an ideal exchanger.

Two hypothesis are possible explanations to this observation. Firstly, the process was operating in a closedloop configuration. As a result, the inlet absorbent may not be fully regenerated before entering the absorption module. On the contrary, the ideal exchanger curves suppose that the inlet liquid is free of dissolved gases. The second assumption is that the experimental set-up used in the DoE campaign (i.e. a X-50 absorption module) is limited by mass transfer and therefore cannot meet the performances of an ideal exchanger. An increase of the exchange surface would improve the process performances. The second hypothesis is tested in the following paragraph.

Influence of the absorption surface

According to the transfer unit methods, a longer exchanger leads to an improved biomethane quality. Therefore, two membranes modules were set in series to double the exchange surface and test the mass transfer limitation hypothesis. An open-loop configuration was set to avoid the question of the desorption efficiency.

The biomethane quality increases by 10-15 points when the exchange surface is doubled (Figure VII) but conserves the dependency on the 𝑄 𝑔 𝑄 𝑙 ⁄ ratio. These results are compared to simulated results from a 1D model describing both CO2 and CH4 absorption. The constant overall mass transfer coefficient 𝐾 𝐿 is an adjustable parameter of the 1D model. The adjusted 𝐾 𝐿 is not constant between the two series but still in the same range (1.9 × 10 -5 𝑚. 𝑠 -1 and 2.5 × 10 -5 𝑚. 𝑠 -1 respectively for one and two membrane modules).

Figure VII: Comparison of biomethane quality at different 𝑄 𝑔 𝑄 𝑙 ⁄ ratio with one or two membrane modules Open-loop configuration

The adjusted values of the model are consistent with the overall mass transfer coefficient derived from the calculation of the average mass transfer potential by Eq. 5 (Table VI). The overall mass transfer coefficient relative to the gas phase is very slightly below literature values (2.4 -4.0 10 -5 𝑚. 𝑠 -1 as reported by Lu [START_REF] Lu | Membrane contactor for CO2 absorption applying amino-acid salt solutions[END_REF]). 𝐾 𝐿 is constant for the 1-module serie (2.02 ± 0.07 10 -5 𝑚. 𝑠 -1 ). It slightly increases during the 2module serie and a higher standard deviation is observed (2.80 ± 0.19 10 -5 𝑚. 𝑠 -1 ).

The increase of the membrane area could improve the biomethane composition. Though, simulations from the 1D model indicates that 7 membrane modules in series . This ratio was chosen to be slightly below the ratio needed in the case of an ideal mass exchanger (Figure VI). Moreover, the inlet liquid is supposed to be free from dissolved gas content.

The important number of modules is explained by a reduced average absorbed CO2 flux 𝐽 𝐶𝑂 2 . Figure VIII shows that this flux is reduced by 20 % from the 1module serie to the 2-module serie. In the second module, CO2 has already been largely removed from the gas flow and therefore the mass transfer potential is lower than in the 1-module serie. In the membrane-wet regime, similar fluxes were reported for the absorption of pure CO2 under atmospheric pressure in water (0.9 -1.4 10 -3 𝑚𝑜𝑙. 𝑚 2 . 𝑠 -1 ) [START_REF] Mavroudi | A study of mass transfer resistance in membrane gas-liquid contacting processes[END_REF]. The methane recovery rate 𝑅 𝐶𝐻 4 was calculated by Eq. 2 assuming that the recycled gas flow would be mixed with process inlet biogas flow. Indeed, the methane content in the recycled gas flow would not significantly affect the inlet biogas quality (𝑦 𝐶𝐻 4 𝑖𝑛 = 60%). 𝑅 𝐶𝐻 4 thus decreased from 93.5 to 88.5% with the increase of 𝑃 𝑖𝑛𝑡 . As 𝑃 𝑖𝑛𝑡 increases, the recycled gas flow

Figure I :

 I Figure I: Process flow diagram for the reference process architecture Gas composition are analysed at the inlet, biomethane outlet and offgas outlet with an Agilent 490 micro Gas Chromatography (micro-GC) equipped with a thermal conductivity detector (TCD) detector. 𝑂 2 , 𝑁 2 and 𝐶𝐻 4 are separated on a 10m Molsieve 5Å column while 𝐶𝑂 2 is separated on 10m PoraPLOT U. The micro-GC equipment was combined with a flow-through selection valve (VICI 6-streams selector valve). The pressure of the biogas inlet and biomethane outlet samples are regulated below 1 bar g to comply with micro-GC specifications. 4 standard gases are used for calibration. Due to the combination of the error from the gas sampling system, the calibation error, the analytical repeatability and the analytical precision, the analytical error is estimated to ±0.3% 𝑣𝑜𝑙.. The liquid flow was operated in a closed-loop configuration during DoE. Two pumps are respectively regulating the pressure and the flowrate. They are controlled by a pressure transmitter PT 11 and a vortex flowmeter FI (Liqui-View) through frequency converters. Liquid pH is measured in a flow-through fitting positioned after CO2 absorption step (Mettler-Toledo InPro 4262i).The DoE investigated the available operating conditions range presented in TableIII.

Figure II :

 II Figure II: Mass balance errors on the DoE experiments 3.5 ResultsThe residual distribution of the quadratic model with interactions fits a normal density function (FigureIII): the error is randomly distributed within the domain. Moreover, no deviation of these residuals were observed neither with the experiment order nor with the measured response. The analysis of the model response in FigureIVshows that the model is in agreement with the measured response in a 95% confidence interval. The precision is though limited to 𝑃 𝑚𝑜𝑑𝑒𝑙 = ±4.0% (Eq. 8).

Figure III :

 III Figure III: Residuals distribution for the methane content response (quadratic model with interactions)

Figure V :

 V Figure V: First order effects and significant second order parameters from the quadratic model

4 STRATEGIES FOR PROCESS IMPROVEMENT 4 . 1

 41 Limitation in the DoE experimentsThe DoE confirmed the gas-to-liquid ratio to be a significant parameter to adjust the process performances. To further investigate the role of this aggregated parameter, the response of the DoE model was plotted against the performances of an ideal mass transfer exchanger (liquid outlet at equilibrium with inlet biogas) in Figure VI. The inlet biogas composition was set to a typical biogas composition (𝑦 𝐶𝐻 4 𝑖𝑛 = 60%) and the absorption pressure was maximised to 𝑃 𝑔 = 5 𝑏𝑎𝑟 𝑔 in order to improve the biomethane quality.

Figure VI :

 VI Figure VI: Influence of 𝑄 𝑔 𝑄 𝑙 ⁄ ratio on the biomethane quality: comparison of the DoE quadratic model with ideal exchanger performances DoE model is applied with 𝑄 𝑔 = 205 𝑁𝐿. ℎ -1

Figure VIII :

 VIII Figure VIII: Absorbed CO2 flux at different 𝑄 𝑔 𝑄 𝑙 ⁄ ratio with one or two membrane modules Open-loop configuration 4.3 Addition of a methane recycling loop The impact of the 𝑄 𝑔 𝑄 𝑙 ⁄ ratio was highlighted both from DoE findings and from Figure VII. The biomethane quality is improved with a lower 𝑄 𝑔 𝑄 𝑙 ⁄ ratio: it reaches 𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 94.6% under 𝑄 𝑔 𝑄 𝑙 = 1.00 𝑁𝐿. 𝐿 -1 ⁄ . In the meantime, the methane recovery rate decreases (𝑅 𝐶𝐻 4 = 84.2% under 𝑄 𝑔 𝑄 𝑙 = 1.00 𝑁𝐿. 𝐿 -1 ⁄). A new process arrangement is suggested and assessed to minimize this methane slip.Commercial biogas water scrubbers are designed with an intermediate step for the regeneration of the CO2-rich solvent. A flash-tank is set to decrease the pressure at an intermediate level in order to recycle part of the dissolved methane content[START_REF] Xiao | CO2 Removal from Biogas by Water Washing System[END_REF]. A new membrane contactor was therefore added in the liquid closed-loop to mimic this architecture (FigureIX).

Figure IX :

 IX Figure IX: Simplified closed-loop configuration with methane recycling loop The process performances were compared against the results obtained in the reference architecture (Figure I) under the same operating conditions. The recycled flow was controlled by a needle valve to investigate the influence of the gas pressure in the intermediate membrane contactor. The biomethane quality was maintained constant (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 95.5 -95.9%) and comparable to the methane content obtained in the reference case (Figure X). By progressively closing the needle valve, the intermediate gas pressure 𝑃 𝑖𝑛𝑡 increased from 0.25 to 1.9 𝑏𝑎𝑟 𝑔. It increased the methane content in the recycled gas flow 𝑦 𝐶𝐻 4 𝑟𝑒𝑐 from 57.2 to 76.9 %.

Figure X :

 X Figure X: Addition of methane recycling loop: influence of intermediate gas pressure

Table I :

 I Geometric features of the 2.5"x8" Extra-Flow module

	Parameter	Value	Source
	External membrane surface	1.4 𝑚 2	[34]
	the fiber bundle 𝑅 𝑖 Inner radius of	× 10 -2 𝑚 1.10	[34], [35]
	the fiber bundle 𝑅 𝑜 Outer radius of	× 10 -2 𝑚 2.32	[34]
	Fiber length	0.203 𝑚	[34], [35]
	Fiber effective length 𝑙 𝑒𝑓𝑓	0.146 𝑚	[34]
	Number of fibers 𝑁 𝑓	10 200	[35]
	Specific area 𝑎	𝑚 2 /𝑚 3 4 089	Calculated
	2 EXPERIMENTAL APPARATUS
	2.1 Hollow fiber membranes and modules

Table II :

 II Hollow fiber membrane data from manufacturer

	Parameter	X-40	X-50
	Applications	Degassing of	CO2
		other gases	degassing
	Material	PP	PP
	Inner fiber diameter	200 µ𝑚	220 µ𝑚
	External fiber	300 µ𝑚	300 µ𝑚
	diameter
	Porosity	20 -25%	40 -45%
	Tortuosity	2 -3	2 -3
	Average pore	0.03 µ𝑚	0.03 µ𝑚
	diameter
	Maximum liquid	7.2 𝑏𝑎𝑟 𝑔	7.2 𝑏𝑎𝑟 𝑔
	operating pressure
	(for 𝑇 < 40 °𝐶)

Table III :

 III Range of the operating conditions

	Operating parameter	Range	Unit
	𝐏 𝐠	0 -7	𝑏𝑎𝑟 𝑔
	𝐐 𝐂𝐎 𝟐	74 -739	𝑁𝐿. ℎ -1
	𝐐 𝐂𝐇 𝟒	76 -530	𝑁𝐿. ℎ -1
	𝐏 𝐥 Water flowrate 𝑸 𝒍	0 -8 30 -240	𝑏𝑎𝑟 𝑔 𝐿. ℎ -1
	2.3 Data treatment
	Experimental results are compared according to two
	key performance indicators: the biomethane quality 𝑦 𝐶𝐻 4

Table IV :

 IV Expected influence of process parameters 𝑃 𝑔 and the liquid temperature T l impact the gas-liquid equilibrium (i.e. the partition coefficient m). According to Henry's law, the liquid pressure 𝑃 l has no impact on this equilibrium. Its impact on the gas absorption process is found through the increase of transmembrane pressure. Indeed, membrane wetting is a limitation to mass transfer which is enhanced by a high transmembrane pressure. A low vacuum level 𝑃 𝑣𝑎𝑐 applied at the desorption step increases the mass transfer potential. Therefore, a more efficient solvent regeneration is expected.The inlet gas composition 𝑦 𝐶𝐻 4

	Parameter	Influence	𝐲 𝐂𝐇 𝟒 𝐨𝐮𝐭	𝐑 𝐂𝐇𝟒
	𝑸 𝒈	Increase the CO2 flow to be removed	-	+
	𝑸 𝒍	Increase the potential gas removal both for CO2 and CH4	+	-
	𝑷 𝒈	Decrease the partition coefficient m, thus the liquid potential gas removal is increased	+	-
	𝑷 𝒍	Increase the membrane wetting and consequently the mass transfer resistance [37]	-	+
	𝑷 𝒗𝒂𝒄	The lower the pressure at desorption, the higher the mass transfer potential to	-	+
		regenerate the absorbent
	𝑻 𝒈	Influence of gas temperature is not known but likely not to be significant	?	?
	𝑻 𝒍	A lower temperature enhances the absorption capacity [36]	+	-
	𝒚 𝑪𝑯 𝟒 𝒊𝒏	Decrease the CO2 mass transfer potential (𝑥 * -𝑥) Lower the CO2 flow to be removed	+	-
	at absorption

Table V :

 V Definition of parameters' levels for the Design of Experiments

	Level	𝑸 𝒍 [𝑳. 𝒉 -𝟏 ]	𝑸 𝒈 [𝑵𝑳. 𝒉 -𝟏 ]	𝑷 𝒈 [𝒃𝒂𝒓 𝒈]	𝑷 𝒗𝒂𝒄 [𝒎𝒃𝒂𝒓 𝒂]	𝒊𝒏 𝒚 𝑪𝑯 𝟒 [%𝒗𝒐𝒍. ]
	Parameter	P1	P4	P3	P2	P5
	-θ	60	205.0	1.10	100.0	30.0
	-1	103.5	355.7	2.2	157.9	47.4
	0	135	465.0	3.0	200.0	60.0
	+1	166.5	574.3	3.8	242.1	72.6
	+θ	210	725.0	4.90	300.0	90.0

Table VI :

 VI Overall mass transfer coefficients calculated from experimental data

	N modules	𝑸 𝒍 [𝑳. 𝒉 -𝟏 ]	𝑸 𝒈 /𝑸 𝒍 [𝑵𝑳. 𝑳 -𝟏 ]	𝑷 𝒈 [𝒃𝒂𝒓 𝒈]	𝒚 𝑪𝑯 𝟒 𝒊𝒏 [%𝒗𝒐𝒍. ]	𝐊 𝑳 [𝟏𝟎 -𝟓 𝒎. 𝒔 -𝟏 ]	𝐊 𝑮 [𝟏𝟎 -𝟓 𝒎. 𝒔 -𝟏 ]
	1	210	0.7	5.0	60.1	2.02	1.85
		210	1.0	5.0	60.0	2.02	1.85
		210	1.7	5.0	59.8	1.96	1.80
		210	2.5	5.0	59.9	1.94	1.78
		210	3.2	5.0	60.0	1.99	1.80
		210	4.0	5.0	60.0	2.15	1.92
	2	210	2.5	5.0	59.8	2.58	2.38
		210	3.2	5.0	60.3	3.03	2.73
		210	4.0	5.0	60.0	2.81	2.59
	are required to reach the targeted biomethane quality
	(𝑦 𝐶𝐻 4 𝑜𝑢𝑡 > 97%) under 𝑄 𝑔 𝑄 𝑙 ⁄ = 4.0 𝑁𝐿. 𝐿 -1

ACKNOWLEDGEMENTS



This work was financially supported by the Marne departmental council, the Champagne-Ardenne regional council and Reims Metropole through the Biotechnology Chair of CentraleSupelec.

4.4 On-site experimental campaign with farm biogas These results obtained with synthetic biogas had to be confirmed during an on-site experimental campaign to upgrade real raw biogas. The farm AD unit co-digests cow slurry with agroindustrial waste (beet pulp, brewing dregs…), crop residues or grass silage. A by-pass flow is diverted from the average 30 𝑁𝑚 𝑏𝑖𝑜𝑔𝑎𝑠

3

. ℎ -1 flow to be upgraded in the pilot process. During the 2-week field campaign, the methane content was stable (𝑦 𝐶𝐻 4 𝑖𝑛 = 56.3 -59.0%). An air injection inside the digester maintains the 𝐻 2 𝑆 content low (< 200 𝑝𝑝𝑚𝑣) by sulfur oxidation. Though this air injection results in the presence of nitrogen and oxygen in the biogas flow (1.7 -2%). This air injection was stopped during 3 days for the experimental campaign purpose.

From the previous values of 𝐾 𝐿 , this experimental campaign was operated with a higher inlet biogas flowrate (𝑄 𝑔 = 229 NL. ℎ -1 ) to intensify the pilot process to its maximum. In addition, due to biogas compressor limitations, the absorption pressure 𝑃 𝑔 was set only to 4.5 𝑏𝑎𝑟 𝑔. The 97% biomethane quality was almost reached with a 𝑄 𝑔 𝑄 𝑙 ⁄ ratio of 1.09 𝑁𝐿. 𝐿 -1 (Table VII). The same experimental process point was repeated with a lower inlet air content. The biomethane quality was thus increased from 94.3 to 96.2 %.

Similar biogas upgrading performances were achieved in the reference architecture both with synthetic and raw farm biogas. The process displayed similar average mass transfer coefficients with synthetic biogas and real biogas (2.62 -2.85 × 10 -5 𝑚. 𝑠 -1 against 2.92 -3.07 × 10 -5 𝑚. 𝑠 -1 ). The slight increase may be due to an increased liquid flowrate during farm experiments. Therefore, it is reasonable to assume that the process optimization can be carried out using synthetic biogas at the lab for the sake of convenience.

Moreover, within the 2 weeks, no significant performance degradation was observed. Therefore, the exposure of the membrane material is not extremely sensitive to the presence of biogas impurities (like 𝐻 2 𝑆 or volatile organic compounds). The polymer compatibility with these impurities must be carefully study for future industrial developments.

CONCLUSION

The use of membrane contactors for biogas upgrading purposes was optimized in an absorption/desorption cycle with a liquid closed-loop. The reference process configuration with equivalent membrane areas for the absorption and desorption steps was optimized by a Design of Experiments. A high biomethane purity (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 97.6%) was obtained with an improved methane recovery rate compared to previous reported results (𝑅 𝐶𝐻 4 = 83.7% against 75%). This moderate recovery rate is a first process point of interest as the calorific content of the nondiluted offgas flow (𝑦 𝐶𝐻 4 𝑜𝑓𝑓 = 20.0 -25.3%) may be used to cover the digester heating demand in a low calorific burner. Doubling the membrane area at the absorption step proved that a significant number of modules [START_REF] Rotunno | Energy and economic analysis of a water scrubbing based biogas upgrading process for biomethane injection into the gas grid or use as transportation fuel[END_REF] are needed to approach an ideal mass exchanger under conditions 𝑄 𝑔 𝑄 𝑙 ⁄ = 4.0 𝑁𝐿. 𝐿 -1 close to the theoretically optimized flowrates ratio. An alternative process arrangement mimicking the conventional WS methane recycling loop revealed a significant improvement of the methane recovery rate without affecting the biomethane quality (𝑦 𝐶𝐻 4 𝑜𝑢𝑡 = 95.9% / 𝑅 𝐶𝐻 4 = 93.5%). This second process point is a very promising result to combine the advantages of membrane contactors and WS into an innovative biogas upgrading process. Additional studies are required to optimize the membrane area ratio between the absorption, partial and final degassing steps (set to 1:1:1 in the present study) and to meet the performances of available biogas upgrading technologies.