Ces trois années de thèse menées entre le CEA de Marcoule et l'IMFT ont été passionnantes. Je remercie encore mes encadrants, Sophie Charton, Eric Climent et Dominique Legendre, pour m'avoir donné cette opportunité et m'avoir accompagné au quotidien dans cette thèse. Merci la team.

Je tiens à remercier les collègues et amis que j'ai côtoyés au CEA et à l'IMFT. Pour leur soutien et leur aide tout au long de mes recherches. Leurs amitiés, leurs encouragements et leurs conseils ont été très enrichissants. Je garderai d'excellents souvenirs des rendezvous "échec et goûter" du pigeonnier ainsi que les "cafés Camichel" et petit-déj Doct'eau à l'IMFT. J'ai eu la chance de participer à quatre conférences dans cette thèse. Ce qui m'a permis d'avoir de nombreuses discussions scientifiques avec les meilleurs chercheurs et doctorants de mon domaine (notamment sur les mailleurs à Kobé). Ces moments sont finalement assez rares et privilégiés. Ainsi, je remercie le LRVE de m'avoir fait confiance et de m'avoir donné ces opportunités.

LIST OF TABLES

the favorable effect of this bifurcation towards 3D flows on mass transfer, by promoting convective transport of the solute in and around the drop.

Thanks to direct numerical simulations, we have described a wide range of uniform or oscillating flows, over a range of Reynolds (likely to be encountered in extraction processes) and viscosity ratios (bubbles, drops, solid particles), and translated their impact on mass transfer, drag coefficient and history force. This has led to correlations directly applicable to processes simulation.

CONTENTS

List of Figures

More instants are shown during the second half of the period to provide a detailed description of the two recirculation vortices that develop both inside and outside the droplet during the deceleration phase. . . . . . . 2.3 (a) Evolution of the normalized tangential velocity U θ /W 0 with the polar angle θ at the interface. (b) Radial evolution of the normalized axial velocity U x (x = 0, r)/W 0 . The times considered are the same as in Fig. 2.2: t = 0T , W (t) = W 0 (Fig. 2.2a), t = 0.554T , W (t) = 0.667W 0 (Fig. 2.2b), t = 0.582T , W (t) = 0.507W 0 (Fig. 2.2c), t = 0.670T , W (t) = 0.124W 0 (Fig. 2.2d), t = 0.800T , W (t) = 0.0489W 0 (Fig. 2.2e), t = 0.883T , W (t) = 0.329W 0 (Fig. Comparison of simulation results (point) and correlations, Hadamard [START_REF] Hadamard | Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux[END_REF] and Rybczynski et al. [START_REF] Rybczynski | On the translatory motion of a fluid sphere in a viscous medium[END_REF] with Mei et al. [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF] and Schiller and Naumann [START_REF] Schiller | Über die grundlegenden berechnungen bei der schwerkraftaufbereitung[END_REF] [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], Parentheses: Feng and Michaelides [START_REF] Zhi | Drag coefficients of viscous spheres at intermediate and high reynolds numbers[END_REF] (parenthesis), and Brackets: Oliver and Chung [START_REF] Oliver | Flow about a fluid sphere at low to moderate reynolds numbers[END_REF]. The differences between our results and those of Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF] obtained with the same code, can be explained by the differences in mesh parameters and steady-state criteria. . . . . . . . . 3.4 2D and 3D drag coefficients, and amplification factor (see Sec. 3.6.3).

List of Tables

Comparison of results obtained from simulations achieved with the "full polar" and the hybrid "polar and streamlines" meshes. . . . . . . . . . Une précédente thèse, Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], a permis de mettre en lumière les limitations causées par les hypothèses simplificatrices des modèles de transfert de matière de l'état de l'art, classiquement utilisés en génie des procédés. Nos travaux consistent à comprendre les sources de ces erreurs, pouvant être significatives, afin de pouvoir améliorer la précision de ces modèles, à partir de considérations physiques. De nouvelles corrélations ont ainsi été proposées dans le but de reproduire fidèlement les résultats de simulations numériques directes, dans une large plage de régime hydrodynamique et de conditions physico-chimiques. La démarche a également été étendue au cas du transfert réactif.

L'étude hydrodynamique ayant conduit à ces nouvelles corrélations de transfert de matière a également permis de progresser dans la description et la modélisation des différentes forces subies par la goutte en écoulement.

Abstract

The chemical industry extensively employs liquid-liquid extraction techniques, and they are now attracting renewed interest in the hydro-metallurgical sector as a means of establishing a recycling sector. They are also at the heart of the treatment-recycling of irradiated nuclear fuels. Faced with these new challenges, numerical simulation has a major role to play in the development or adaptation of processes. In this context, it is essential to have predictive models representing the impact of all the physical phenomena involved, at different scales. This thesis is devoted to the study of couplings between hydrodynamics and mass transfer kinetics at the level of a spherical droplet in a uniform flow of an immiscible phase.

A previous thesis, Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], highlighted the limitations caused by the simplifying assumptions of state-of-the-art mass transfer models, classically used in process engineering. The aim of our work is to understand the sources of these errors, which can be significant, to improve the accuracy of these models, based on physical considerations. New correlations have been proposed with the purpose of faithfully reproducing the results of direct numerical simulations, over a wide range of hydrodynamic regimes and physicochemical conditions. The approach has also been extended to the case of reactive transfer.

The hydrodynamic study that led to these new mass transfer correlations also enabled progress to be made in the description and modeling of the various forces experienced by the drop. These forces have a decisive influence on the residence time of droplets in liquid-liquid contactors, during which the two phases are likely to exchange. A review of the literature led us to identify two terms that are still poorly described in the literature: the history force and the steady drag in 3D bifurcated flow.

The study of the history force is very complex in the case of drops, due to the coupling between internal and external flows. We have proposed and validated a model that takes advantage of analytical solutions for a partial slip-type boundary condition, and is therefore applicable to bubbles and solid particles, to predict the evolution of the history force for the case of drops. Finally, we studied the 3D flow bifurcation in and outside a spherical drop. Thanks to an extensive parametric study, corroborated by original experiments, we have identified a bifurcation criterion and modeled the impact of the 2D/3D transition on the stationary drag coefficient. We have also highlighted

An introduction to PUREX process.

The PUREX process consists of an initial stage of putting nuclear combustible waste composed of uranium and plutonium, as well as the other actinides and most of the fission products, into a nitric solution. Then, uranium and plutonium are separated and purified, allowing the recycling of plutonium as a mixed uranium-plutonium oxide nuclear fuel, MOX. The separation is achieved by liquid-liquid extraction techniques using TBP (Tri-Butyl Phosphate) diluted in an organic solvent: TPH (Tetra-Propylene). Associated with these main operations are operations to recycle the nitric acid and the solvent (TBP and TPH), and to treat liquid and gaseous effluents and solid waste.

The operations of extraction of uranium and plutonium by TBP constitute the heart of the PUREX process. In the nitric solution of dissolution, the elements are found in the following forms:

• uranium is present in the + VI oxidation state in the form of the uranyl ion UO + 2

• plutonium is present essentially in the +IV oxidation state

• the minor actinides, neptunium, americium, and curium, are present at different oxidation states: +V and +V I for Np, +III for Am and Cm

• the fission products of the lanthanide family are in the +III oxidation state, the alkalis (Cs, Rb) are in the +I oxidation state, the alkaline earths (Sr, Ba) are in the +II oxidation state, zirconium is in the +IV oxidation state, molybdenum, and technetium are in the anionic form in the +V I and +V II oxidation states respectively.

TBP is generally found to have a much higher affinity for elements in the oxidation state +IV and +V I than for elements in other states. Redox reactions hence play a major role in the PUREX process. The extraction mechanism is of the solvation type of uncharged complex:

M n+ + nA + xT BP ⇔ MA n ,xTBP (1.1)
Today, active research is dedicated to treatment and multirecycling of MOX. In this aim, new solvents, allowing simplification of the processes, i.e. no Redox, are under study. More generally, there is a crucial need for the recycling industry to develop solvent extraction processes with new solvents with increased selectivity and less environmental impact. TBP succession candidates are generally more viscous and therefore exhibit different hydrodynamic and mass kinetic behaviors.

IMPLEMENTATION AND MODELLING OF LIQUID-LIQUID EXTRACTION27

Implementation and modelling of liquid-liquid extraction

For nuclear application, the contactor, where the liquid-liquid extraction takes place, is generally a pulsed column, because of its low maintenance need and efficiency. The operation principle of a pulsed column consists of a circulating solvent (the light phase), moving from bottom to top in a vertical cylinder against the flow of the immiscible aqueous (heavy) phase.

Two modes of operation are possible: if the column is initially filled with solvent, the aqueous phase disperses and falls in the form of droplets that coalesce at the bottom of the column: this is the water in oil (W/O) mode. Conversely, if the column is filled with aqueous phase, it is the organic phase that disperses and rises in the form of droplets that coalesce at the top of the column: this is the oil in water (O/W) mode.

To produce and maintain the emulsion in the column, a back and forth movement of the liquids is produced via a pulsation leg connected to the column at the bottom. This movement is superimposed on the flow of the phases and creates turbulence at obstacles placed in the column. These obstacles are usually made of perforated discs and are usually referred to as "doughnuts" and "disks". Mass transfer takes place at the aqueous-organic interface of the droplet during liquid-liquid extraction.

A dedicated code, PAREX, has been developed at CEA to design and simulate complex solvent extraction processes such as PUREX. Like most process simulators, the model used in the PAREX code assumes mass transfer to occur thanks to a three steps transfer, according to the double film theory where each phase is divided into two distinct regions (Fig. 1.1):

• the core (or "bulk") of the phase where the concentration is uniform.

• a thin layer around the interface, called the "diffusion boundary layer". In this zone, the transport of the solute is assumed to be governed by diffusion, at steady state (Fick's law). The thickness of these boundary layers is fixed by convection effects.

• the interface, where thermodynamic equilibrium is assumed to prevail. The main phenomena encountered are introduced here, to guide the reader to the subject of this thesis and clearly define its scope in the global framework of current liquid-liquid extraction research.

FOCUS ON ISOLATED SPHERICAL DROPLET

29

A significant enhancement in the predictive performance of solvent extraction models cannot be achieved without a comprehensive description of the flow in the column and its interactions with the droplet population.

One of the keys is to be able to predict the amount of interfacial area in the column, regardless of the solvent used, the column geometry and size, and the operating conditions. Interfacial area depends on both the droplet size distribution (DSD) and on the holdup of the dispersed phase.

• Coupling a population balance equation (PBE) with a proper description of the turbulent flow is today the preferred method to predict the DSD in complex apparatus such as the pulsed column Amokrane et al. [START_REF] Amokrane | Single-phase flow in a pulsed column: Particle Image Velocimetry validation of a CFD based model[END_REF], Castellano [START_REF] Castellano | Multiscale study and modeling of dispersion properties relevant for liquid-liquid extraction : adaptation of breakup and coalescence kernels to industrial processes[END_REF].

• The prediction of the holdup is related to the DSD, and to the relative (or "slip") velocity of each class of droplets, that fix their residence time in the apparatus. This requires particular insight on the forces acting on the droplets, which is one major issue focused in this study.

Focus on isolated spherical droplet

Here, we focus on the simple system that is an isolated spherical droplet in a liquid-liquid extraction column, and especially on three features: its shape, velocity, and composition.

Regarding the shape, droplets moving in an immiscible liquid have been characterized by three different shapes, in Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF]:

• Spherical: if interfacial tension and/or viscous forces are much more important than inertia forces, the droplet takes the shape of a sphere.

• Ellipsoidal: Oblate droplet with a convex interface, generally the small semi axis is parallel to the flow direction, but it has been observed that oblate droplets can have an oscillation regime, making their identification difficult.

• Spherical-cap or ellipsoidal-cap: Large droplets tend to have flat or indented bases and lack any kind of symmetry. If the liquid particle has a concavity at the back, it's called "dimpled". Large spherical or ellipsoidal caps may also have a trail of dispersed fluid, referred to as "skirts".

In their book, Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF] proposed a widely adopted graphical correlation relating the Eötvös number (1.4), Morton number (1.2) and Reynold number (1.3) to droplet shape: The behavior of ellipsoidal droplets is well predicted by Grace correlation Grace et al. [START_REF] Jr Grace | Shapes and velocities of single drops and bubbles moving freely through immiscible liquids[END_REF] and has also been studied in Bozzano and Dente [START_REF] Bozzano | Shape and terminal velocity of single bubble motion: a novel approach[END_REF] or Gibbons et al. [START_REF] Joseph H Gibbons | Effect of a surface active agent on the velocity of rise of benzene drops in water[END_REF]. For cap's case, shape study can be found in Davenport et al. [START_REF] Davenport | Spherical cap bubbles in low density liquids[END_REF], Wairegi and Grace [START_REF] Wairegi | The behaviour of large drops in immiscible liquids[END_REF] or Collins [START_REF] Collins | A second approximation for the velocity of a large gas bubble rising in an infinite liquid[END_REF]; Terminal velocity study can be found in Davies and Taylor [START_REF] Davies | The mechanics of large bubbles rising through extended liquids and through liquids in tubes[END_REF] or Wu et al. [START_REF] Bjc Wu | Rise speed of spherical cap bubbles at intermediate reynolds number[END_REF]. In this thesis, we will mainly consider spherical droplet, as it is the most encountered shape in liquid-liquid extraction columns.

M o = gµ 4 e ∆ρ ρ 2 e σ 3 (1.2) 30CHAPTER 1. INTRODUCTION Re = 2ρ e U r µ e (1.3) Eo = 4g∆ρr 2 σ (1.
The droplet velocity results from the balance of the forces acting on it (second Newton law). We can decompose the force experienced by the droplet as a sum of different phenomena contribution, this is written in Eq. 1.5.

1.3. FOCUS ON ISOLATED SPHERICAL DROPLET m - → a = ----→ F weight + -----→ F buoyancy + ---→ F drag + -------→ F added-mass + ----→ F history (1.5)
With:

• ----→ F weight : Weight is equal to the fluid density time the droplet volume time the gravity, its direction points downward.

• -----→ F buoyancy : Buoyancy is equal to the difference between Archimedes and weight force because the internal fluid density in liquid-liquid extraction process is generally lower than the external one, its direction generally points upward.

• ---→ F drag : Steady drag force represents the force that opposes the motion of a droplet in a liquid under steady (no acceleration) condition, and therefore its direction is opposed to the droplet velocity.

• -------→ F added-mass : Added mass force is caused by the inertia of the continuous phase liquid that is displaced by the droplet during its acceleration and is therefore opposed to the droplet acceleration.

• ----→ F history : History force is the consequence of unsteady viscous effects when a droplet is accelerating and its direction is opposed to the droplet acceleration.

Weight and buoyancy forces are an easy task to compute, but for the three others they depend on the flow inside and outside the droplet, and their coupling at the interface, thus computing them can be more complicated. The fluid flow inside a spherical droplet is often described as a single toroidal vortex called a "Hill vortex" Hill [START_REF] Hill | On a spherical vortex[END_REF], shown in Fig. 1.4. This configuration is encountered over a wide range of parameters in the case of non-deformable spherical drops. The Hill vortex is a Stokes flow axisymmetric solution: independently of the cross-sectional plane along the x-axis of the droplet (the axis in the direction of flow and passing through its center). We observe two regions on either side of this axis comprising nested elliptical streamlines with a semi-major axis parallel to the x-axis, twice as large as the semi-small axis.

Using this solution obtained under Stokes flow condition, many correlations have been derived, and used to compute steady drag, added mass and history forces. Some have extended the Hill's vortex assumption to application at moderate Reynolds number. For steady drag force, this is well described in Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], and summarized in C D = 24

Re

Re ≪ 1 (solid particle) Schiller and Naumann [START_REF] Schiller | Über die grundlegenden berechnungen bei der schwerkraftaufbereitung[END_REF] C D = 24 Re (1 + 0.15Re 0.687 ) Re < 800 (solid particle)

Mei and Adrian [START_REF] Mei | Flow past a sphere with an oscillation in the free-stream velocity and unsteady drag at finite reynolds number[END_REF] C D = 16 Re 1 + 8 Re + 0.5 (1 + 3.315Re -0.5 )

-1 0.1 ≤ Re ≤ 200 (buble)
C D = 2-µ * 2 C D (Re, 0) + 4µ * 6+µ * C D (Re, 2) 5 < Re < 1000; 0 ≤ µ * ≤ 2 C D = 4 µ * +2 C D (Re, 2) + µ * -2 µ * +2 C D (Re, ∞) 5 < Re < 1000; 2 ≤ µ * ≤ ∞ C D (Re, 0) = 48 Re 1 -2.21 √ Re + 2.14

Re

With: Regarding steady drag force at larger Reynolds number, we know from Edelmann et al. [START_REF] Edelmann | Numerical investigation of different modes of internal circulation in spherical drops: Fluid dynamics and mass/heat transfer[END_REF] and Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], that the axisymmetric assumption is no longer respected. Therefore, the correlations shown in Table 1.1 obtained numerically are no longer valid for some conditions.

C D (Re, 2) = 17.0Re 2/3 C D (Re, ∞) = 24 Re 1 + 1 6 Re 2/3
The complexity of the coupling of the flow between both phases makes the expression of the history force very difficult to obtain in a liquid-liquid system. Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] proposed a first expression and validated it under creeping Stokes flow conditions, but more studies are needed on this subject, especially to be able to represent pulsed flux conditions, where the flow is governed by the oscillatory forcing.

THESIS PLAN

Understanding 3D flow bifurcations and their consequences on the velocity of an isolated spherical droplet, and improving the modelling of the history force in a pulsed column, are the main contributions on hydrodynamics of this thesis.

As previously mentioned, mass transfer from a droplet to a surrounding immiscible liquid is commonly described by the double-film model. This model considers that the transfer occurs through a series of films of width δ with a constant concentration gradient inside this film, see Fig. 1.1. The resistances of both films should be considered to properly quantify the apparent flux of solute exchanged between the two phases (coupled problem). Only when the transfer from one phase is significantly fast compared to the other one the corresponding resistance can be neglected (isolated internal/ external problems).

The challenge for solvent extraction is to find proper and robust correlations for the mass transfer coefficient, regardless of the physical-chemical settings and flow conditions. The main studies on mass transfer from an isolate droplet are assuming an isolated problem governed by either the inner our the outer fluid phase. They have been reviewed and discussed recently by Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF]. It was evidenced that using those correlations to predict the behavior of a coupled problem with the mass transfer resistance additivity law can lead to important errors and is therefore not recommended. In this work, we will later discuss how to correct the state-of-the-art correlation for a coupled mass transfer configuration.

Additionally, mass transfer between two liquids is often accompanied by a chemical change, as e.g. the redox reactions in the PUREX process. The chemical reaction consumes a portion of the solute extracted and therefore modifies the driving force of mass transfer. It is common to physically model this phenomenon as a mass transfer enhancement factor. In this thesis, we will consider the enhancement provided by a first order reaction occurring in the continuous phase.

Thesis plan

The thesis document is written as a collection of scientific contributions where each chapter is a published, under review or soon submitted article. Therefore, they are all self-sustained and do not require a global bibliography or material and methods chapter.

It follows this plan:

Chapter 2 is a published article that discusses the contribution of the history force to the total force acting on a droplet embedded in a uniform oscillating flow as a simplified configuration for a pulsed column. The history force is a transient force often neglected in Lagrangian simulations, due to the difficulty of integrating it in time and the lack of an analytical expression for it. We used direct numerical simulations (DNS) to study the contribution of the history force to the total force acting on a spherical droplet in time periodic oscillating flow. The results of the DNS simulations indicated that the 34CHAPTER 1. INTRODUCTION contribution of the history force to the total force is significant for systems with high viscosity ratios, such as droplets in liquids, and medium flow oscillation frequencies. We also showed that the kernel expression proposed by [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] can be used to accurately model the history force. This work has been published in Physical Review Fluids Godé et al. [START_REF] Godé | Bassetboussinesq history force acting on a drop in an oscillatory flow[END_REF].

Chapter 3 is an article that will be further sent for review. This study investigated the 3D bifurcation phenomena that occur inside and outside a spherical droplet in a uniform steady flow. The study used both experimental and numerical approaches to evidence two types of bifurcations: internal and external. The bifurcations can be either steady or unsteady, and they affect the forces on the droplet. The study also proposes a method to predict the 3D drag based on the bifurcation criterion and a drag coefficient amplification factor. The experimental results confirmed the simulation results.

Chapter 4 is an article currently under review. Direct numerical simulations (DNS) were used to study mass transfer from an isolated droplet in a liquid-liquid solvent extraction process. We conducted a wide parametric study for an isolated droplet in various liquids using DNS. The results of simulations were used to develop new correlations for the Sherwood number which is a dimensionless number that characterizes the rate of mass transfer from the droplet to the surrounding fluid. The new correlations are more accurate than existing correlations and are applicable to a wider range of operating conditions. They are also applicable to heat transfer problems.

Chapter 5 is an article that will be soon finalized for submission. It discusses, using direct numerical simulation, the enhancement of mass transfer provided by a first order chemical reaction. The approach aims at proposing a comprehensive method for predicting the Sherwood number, based on chemical and physical properties relevant to the emergence of new processes.

At last, chapter 6 summarizes all the results and findings of the thesis, and provides insight into possible applications and outlooks for this work.

Chapter 2

The Basset- This chapter is a published article, see Godé et al. [START_REF] Godé | Bassetboussinesq history force acting on a drop in an oscillatory flow[END_REF]. It discusses the contribution of the history force to the total force acting on a droplet embedded in a uniform oscillating flow. The total force experienced by the droplet can be divided into four contributions: steady drag force, inertial or pressure gradient force, added-mass force, and Basset-Boussinesq history force. The history force is the force resulting from the diffusion of a flow state to another when the droplet undergoes unsteady motion. The effect of the history force is often neglected, due to numerical difficulties for its time integration, but also because no analytic expression is available in the time domain for droplets.

The contribution of the history force acting on a spherical droplet in an oscillatory flow was determined using direct numerical simulation (DNS). Variation of the viscosity ratio makes the analysis relevant to bubbles, droplets, and solid particles. By changing the flow oscillation frequency, we determined the range of physical parameters for which the contribution of the history force is significant.

Additionally, the relevance of the kernel expression recently proposed by Legendre et al. (Physical Review Fluids 4, 073603, 2019) is discussed and compared to theoretical predictions in the frequency domain. We demonstrate that this kernel can be used to model history effects for droplets, especially since the contribution of the history force to the total force can be important. The proposed history force expression provides a significant improvement for droplet trajectory prediction when the history force contribution cannot be neglected.

Introduction

Predicting droplet dynamics is of primary importance in many applications related to natural and industrial situations, where different flow conditions and liquid/liquid systems can be encountered. Predicting droplet dynamics is particularly challenging in strongly unsteady situations. In this study, we focus on the history force contribution to the force balance. Depending on the viscosity ratio, the physical configuration may vary from liquid droplets in a gas for cloud microphysics Ravichandran et al. [START_REF] Ravichandran | Fluid Dynamics in Clouds[END_REF], gas bubbles in a liquid as encountered in bubble columns for petrochemical, and chemical industries or liquid-liquid systems related to solvent extraction, petroleum, or other processes involving emulsions.

Under such operating conditions, the flow of the carrier fluid is commonly turbulent, which means that droplets experience a wide range of temporal fluctuations. Modelling approaches based on Euler-Lagrange formalism are well suited to account for the effects of temporal and spatial turbulent structures of the fluid flow on particle trajectories obtained through time integration of the force balance. Analytic formulations of the force balance are built by integrating specific forces: i) the drag force related to the resistance due to the relative (or slip) velocity between the droplet and the surrounding fluid, ii) the history force which is a correction to steady drag for fluctuating slip velocities, accounting for the unsteady viscous diffusion of vorticity produced at the interface, iii) the added mass force that originate from inertial effects due to both the droplet and fluid accelerations, iv) the inertial or pressure gradient force due to the mean flow oscillations, and v) the lift force induced by the presence of velocity gradients in the flow field whose contribution is orthogonal to the relative motion.

History force has often been neglected because its integration over viscous diffusion times of several order of magnitudes requires tedious algorithms (high computational cost and large memory storage of the relative acceleration). But also because its expression has been an open question for viscous fluid particles. However, recent literature has demonstrated the need to consider this force, as in Haller [START_REF] Haller | Solving the inertial particle equation with memory[END_REF], Daitche [START_REF] Daitche | On the role of the history force for inertial particles in turbulence[END_REF], Candelier et al. [START_REF] Candelier | On the effect of the boussinesq-basset force on the radial migration of a stokes particle in a vortex[END_REF]. The history force has been found to have a significant impact on the particle's slip velocity and acceleration, resulting in modified preferential accumulation and collision rates. It is now possible to consider the history force in CFD because of the continuous increase in computing resources, but also thanks to numerical algorithms dedicated to the temporal integration of the history kernel Daitche [START_REF] Daitche | Advection of inertial particles in the presence of the history force: Higher order numerical schemes[END_REF], Prasath et al. [START_REF] Prasath | Accurate solution method for the maxey-riley equation, and the effects of basset history[END_REF], Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF]. In addition, a history kernel for the history force has recently been proposed by Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF].

In this paper, we study the hydrodynamics of a fluid sphere in a uniform oscillatory flow. We aim at specifying with direct numerical simulations the conditions under which the history force is significant. The impulsive flow conditions studied in Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] are extended to the periodic flow encountered in applications e.g., in oscillatory flow reactors and pulsed columns.

The paper is organized as follows. First, we describe the flow configuration and the numerical methods. Second, we describe the flow structure, both inside and outside the droplet. Then, we discuss the quantitative contribution of the history force to the total momentum balance and conclude on the accuracy of the proposed analytic history force kernel over a wide range of the relevant dimensionless parameters, i.e. the viscosity ratio and the droplet Stokes number.

Problem statement

We consider a fluid sphere of radius R, density ρ i and viscosity µ i moving at a velocity V (t) in an unbounded continuous fluid phase of density ρ e , viscosity µ e and velocity U (t). We define W (t) = U (t) -V (t) as the relative velocity of the fluid with respect to the moving sphere.

Under Stokes flow conditions, the fluid sphere experiences a total force F (t) which breaks down into four contributions: the steady drag force F SD Hadamard [START_REF] Hadamard | Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux[END_REF], Rybczynski et al. [START_REF] Rybczynski | On the translatory motion of a fluid sphere in a viscous medium[END_REF], the added-mass force F AM Gatignol [START_REF] Gatignol | The faxen formula for a rigid particle in an unsteady non-uniform stokes flow[END_REF], Maxey and Riley [START_REF] Maxey | Equation of motion for a small rigid sphere in a nonuniform flow[END_REF], the inertial or pressure gradient force F I , and the Basset-Boussinesq history force F H described in Boussinesq [START_REF] Boussinesq | Sur la resistance quoppose un liquide indefini en repos, sans pesanteur, au movement varie d'une shere solide[END_REF], Basset [START_REF] Barnard | On the motion of a sphere in a viscous liquid[END_REF] for a solid sphere, in Mei et al. [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF], Yang and Leal [START_REF] Man | A note on memory-integral contributions to the force on an accelerating spherical drop at low Reynolds number[END_REF] for spherical bubbles, and recently revisited by Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] for fluid viscous spheres.

38CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW

These contributions write respectively from left to right:

F (t) = 6πµ e R 2 + 3µ * 3 + 3µ * W (t) + C m ρ e 4 3 πR 3 dW dt + ρ e 4 3 πR 3 dW dt +6πµ e R t 0 dW dt ′ K µ (t -t ′ , µ * ) dt ′ , ( 2.1)
where µ * = µ i /µ e is the viscosity ratio, C m is the added mass coefficient which is equal to 1 2 for spheres and K µ the history force kernel whose expression depends on µ * . In the following, the added mass force and the inertial (pressure gradient) force are combined into a single term.

Boussinesq [START_REF] Boussinesq | Sur la resistance quoppose un liquide indefini en repos, sans pesanteur, au movement varie d'une shere solide[END_REF] and Basset [START_REF] Barnard | On the motion of a sphere in a viscous liquid[END_REF] proposed an expression for K µ for solid particle in the time domain. This kernel is valid for systems with high µ * , i.e., when the viscosity of the fluid sphere is much higher than that of the surrounding fluid yielding no-slip at the surface:

K µ (t, µ * → ∞) = 1 πt/t ν , ( 2.2)
where t ν = R 2 /ν e is the characteristic diffusion time in the surrounding fluid. For a bubble (shear free interface), or when the viscosity of the fluid sphere is much lower than the surrounding fluid, the kernel writes Yang and Leal [START_REF] Man | A note on memory-integral contributions to the force on an accelerating spherical drop at low Reynolds number[END_REF]:

K µ (t, µ * → 0) = 4 3 exp 9 t t ν erfc 3 t t ν . (2.3)
Considering the transient evolution of a droplet suddenly submitted to a uniform and constant flow under creeping flow conditions, Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] have reported that the slip length along the interface reaches a uniform constant value:

λ(µ * ) = R 3µ * . (2.4)
Based on this result, the kernel expression determined for a slip sphere Gatignol [START_REF] Gatignol | The faxen formula for a rigid particle in an unsteady non-uniform stokes flow[END_REF], Michaelides and Feng [START_REF] Michaelides | The equation of motion of a small viscous sphere in an unsteady flow with interface slip[END_REF], Premlata and Wei [START_REF] Premlata | The Basset problem with dynamic slip: Slip-induced memory effect and slip-stick  transition[END_REF]

K µ (t, λ) = (1 + 2λ/R) 2 λ/R(1 + 3λ/R) exp (1 + 3λ/R) 2 λ 2 /R 2 t t ν erfc (1 + 3λ/R) λ/R t t ν . (2.5)
has been extended to the case of droplet with a finite viscosity ratio µ * as

K µ (t, µ * ) = (2 + 3µ * ) 2 3(1 + µ * ) exp 9(1 + µ * ) 2 t t ν erfc 3(1 + µ * ) t t ν . (2.6)
This expression has been found to properly reproduce direct numerical simulations of a fluid sphere in an impulsive uniform flow. The accuracy of this analytic kernel for an oscillating creeping flow is investigated here by considering the unsteady periodic flow

W (t) = W 0 [1 + sin(2πf t)] , (2.7)
where f is the oscillation frequency.

For such an oscillating flow, an exact solution of the total force 2.1 can be expressed in the frequency domain considering the solution derived by Gorodtsov [START_REF] Gorodtsov | Slow motions of a liquid drop in a viscous liquid[END_REF] and Yang and Leal [START_REF] Man | A note on memory-integral contributions to the force on an accelerating spherical drop at low Reynolds number[END_REF] investigated theoretically the motion of viscous spheres. A small error was present in those papers, see the discussion by Pozrikidis [44]. This has been corrected by Galindo and Gerbeth [START_REF] Galindo | A note on the force on an accelerating spherical drop at low-Reynolds number[END_REF]. Following the notation of Abbad and Souhar [START_REF] Abbad | Experimental investigation on the history force acting on oscillating fluid spheres at low Reynolds number[END_REF], the corresponding force expression in the frequency domain is

F (f * ) = 6πµ e R W (f * ) 2 + 3µ * 3(1 + µ * ) + K 2 o 3 + µ * 1 + µ * K o + L(K i , K o , µ * ) , ( 2.8)
where W (f * ) is the Fourier transform of the flow oscillation, and

L(K i , K o , µ * ), Q, K o and K i are given by L(K i , K o , µ * ) = 1 + 3K o 3(1 + µ * ) - (1 + K o ) 2 3 + K o + µ * Q(k i ) , ( 2.9
)

Q = K i (6 + K 2 i ) -3(2 + K 2 i )tanh(K i ) (3 + K 2 i )tanh(K i ) -3K i , ( 2.10
)

K o = -2πif * , K i = K o ρ * µ * . (2.11)
where i is the complex number. In this work, the oscillation frequency f is normalized by the diffusive time t ν = R 2 /ν e using the kinematic viscosity of the external fluid:

f * = f t ν = f R 2 ν e .
(2.12)

Similarly, the flow frequency can be normalized with the internal diffusion time:

f * i = f R 2 /ν i = f * ρ * /µ * .
In expression 2.8, the four contributions from left to right are the steady drag force, the added mass force and the inertial (pressure gradient) force combined in a single term, and the two last terms give the history force for viscous spheres as

FH (f * , µ * ) = 6πµ e R W (f * ) µ * 1 + µ * K o + L(K i , K o , µ * ) .
(2.13)

K o µ * /(1 + µ * ) is the extension of the Basset-Boussinesq term K o to fluid spheres while L(K i , K o , µ *
) is an additional term specific to fluid spheres. The limit µ * → ∞ gives the solid sphere behavior (Basset-Boussinesq history force)

FH (f * , µ * → ∞) = 6πµ e RK o W (f * ) (2.14)
while the limit µ * → 0 gives the spherical bubble behavior (Yang and Leal history force) 2.15) will be used for a direct comparison with our numerical simulations and to discuss the relevance of the kernel proposed in this study for fluid spheres.

FH (f * , µ * = 0) = 8πµ e R K o 3 + K o W (f * ) (2.
Under the flow conditions considered in this paper, the flow field around a fixed fluid sphere is a priori controlled by four dimensionless numbers: the droplet Reynolds number Re based on the physical properties of the external fluid, the viscosity ratio µ * , the density ratio ρ * respectively defined as

Re = 2Rρ e W 0 µ e , µ * = µ i µ e , ρ * = ρ i ρ e , ( 2.16)
and the normalized oscillation frequency f * . For low values of both f * and f * i , momentum diffusion dominates both inside and outside the droplet. The local velocity field is hence expected to quickly adjust to the flow oscillations with a quasi-static evolution corresponding to a weak history force contribution. For large values of f * or f * i important delays in diffusion are expected, yielding a significant contribution of the history force on the momentum balance.

In the following, we report on an extensive analysis of the quantitative contribution of the history force both in the frequency and time domains by varying the viscosity ratio µ * and the normalized frequency f * considering Reynolds number much smaller than 1. To reduce the range of parameters explored, the density ratio ρ * is set to 1, as it has been shown to have little influence on the history force Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], Feng and Michaelides [START_REF] Zhi | Heat and mass transfer coefficients of viscous spheres[END_REF]. The parametric study was carried out within the range of parameters reported in Table 2.1.

Numerical method

The simulations were carried out using the JADIM code developed at IMFT. This code, initially dedicated to the direct simulation of unsteady 3D flows around solid spheres and spherical bubbles Magnaudet et al. [START_REF] Magnaudet | Accelerated Flows Past a Rigid Sphere, Part 1[END_REF] Legendre and Magnaudet [START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF] has recently been extended to study the flow dynamics related to fluid spheres Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF] Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. We refer the reader to these references for a detailed presentation of the code, as well as its validation for the simulation of such 2.3. NUMERICAL METHOD problems. In this study, an axisymmetric polar mesh is used to solve the flow both inside and outside the droplet, with typically hundreds of cells regularly distributed along the interface. The droplet is set fixed at the origin of a domain extending from x = -50R to x = 50R, to reduce confinement effect for the considered range of Reynolds number. The Navier-Stokes equations are integrated in space using a finite volume method using a second order centered scheme, while time advancement is achieved through a three steps Runge-Kutta scheme (the nonlinear terms are computed explicitly while the diffusive terms are computed using a Crank-Nicolson scheme). The time step was chosen to satisfy both the CFL stability criterion of the numerical method and to guarantee a minimum of 100 time steps per flow oscillation period (1/f ). The simulations are conducted by first simulating the steady flow induced by a constant inlet velocity W 0 . Once a steady-state is reached, the periodic flow given by Eq. 2.7 is imposed as a boundary condition at the domain inlet (at x = -50R). Then we let the simulation run over several characteristic times t ν to reach the periodic oscillating solution.

The unsteady total force F (t) acting on the droplet is calculated by the direct integration of the pressure and viscous stresses over the interface separating the two fluids. The force decomposition into the steady drag, the added-mass force and the Basset-Boussinesq history force is described below. Considering the analytic expressions of both the steady drag force:

F SD (t) = 6πµ e R 2 + 3µ * 3 + 3µ * W 0 (1 + sin(2πf t)), (2.17)
and the sum of the inertial (pressure gradient) force and the added-mass force:

F I&AM (t) = 4π 2 ρ e R 3 f W 0 cos (2πf t) , (2.18)
the history force is determined as:

F H (t) = F (t) -F SD (t) -F I&AM (t). (2.19)
F 0 = 6πµ e RW 0 2 + 3µ * 3 + 3µ * . (2.20)
Each contribution follows a periodic oscillation due to the sinusoidal velocity forcing. When the fully periodic regime is reached, each force can be decomposed into an average and a fluctuating part.

F (t) = F + F ′ (t), (2.21)
with F the mean value of F and F ′ (t) the unsteady fluctuation over an oscillation period. This enables to characterize the intensity of the unsteady contribution F ′ (t) using the 42CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW Root-Mean-Square (RMS) metrics, defined as:

RM S(F ) = 1 n n i=1 (F ′ (i)) 2 . (2.22)
where n is the number of time steps over the considered period. ). All forces are scaled by the mean steady drag force over the period F 0 given by Eq. 2.20.

Flow field description

The streamlines inside and around the droplet are reported in Fig. 2.2 at different times of the oscillation period, in the case Re = 0.1, µ * = 0.2 and f * = 1. In this configuration, the external fluid is more viscous than the fluid inside the droplet. Over the time period, the structure of the flow has been observed to evolve drastically over the time period. In particular, counter-rotating vortices appear both inside and outside the droplet.

At the beginning of the oscillation period (Fig. 2.2a), the streamlines are similar to the streamlines of a steady uniform Stokes flow around a viscous drop: a Hill's spherical vortex turns clockwise inside the droplet entrained by the viscous stress at the interface exerted by the external flow moving from left to right. Linearity of Stokes equations yields fore-aft symmetry of the creeping flow streamlines. This streamlines structure is observed until the inlet velocity reaches its maximum W (t) = 2 W 0 at t = T /4 (Fig. 2.2b), and also when it decelerates down to W (t) = W 0 at t = T /2. Just after the half period at time t = 0.554T when W (t) = 0.667 W 0 , Fig. 2.2d shows the appearance of a small recirculation zone located in the near wake of the drop.

Under oscillating flow conditions around a solid sphere at Reynolds varying from 0.1 to 40, similar phenomena have already been observed during the phase of the flow deceleration Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF], Chang and Maxey [START_REF] Chang | Unsteady flow about a sphere at low to moderate reynolds number. part 1. oscillatory motion[END_REF]. The flow separation has been observed in the decelerating phase of the oscillation: under strong decelerating conditions the flow can completely detach itself from the sphere and a reverse flow surrounds the sphere.

Then, this recirculation region grows along the upstream direction, its size increases and finally completely covers the droplet(Fig. 2.2e). Meanwhile, a second vortex develops inside the droplet to match this clockwise external recirculation with the internal flow, while the size of the external vortex increases along the radial direction (Fig. 2.2f).

The second internal recirculation vortex increases from t = 0.582 T to t = 0.8 T (Fig. 2.2h), while the first one vanishes, resulting in a Hill spherical vortex now rotating counter-clockwise. At t = 0.883 T , a new recirculation appears close to the front stagnation point (Fig. 2.2i), and develops between the droplet and the clockwise oriented external recirculation, which eventually vanishes at the end of the period. The evolution of the tangential velocity U θ along the droplet interface, normalized by W 0 is reported in Fig. 2.3 at the same times as in Fig. 2.2. We can observe that the sign of U θ is changing during the period, between t = 0.582 T and t = 0.8 T , in agreement with the development of the external clockwise and internal counter-clockwise vortices. The evolution of the longitudinal velocity U x (x = 0, r) is reported as a function of the radial distance r to the symmetry axis. Due to the development of the two vortices, the direction of U x changes considerably over the period. The affected region extends for a distance of about 3R around the droplet interface. Beyond this region, the flow direction is controlled by the upstream velocity W (t) prevailing far from the drop. The influence of the viscosity ratio on the flow structure is illustrated in Fig. 2.4 for Re = 0.1, µ * = 5 and f * = 1. The viscosity ratio is the inverse of that considered in Fig. 2.2 which means that the external fluid is now less viscous than the fluid inside the droplet. A flow structure similar to Fig. 2.2 is observed over the period, with the development of external clockwise and internal counter-clockwise vortices. However, recirculations appear earlier and grow faster for µ * = 5. One explanation is that the presence of a more viscous liquid inside the droplet causes the continuous phase motion to remain controlled by the inner fluid motion as a consequence of the continuity of the tangential velocity at the droplet interface. To consider now the effect of the oscillation frequency, Fig. 2.5 and Fig. 2.6 compare the streamlines obtained inside and around the droplet at f * = 10 -2 and f * = 10 3 , respectively. The Reynolds number and viscosity ratio are the same as in Fig. 2.2 (Re = 0.1 and µ * = 0.2). A significant effect of f * can be observed on the development of the flow. For f * = 10 -2 , a quasi-static evolution of the flow was expected, with a flow field corresponding at each time to that of a droplet under steady condition. However, the external vortices are still observed, revealing a strong coupling between the flow oscillation and the flow field response. An evolution similar to the one reported for f * = 1 is observed, with the development of two vortices inside and outside the drop, but over a shorter part of the period (Fig. The fraction of the period during which the external recirculation is observed is a relevant indicator to quantify the impact of the oscillation on the flow field change. When the imposed velocity oscillates at moderate frequency, f * = 1, the external recirculation is observed during 43% of the period duration. When reducing the frequency, the flow seems to evolve toward a mostly quasi-static pattern, a recirculation being observed during 13% of the period duration for f * = 0.01. At last, in the case of high-frequency fluctuations, f * = 10 3 , the recirculation is observed up to 50% of the period.

a) t = 0T -W = W 0 (b) t = 0.25T -W = 2W 0 (c) t = 0.5T W = W 0 (d) t = 0.554T -W = 0.667W 0 (e) t = 0.582T -W = 0.507W 0 (f) t = 0.670T -W = 0.124W 0 (g) t = 0.75T -W = 0W 0 (h) t = 0.8T -W = 0.0489W 0 (i) t = 0.883T -W = 0.
(a) t = 0T -W = W 0 (b) t = 0.663T -W = 0.146W 0 (c) t = 0.665T -W = 0.139W 0 (d) t = 0.666T -W = 0.134W 0 (e) t = 0.75T -W = 0W 0 (f) t = 0.803T -W = 0.0549W 0

Discussion on force contributions

We now consider the force decomposition as introduced in the section Numerical method. Due to the flow reversibility and to the linearity of Stokes equations, the mean history 2.5. DISCUSSION ON FORCE CONTRIBUTIONS force (2.17) is zero. At low Re, there is no coupling between the effects of the steady component of the flow and its unsteady periodic oscillation. Therefore, in the fully periodic regime, the average total force is equal to the steady drag force F 0 corresponding to the mean flow W 0 . From the analytical expression of the pressure gradient and addedmass force, it is obvious that the time average of this force is zero as well, such that the mean contributions of the force are:

F tot = F SD = 6πµ e RW 0 2 + 3µ * 3 + 3µ * , F I&AM = 0, F H = 0. (2.23)
According to our simulations, these conditions are satisfied, the measured values of F I&AM and F H being both lower than 3 × 10 -3 F 0 .

To compare the relative contributions of each force, we consider their root-mean-square (RMS) magnitude over one oscillation period. For given Reynolds number and viscosity ratio, the RMS values of the history force, inertial and added-mass force and therefore the total force amplitude depend on f * while the RMS of the drag force has a constant value because it does not depend on the oscillation frequency f . The results are shown in figure 2.8, where we use the mean value of the steady drag force F 0 to scale all the RMS contributions. The history force RMS is observed to increase with f * and to become larger than the steady drag force for f * ≥ 1. The inertial (pressure gradient) and added-mass forces RMS are also increasing with f * but with a faster rate so that it controls the total force magnitude for f * > 30. As a consequence, the fluctuating contributions of both the steady drag and history forces can be neglected at large frequencies. Indeed, it is possible to make the forces dimensionless using 6πµ e RW 0 and time using f -1 in equation (2.1). Therefore, the drag force becomes of order unity and the added mass force and pressure gradient term scale with f * dominating the total force for rapid 50CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW oscillating flows. Depending on the asymptotic behavior of the kernel at high frequency the history force scales differently. For solid particles, the Basset history force scales with (f * ) 1/2 whereas it becomes constant for bubbles. For droplets with finite viscosity ratio, the high frequency regime evolves between these two limits.

To gain a quantitative analysis of the history force contribution to the total force balance, we seek the conditions yielding the history force RMS to be more than 10% of the total force RMS. To do so, we plot the relative contribution (to the total force RMS) of each force RMS in Fig. 2.9 as a function of f * . As observed in the figure, the oscillating contribution of the history force appears to be significant when 0.01 ≤ f * ≤ 10.

History force in the frequency domain

In this section, we discuss the evolution of the history force FH in the frequency domain. For that purpose we introduce the normalized transfer function F defined as FH = 6πµ e R W F. The magnitude of the transfer function F obtained from our numerical simulations is compared in Fig. 2.10 with the expressions (2.14) (2.15) and (2.13) for a solid sphere, a spherical bubble and a spherical fluid particle, respectively. The viscosity ratios µ * = 0.05, µ * = 0.2 and µ * = 20 are shown. As observed, our numerical simulations are in very good agreement with the exact solution (2.13) for a fluid sphere over a wide range of viscosity ratio and oscillation frequency. At low frequency, F evolves as √ f * in all cases. However, at high frequency, the transfer function for a spherical bubble tends to the value 4/3 while both our simulations and the exact relation (2.13) still evolve as √ f * . In fact, much lower values for µ * are required (typically lower than 10 -3 ) to be able to use the Yang & Leal history force relation at high frequency.

The transfer function of the history force expression (2.6) proposed by Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] is reported in Fig. 2.10. The agreement is found to be very good at low frequency but a clear difference is observed at large frequency where expression (2.6) follows the behavior observed for a bubble. This result is consistent with the numerical simulations of Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] that show that the the early-time transient evolution of the slip length (corresponding to large values of f * ) can not be described by λ = R/3µ * (relation 2.4). Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] have proposed a correlation for the unsteady slip length which matches the numerical results

λ(µ * , t) = R 3µ * 1 -exp - 60 t t νi 0.55 cos 20 t t νi . (2.24)
where t νi is based on the droplet kinematic viscosity ν i . As shown in Fig. 2.10, the evolution for large f * is drastically improved when considering Eq. 2.24 instead of λ = R/3µ * in the Kernel expression 2.6 and agrees very well with both the theoretical prediction and direct simulation results.

Kernel for history force in the time domain

In the context of particulate turbulent flows, the dynamics of a dispersed phase made of droplets is usually characterized in terms of the Stokes number St = τ p /τ f that compares the droplet relaxation time τ p to a relevant fluid flow timescale τ f . If St >> 1, the trajectory of the droplet is only weakly affected by the flow structures of characteristic time τ f . Otherwise, if St << 1, the droplet evolves as a fluid tracer closely following the flow streamlines and can easily be captured by any vortex. In the intermediate regime of finite Stokes number, preferential accumulation is observed either within vortices for bubbles or in high strain regions for inertial solid particles. At small particle Reynolds number Re = ρ e dW 0 /µ e ≪ 1, τ p is given by:

τ p = 2 9 (ρ i + C M ρ e ) µ e 2 + 3µ * 3 + 3µ * R 2 , (2.25)
with C M the added-mass coefficient equal to 1/2 for a sphere. In the flow configuration considered here (see Eq. 2.7), the flow timescale is the oscillation period τ f = T and the Stokes number reads: The history force RMS contribution is found to be more than 10 %, for Stokes number in the range 0.01 < St < 5. Similar results have been observed for radial migration of particles in a vortex for St < 0.01 with a significant history force contribution Candelier et al. [START_REF] Candelier | On the effect of the boussinesq-basset force on the radial migration of a stokes particle in a vortex[END_REF]. At lower Stokes numbers, the history force RMS is almost equal to the steady drag RMS contribution, while at high frequency the pressure gradient and added-mass force are the dominant forces. However, the range of parameters (Stokes number and viscosity ratio) over which the history force must be considered is wide and corresponds to the regime where preferential accumulation is expected (0.05 < St < 5) for droplets (0.5 < µ * < 10). Therefore, an accurate expression for the history kernel able to calculate the unsteady force for Lagrangian tracking is highly desirable. For that purpose, the expression of the history force

St = τ p T = 2 9 (1 + ρ * /2) 2 + 3µ * 3 + 3µ * f * . ( 2
F H = 6πµ e R t 0 dW dt ′ K λ (t -t ′ , λ) dt ′ , ( 2.27)
is computed using Simpson's method to guarantee second order accuracy. In the kernel expression given by Eq. 2.5, multiplications between exponential and complementary error functions can cause computational divergence because each term can easily overpass the maximum number of digits a computer can store, even if they compensate each over (the first one tends to infinity while the second one tends to zero). A solution to overcome this problem proposed by Daitche [START_REF] Daitche | Advection of inertial particles in the presence of the history force: Higher order numerical schemes[END_REF] is described in the appendix, and used here for the computation of the history force.

The interfacial slip length is calculated by

λ = U θ,e r ∂ ∂r U θ,e r (2.28) 54CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW
with U θ,e the tangential velocity of the external fluid at the interface and r the radial position from the droplet center.

The interfacial slip length measured in the oscillatory flow considered here is reported as a function of the polar angle θ in Fig. 2.12a. We can observe that the interfacial slip length 2.28 is uniform along the interface all over the oscillation period for the different cases considered, i.e. when varying both the viscosity ratio and the normalized frequency. We find a maximum deviation of λ along the interface of 1.5%, except when the flow changes direction because division by a zero shear stress occurs, see Fig. As Eq. 2.4 is relevant in most considered cases, it is introduced in Eq. (2.27) to calculate expression 2.5. The resulting history force is compared in Fig. 2.13 with direct numerical results for some reference cases. We observe a good agreement as long as f * ≤ 100. Otherwise, for f * > 100 it appears that the model prediction is no longer in phase, and underestimates the amplitude of oscillation. The corresponding difference is quantitatively measured for all the cases with the Mean Absolute Error (MAE):

MAE * (F H model , F H ) = 1 A(F h )n n i=1 |F H model i -F H i | . (2.29)
where n is the number of time steps.

The resulting difference map is shown in Fig. 2.14 (left). As expected, the difference between the model and the numerical simulations increases when increasing the Stokes number (i.e. the frequency f * ) for model with slip length expression 2.4, specially at low viscosity ratio. Using expression (2.24) for the slip length gives much better agreement with the simulation results at high frequency, see Fig. 2.10 and Fig. 2.15, over all the range of Stokes number and viscosity ratio. The discrepancy is only 4% at low viscosity ratio when the flow recirculation appears and develops more in the continuous phase. When comparing the plot 2.14, to the plot reporting the relative contribution of the history force to the total force (see Fig. 2.11) we observe that the error in calculating the history force collapses to the region where its contribution to the total force is not relevant. Thus, the history force expression given by Eq. (2.27-2.5) with slip length given by ( 2.4) is accurate when its contribution is of importance for the total force, as indicated by Fig. 2.14 (right). Interestingly, even if we cannot correctly predict the history force for some range of parameters, in those cases the history force contribution to the total force balance is negligible. As a consequence, the proposed history force expression given by Eq. (2.27-2.5) with slip length (2.4) can be used for any set of parameters and will provide a significant improvement of droplet trajectory prediction when the history force contribution cannot be neglected. An even better accuracy can be achieved over the full range of physical parameters we investigated by using the history force kernel given by Eq. (2.27-2.5) with an unsteady slip length ( 2.24). ) and equation (2.5) ( ) with the kernel given by Eq. (2.5), the slip length by relation( 2.4) and ( ) with the kernel given by Eq. (2.5) and the slip length by relation ( 2.24) . (a) Re = 0.01, Contour is created using (Delaunay et al. [START_REF] Delaunay | Sur la sphere vide[END_REF]) triangulation between data points.

µ * = 5, f * = 0.1 (f * i = 0.02), r 2 = 0.9619, (b) Re = 0.1, µ * = 0.2, f * = 1 (f * i = 5), r 2 = 0.9452, (c) Re = 0.01, µ * = 0.05, f * = 100 (f * i = 2000), r 2 = 0.1970.

Conclusion

Using direct numerical simulations, we have carried out a parametric study to characterize the history force acting on a droplet held fixed in an oscillatory uniform flow at low Reynolds number. To place such results in the context of turbulent dispersed two-phase flows, we analyzed our results in terms of the Stokes number and the viscosity ratio. Based on simulation results, the history force contribution found to be significant (typically larger than 10% of the RMS total contribution) for a wide range of Stokes number (typically between 0.5 and 5) and viscosity ratio (between 0.5 and 10) where all important dynamics of turbulence induced clustering is expected. Neglecting history force could lead to significant errors in the force balance, which makes it highly desirable to get an analytic history kernel for droplet trajectory equation.

58CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW

Although an analytical expression of the history force kernel exists in the frequency domain, it cannot be transformed in the time domain for the general case of droplets (finite viscosity ratio). A model of the history force for droplets in the time domain has already been presented by Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF], and validated under creeping impulsive flow conditions. This model takes advantage of the analogy of the boundary condition existing at the droplet interface and the slip length boundary condition for superhydrophobic surface. We have tested this model, assuming a uniform and constant slip length at the interface in an oscillatory flow. We have obtained good agreement between the model and the DNS results for Stokes number below 10. It shows good agreement with the analytical expression in the frequency domain for low to moderate frequency oscillations. Over a wide range of parameters relevant to most applications, we observe that the discrepancy between the history force model with constant slip and uniform slip length and the DNS results is small when the history force contribution is significant. Using a correlation fitted on simulation data to model the early time evolution of the slip length provides an even better accuracy over the whole frequency domain. Therefore, we propose a Basset-Boussinesq history force expression for droplets that can be used for any set of parameters, making possible an accurate calculation of droplet trajectory when the history force contribution cannot be neglected.

In addition, we report an appropriate method for the Kernel calculation, which allows easier computations and avoids overflow error from exponential and complementary error function terms.

Outlook, accuracy of history force kernel in a turbulent flow

The flow in a pulsed liquid-liquid extraction column is indeed oscillatory, but based on the characteristic length and velocity scales, turbulence is most likely prevailing. However, the flow conditions surrounding the droplet surface are related to low to moderate Reynolds number based on the slip velocity. At the local scale, it is reasonable to conclude that the droplet Reynolds number is relatively low, and the external flow can be assumed as a uniform flow that fluctuates only in time. In this section, we aim at testing the accuracy of the history force kernel for droplets in a turbulent flow with a fluctuating velocity time signal containing a broad range of frequencies.

The time signal is obtained from the results of a Lagrangian simulation of particles seeded in a homogeneous isotropic turbulent flow simulated through DNS from Zhang et al. [START_REF] Zhang | Model for the dynamics of micro-bubbles in high reynolds number flows[END_REF], Zhang et al. [START_REF] Zhang | Fluid inertia effects on the motion of small spherical bubbles or solid spheres in turbulent flows[END_REF] and Gorokhovski and Zamansky [START_REF] Gorokhovski | Modeling the effects of small turbulent scales on the drag force for particles below and above the kolmogorov scale[END_REF]. The slip velocity time signal is used as the velocity inlet condition for our simulations, instead of the periodic oscillatory signal previously used. The time evolution of the velocity and its Fourier energy spectrum (to compute the Stokes number) are shown in Fig. 2.16.

OUTLOOK, ACCURACY OF HISTORY FORCE KERNEL IN A TURBULENT FLOW59

(a) Time evolution of the slip velocity.

(b) Velocity signal in the Fourier frequency domain.

Figure 2.16

Transforming the signal in the Fourier frequency domain, see Fig. 2.16b, yields the range of variation for the Stokes number of droplets. We consider the relevant range of Stokes number for all frequencies corresponding to at least one tenth of the maximum energy content. It corresponds to a frequency range from one to twenty-five in arbitrary units.

From the Stokes number definition in Eq. 2.26, we obtain the following ranges for each viscosity ratio we tested:

• For µ * = 0.2, St ∈ [0.120, 3.01] • For µ * = 0.5, St ∈ [0.130, 3.24] • For µ * = 2, St ∈ [0.148, 3.70] • For µ * = 5, St ∈ [0.157, 3.94]
In Fig. 2.17, we can observe the comparison between the history forces obtained from the numerical simulations and predicted by the model with equations (2.27-2.5) using the varying slip length given by equation (2.24). The results are extremely accurate in accordance with Fig. 2.18 considering the range of variation of the Stokes number for each viscosity ratio. We obtain an excellent agreement between the history force from simulations and the model (see zoom sections of the time evolution). Similarly to Fig. 2.17, we compare in Fig. 2.19 the total force obtained from simulations with the model (Eq. 2.1) using the history force kernel. The model results are again in excellent agreement with simulations.

INTRODUCTION

65

In this chapter, we study the 3D bifurcations of the flow occurring inside and outside the interface of a spherical droplet in a uniform flow. The investigation uses both experimental study and numerical simulations to find out how the hydrodynamic configuration affects the 3D bifurcation and therefore the drag coefficient.

We were able to distinguish two distinct types of bifurcation, depending on whether they originate from internal or external flow instabilities. By analyzing the amplification rate of a flow disturbance, we have interpolated critical bifurcation thresholds, leading to a precise bifurcation criterion that depends on the Reynolds number and on the viscosity and density ratios.

These bifurcations can be characterized as steady or unsteady, each having a different influence on the drag force. Further analysis of the impact of the bifurcation on drag coefficient leads to a comprehensive drag prediction method for a 3D regime, based on the bifurcation criterion and a drag coefficient amplification factor.

Experimental results, obtained on droplets meeting the sphericity criterion, confirm the simulation results, and demonstrate that a 3D bifurcation can be observed under practical operating conditions. This chapter will be submitted soon as a research article.

Introduction

The flow patterns and hydrodynamic instabilities related to moving solid spheres, gas bubbles, or liquid droplets in a continuous phase have been extensively studied due to their relevance in many engineering systems and natural phenomena. These studies have focused on fixed objects or freely rising or falling particles, which exhibit different types of flow instabilities depending on the geometrical configuration and the number of degrees of freedom (translation, rotation).

The configuration of solid spheres held fixed in a uniform stream has been studied precisely from the early nineties by Johnson and Patel [START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF] who conducted both experiments and numerical simulations. A series of simulations by Tomboulides and Orszag [START_REF] Tomboulides | Numerical investigation of transitional and weak turbulent flow past a sphere[END_REF], Natarajan and Acrivos [START_REF] Natarajan | The instability of the steady flow past spheres and disks[END_REF] and Fabre et al. [START_REF] Fabre | Bifurcations and symmetry breaking in the wake of axisymmetric bodies[END_REF] complement the scenario of the successive flow instabilities. The wake flow remains axisymmetric up to a Reynolds number of approximately Re = 212. Beyond this critical threshold, the flow loses its axial symmetry and bifurcates to a steady flow with a single plane symmetry. As Re is further increased, a periodic flow is observed in the wake due to a Hopf bifurcation, which finally becomes unstable for Re greater than 300.

This transition scenario is similar for fixed or freely moving solid spheres. However, when particles are freely moving, the loss of axisymmetry in the wake flow at the first bifurcation generates a lift force which yields oblique trajectories. Jenny et al. [START_REF] Jenny | Instabilities and transition of a sphere falling or ascending freely in a Newtonian fluid[END_REF] found that this first bifurcation is weakly influenced by the solid to fluid density ratio (ρ * ). The 66CHAPTER 3. 3D BIFURCATIONS STUDY second bifurcation resembles that of a fixed sphere only if the density ratio is sufficiently high. Trajectories of light particles can alternate between in-plane zigzag patterns and fully three-dimensional trajectories, while the trajectories of denser particles directly exhibit three-dimensional complex trajectories right after the oscillatory oblique regime.

The specific case of fluid particles, such as bubbles and drops, is more complex due to their deformable shape and the need to consider the velocity fields on both sides of the interface. Experimental studies of Magnaudet and Eames [START_REF] Magnaudet | The motion of high-reynolds-number bubbles in inhomogeneous flows[END_REF], Prosperetti [61] have shown that small gas bubbles rise on rectilinear paths aligned with gravity, while larger ones can exhibit spiral or zigzag trajectories. In clean water, the wake instabilities for solid spheres are remarkably similar to those observed for bubbles. However, the stability behavior of bubbles depends on more than just the Reynolds number when considering free-slip boundaries. Deformations of the bubble shape is also a possible cause for the regular mode to become unstable. In Bonnefis et al. [START_REF] Bonnefis | When, how, and why the path of an air bubble rising in pure water becomes unstable[END_REF] and Herrada and Eggers [START_REF] Herrada | Path instability of an air bubble rising in water[END_REF], the path of a millimeter-sized air bubble rising in water becomes unstable when it reaches a critical size of approximately 0.93 mm (same value for both articles, at this point the bubbles are no more spherical). The authors state that the instability is caused by the coupling of the body and fluid motions. These instabilities create a time-varying force on the bubble that can cause it to deviate from its straight path. Whereas for a perfectly spherical bubble, simulations revealed that the flow remains axisymmetric at all Reynolds numbers Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF].

The case of viscous droplets has been much less studied by rigorous linear stability analysis and identification of successive bifurcations. At low Reynolds numbers (small spherical droplet) the flow solution is given by Stokes equations and allowed Schiller and Naumann [START_REF] Schiller | Über die grundlegenden berechnungen bei der schwerkraftaufbereitung[END_REF], Rivkind and Ryskin [START_REF] Ya Rivkind | Flow structure in motion of a spherical drop in a fluid medium at intermediate Reynolds numbers[END_REF] and Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF], to derive theoretical prediction for the drag coefficient. Some studies dedicated to measuring the terminal drop velocity (to assess correlations for the drag coefficient over a wider range of Reynolds number and viscosity ratio), have observed oscillatory behavior for large spherical or deformed droplets, as in Thorsen et al. [START_REF] Thorsen | On the terminal velocity of circulating and oscillating liquid drops[END_REF]. There are often questions about the effect of surface contamination that may significantly alter droplet dynamics. At low Reynolds number, the fluid flow inside a droplet is characterized by the presence of a single toroidal vortex, called Hill's vortex (1894) Hill [START_REF] Hill | On a spherical vortex[END_REF]. This configuration, while being an analytical solution under Stokes flow condition, has been observed over a large range of Reynolds number. Hill's vortex stability has already been studied in an unbounded fluid by Moffatt and Moore [START_REF] Moffatt | The response of Hill's spherical vortex to a small axisymmetric disturbance[END_REF], Amick and Fraenkel [START_REF] Amick | The uniqueness of hill's spherical vortex[END_REF] and Pozrikidis [START_REF] Pozrikidis | The nonlinear instability of Hill's vortex[END_REF]. The initial spherical shape of the vortex is therefore free to deform which is prone to hydrodynamic instabilities. Those findings do not directly apply to our investigation since Hill's vortex is confined within the droplet surface.

The interest for a better understanding of the dynamics of liquid-liquid systems is growing related to liquid-liquid extraction processes (CO2 sequestration and storage, new solvents for the nuclear power industry), oil spills, agro-food industry... Recent studies addressed new issues concerned with droplet hydrodynamics at moderate to high Reynolds number. For high viscosity ratios, Albert et al. [START_REF] Albert | Dynamic behaviour of buoyant high viscosity droplets rising in a quiescent liquid[END_REF] observed by means of VOF numerical simulations that spherical droplets rising in oblique paths have a bifid wake similar to those observed for solid sphere first flow bifurcation. For a wider range of viscosity ratio, we can note in particular the work of Edelmann et al. [START_REF] Edelmann | Numerical investigation of different modes of internal circulation in spherical drops: Fluid dynamics and mass/heat transfer[END_REF] and Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF]. In both those studies, the structure and instability of the axisymmetric Hill's vortex were discussed. For some specific Reynolds number and viscosity ratio, a very different flow dynamics is observed inside the drop: two counter-rotating vortical structures are observed in the xy plane (y and z defining any plane perpendicular to the main flow x-direction). In the yz plane, four counter-rotating vortical structures are also observed.

Those flow structures can be stable, and therefore reach a steady state or unstable, and then flow oscillations are observed reminiscent of a Hopf flow bifurcation. The external wake flow may have its own hydrodynamic instabilities, leading to a complex interplay between inner and outer flow dynamics.

The objective of this work is to investigate the stability of the flow inside and outside the interface of a steadily moving spherical droplet. We will be considering both experiments and numerical simulations to disentangle the origin of the hydrodynamic instabilities. The chapter is organized as follows: we present first our experimental set-up and the numerical methods for the simulation. Then, experiments will be reported and discussed compared to numerical simulations to support and complement the experimental results.

Problem statement

For the experiment, we consider a fluid sphere of radius R, density ρ i and viscosity µ i moving at a constant velocity -→ U in a quiescent infinite fluid of density ρ e and viscosity µ e . Note that the simulations were conducted in the reference frame moving with the droplet, and a constant velocity was imposed at the inlet of the numerical domain. This problem, is characterized by several dimensionless numbers. First, the instantaneous Reynolds number Re, by convention the external Reynolds number, based on the magnitude U of the droplet velocity.

Re = 2ρ e U R µ e (3.1)
In some discussions, the internal Reynolds number Re i can be used:

Re i = 2ρ i U R µ i (3.2)
In the case of unsteady motion, an average Reynolds number Re based on U the mean value of U is introduced:

Re = 2ρ e U R µ e . ( 3.3)
The viscosity and density ratios, µ * and ρ * , compare the viscosity and density inside 68CHAPTER 3. 3D BIFURCATIONS STUDY and outside the droplet, respectively.

µ * = µ i µ e (3.4) ρ * = ρ i ρ e (3.5)
This study focuses on the characterization of the droplet vertical motion and in particular the drag force. The latter is directly obtained from the DNS results, by integration of the total stress, composed by the pressure p e and the viscous stress Σ e of the steady flow solution for the external fluid, on the droplet interface:

F d = S (-p e - → n + Σ e . - → n ) dS (3.6)
where n is the local unit vector normal to the droplet's interface.

Regarding experimental results, the drag force is determined from post-processing the droplet 2D trajectory and in particular its rising velocity. The illustration of a nonrectilinear path is given in Fig. 3.1, highlighting the oscillations of the droplet velocity in the vertical direction. The fluid sphere motion is described by the following force balance equation:

ρ i V dU dt = ρ i V g -ρ e V g -C m ρ e V dU dt + F d (3.7)
where dU dt is the droplet acceleration, g the gravity acceleration, C m the added mass coefficient (C m = 1/2 for a sphere), and V = 4 3 πR 3 the volume of the droplet. The forces on the right-hand side of Eq. (3.7) are respectively the weight and Archimedes forces, the added mass effect, and F d (t) that contains both the classical steady drag force and some history effects resulting from unsteady motion. In the following, for simplification, this force will be called the drag force and normalized to introduce the drag coefficient as follows:

F d = -C d πR 2 2 ρ e |U |U, (3.8)
Eq. 3.7 then yields:

(ρ i + C m ρ e )a = (ρ i -ρ e )g -C d 3π 8R ρ e |U |U (3.9)
whose various contributions become, after projection along the vertical direction e y :

(ρ i + C m ρ e )    a x a y a z    = -(ρ i -ρ e )    0 g 0    -C d 3π 8R ρ e |U |    U x U y U z    (3.10)
with a = dU dt . For rectilinear motion (or a fixed sphere in a uniform flow as considered in the simulations), all terms are aligned along the vertical direction of motion e y , whereas for an oscillating trajectory, the drag force may contain contributions along e x and e z . By projecting Eq. (3.10) along the vertical y-direction, we can express the drag coefficient C d as a function of the velocity, acceleration and physical properties:

C d = 8R [(ρ e -ρ i )g + (ρ i + C m ρ e )a y ] 3ρ e |U |U y (3.11)
According to our experimental data, the (ρ e -ρ i )g term in Eq. (3.11), is on average 65 times greater than the (ρ i + C m ρ e )a y term. It is therefore possible to neglect the acceleration term and consider the classical expression, Eq. (3.12), with an average median error of less than 1.2%.

C d = 8R(ρ e -ρ i )g 3ρ e |U |U y (3.12)
Droplet deformation is quantified using the aspect ratio χ, that compares the major axis α to the minor axis β of the droplet:

χ = α β (3.13)
To assess droplet sphericity, we also consider the degree of true sphericity Φ defined by Wadell [START_REF] Wadell | Shape, and Roundness of Quartz Particles[END_REF] as the ratio between the droplet projected surface area measured in the images and the area of a sphere of the same volume: The droplet is considered spherical if Φ > 0.99 and χ < 1.1.

We also introduce two other dimensionless numbers, the Morton number, M o, and the Eötvös number, Eo, to characterize droplet motion in experiments.

M o = gµ 4 e ∆ρ ρ 2 e σ 3 , (3.15)
where σ is the surface tension in (N/m) and ∆ρ = ρ e -ρ i .

Eo = 4g∆ρ R 2 σ , ( 3.16)
The frequency f of the 3D oscillations is made dimensionless by the Strouhal number St:

St = 2Rf U , ( 3
.17)

Experimental methods

Experimental setup

Droplet trajectories were observed in a rectangular tank with external dimensions 20 × 10× 40 cm 3 and a capacity of 8 liters, fitted with two optical quality windows (see Fig. 3.2). Two-dimensional projections of the droplet were acquired using a backlight setup to avoid perspective effects. A Photron FASTCAM Mini UX100 CMOS camera was attached to an Opto-engineering TC16M096 bi-telecentric lens placed at around 25 cm from the center of the tank. The tank and sensor were uniformly illuminated with an Opto-engineering collimated green light source (peak wavelength 525 nm; beam diameter 120 mm). The images are composed of 1280 × 1024 squared pixels, corresponding to a field of view of 33.7 × 29.9 mm 2 (scale factor 2.6 10 -2 mm/px), acquired with a frame rate of 2,000 frames per second and a shutter speed of 1/51,200 s. The droplets were injected at the bottom of the tank using a needle connected by a T-junction to two precision syringe-pumps, one for the continuous phase and one for the dispersed phase.

The precision syringe-pumps allow for accurate adjustment of the droplet volume. The droplets are observed at 30 cm above the injection point, once they have reached their terminal rising velocity.

Fluids selection

The fluids were selected so that high Re could be achieved while maintaining spherical droplets, i.e. in the range of M o, Eo and Re delimiting spherical droplets in the Clift diagram of Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF]. Indeed, as discussed later in Sec. 3.6.3, high Reynolds, typically in the order Re = O(100), are required to observe 3D effects.

For simplicity, filtered deionized water was chosen for the continuous phase, while for the droplets candidate fluids were selected from the Pubchem database Kim S [START_REF] Cheng T Kim | Pubchem in 2021: new data content and improved web interfaces[END_REF]. All the fluids in the database, with a solubility lower than 3% in water, have been tested in terms of maximum achievable Reynolds number for a spherical droplets (Clift diagram) and appearance of 3D flow oscillations in the simulations (see Sec. 3.6.3).

Three organic phases were selected: n-Pentane, Heptane, and Cyclohexane. Their properties were measured using the equipment presented previously, and are summarized in Tab. 3.1. For each fluid, the maximum droplet radius allowing a spherical shape to be maintained was deduced from the Clift diagram (see Fig. 3.3), and the Tadaki criterion, T a = ReM o 0.23 , was calculated. For rising droplets, Myint et al. [START_REF] Myint | Shapes of Single Drops Rising Through Stagnant Liquids[END_REF] stated that the droplets remain spherical if T a < 2. A total of 450 droplets were recorded, with 422

Droplet detection and post-processing

Post-processing is carried out on images in "tif" format. The droplet is detected thanks to the Python scikit-image module (see Van der Walt et al. [START_REF] Van Der Walt | scikit-image: image processing in python[END_REF]) based on a light threshold. This function returns many properties, including the coordinates x and y of the droplet's centroid, its major α and minor β axis lengths, orientation, and projected area, which are of interest in this study.

From these properties, and given the image acquisition frequency, the velocity, and acceleration components and their magnitude can be derived at each instant in the xy plane of the images. Note that we cannot access information in the z direction with the setup used.

The estimate of post-processing errors is made following Cai et al. [START_REF] Cai | Hydrodynamic behavior of a single bubble rising in viscous liquids[END_REF], Celata et al. [START_REF] Piero Celata | Measurements of rising velocity of a small bubble in a stagnant fluid in one-and two-component systems[END_REF] and Tan et al. [START_REF] Tan | Experimental study of unsteady drag coefficient of droplets in a liquid-liquid system[END_REF]. Relations to estimate uncertainties are based on the pixel size p and the time step ∆t. They are respectively given below for the radius uncertainty ∆R: With:

∆R = 2p (3.
∂U ∂x = - p ∆t , ∂U ∂y = p ∆t , ∂U ∂∆t = y -x ∆t 2 p, ∂U ∂p = y -x ∆t
The uncertainty on the drag coefficient, estimated from Eq. (3.11) is estimated by taking the worst-case error. We consider the density, viscosity, radius, and velocity errors that result in the largest numerator and the smallest denominator when computing the experimental drag coefficient. Then we take the absolute value of the difference between this value and the actual computed experimental drag coefficient. We obtain for ∆Cd exp :

∆C d,exp =|C d,exp - -8(R + ∆R)9.81(ρ i + ∆ρ i -ρ e -∆ρ e ) 3(ρ e -∆ρ e )(U y -∆U ) 2 | (3.20)
If during its ascension, it appears that a droplet does not satisfy the sphericity criteria defined by Eqs. (3.13) and (3.14), this droplet is deleted from the database of experimental results.

Numerical simulations

Numerical methods

Direct numerical simulations (DNS) were carried out using JADIM, the in-house IMFT code. In the configuration studied, the droplet is set fixed at the origin of a domain extending from x = -50R to x = 50R, that allows to reduce the confinement effect for the considered range of Reynolds number.

The computational domain is discretized using a polar mesh inside the sphere, whereas for the external fluid, a mesh based on the equipotentials and streamlines of the potential flow around a cylinder is used Legendre and Magnaudet [START_REF] Legendre | The lift force on a spherical bubble in a viscous linear shear flow[END_REF]. The mesh is refined at the interface.

A 2D mesh is first produced, then rotated around the streamwise x-axis to become 2D axisymmetric (one single cell in the azimuthal direction) and 3D (N ψ cells in the azimuthal direction, see Fig. 3.5). In total, twelve parameters are necessary to characterize the mesh (see Tab. 3.2 and Fig. 3.5c). A convergence study was carried out to determine the best set of parameters to ensure that the results are mesh-independent while minimizing CPU resources requirements, important for these 3D simulations. In The Navier-Stokes equations, Eqs. (3.21) and (3.22), are integrated in space using a finite volume method with a second order centered scheme. Time-advancement is achieved 76CHAPTER 3. 3D BIFURCATIONS STUDY through a three-step Runge-Kutta scheme. The nonlinear terms are computed explicitly, while the diffusive terms are computed using a Crank-Nicholson scheme.

∇ • U δ = 0 (3.21) ρ δ ∂U δ ∂t + U δ • ∇U δ = -∇P δ + 2µ δ ∇ • S δ (3.22)
where

S δ = ∇U δ + ∇U δ ⊤ /2
is the rate of strain tensor, and P the pressure field.

At the liquid-liquid interface, the Navier-Stokes equations relating to the inner and the outer domains are linked by the 2 following boundary conditions.

The continuity of the normal velocity at the interface at r = R:

U e R = U i R = 0, (3.23)
And the continuity of the tangential velocity and viscous stress along the polar, θ, and the azimuthal, ψ, directions:

U e θ = U i θ ,

Validation tests

Numerical results on the drag coefficient obtained from axisymmetric simulations have been compared with the literature data (Tab. 3.3). The maximum relative difference observed for the Re = 1, µ * = 1 and ρ * = 1 case is equal to 4.27 10 -2 . Note that the simulations at the smallest Reynolds number were carried out with small-time steps (0.4 × 10 -4 for Re = 1 compared to 0.4 × 10 -3 for Re = 250 at µ * = 0.5 and ρ * = 1) to ensure convergence and mitigate numerical instabilities arising from the calculation of the viscous terms at the center of the droplet center with the polar mesh. However, the relative error for these cases remains very low, and we consider that satisfactory agreement is achieved between our simulations and literature results.

Simulation results were also compared to the expression proposed by Hadamard and Rybczynski, Eq. (3.26), Hadamard [START_REF] Hadamard | Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux[END_REF] to relate the drag coefficient of a translating droplet to that of a solid sphere and a bubble of the same size and velocity. Here we used Eq. ( [START_REF] Zhi | Drag coefficients of viscous spheres at intermediate and high reynolds numbers[END_REF] (parenthesis), and Brackets: Oliver and Chung [START_REF] Oliver | Flow about a fluid sphere at low to moderate reynolds numbers[END_REF]. The differences between our results and those of Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF] obtained with the same code, can be explained by the differences in mesh parameters and steady-state criteria. (3.27)

Cd(Re, µ * = ∞) = 24 Re 1 + 0.15Re 0.687 (3.28)
Finally, the non-influence of the outer mesh on the drag coefficient and onset of 3D bifurcation was verified. Simulations with a fully polar mesh (inner and outer) were carried out their results compared with those obtained on the "hybrid" mesh Fig. 3.5b, with the computational domain discretized with an inner polar mesh and an outer mesh based on the streamlines. The results, presented in Tab. 3.4, are very similar. However, as the polar mesh revealed to be less computationally efficient (the computation time required to reach steady state was longer), we have chosen to keep the "hybrid" mesh for this study. 3.4: 2D and 3D drag coefficients, and amplification factor (see Sec. 3.6.3). Comparison of results obtained from simulations achieved with the "full polar" and the hybrid "polar and streamlines" meshes.

Outline of the parametric study

The range of Reynolds number, viscosity, and density ratios considered in the numerical investigation is the following: [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF][START_REF] Davies | The mechanics of large bubbles rising through extended liquids and through liquids in tubes[END_REF][START_REF] Chang | Unsteady flow about a sphere at low to moderate reynolds number. part 1. oscillatory motion[END_REF][START_REF] Kleinman | Unsteady conjugate mass transfer between a single droplet and an ambient flow with external chemical reaction[END_REF]200,500] • µ * ∈ [0.05, 0.1, 0.2, 0.5, 1, 2, 5, 10, 20]

• Re ∈ [1,
• ρ * ∈ [0.1, 0.2, 1, 5, 10]
These values and their combinations ensure that the droplet sphericity condition is met (see Sec. [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF]

.3.2).

Although the density ratio has a small impact on the drag coefficient Feng and Michaelides [START_REF] Zhi | Drag coefficients of viscous spheres at intermediate and high reynolds numbers[END_REF], Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF] in 2D, we chose to explore a sufficiently wide range of ρ * to check whether this conclusion applies to 3D bifurcations.

To limit computational time, the following strategy was adopted. First, we run each case of the parametric study in axisymmetric conditions until steady-state is reached. The resulting flow field is then used as the initial condition for the 3D simulations. A weak random disturbance (typically 10 -5 U ) is added to the initial condition. If this initial perturbation vanishes during the calculation, the case refers to a stable axisymmetric flow. On the other hand, the amplification of the perturbation reveals the development of a 3D regime. To discuss the origin of the instability, simulations were also carried out by forcing a 2D flow inside or outside the droplet.

Experimental observations

Droplet trajectories

Depending on the droplet size and on the fluid system considered, different droplet trajectories are observed. Some droplets follow a straight vertical trajectory and show no 3D wake characteristics (see Fig. 3.7a and 3.7b), while other droplets follow a zigzag trajectory, as shown in Fig. 3.7c. The zigzag ascent resembles that of ellipsoidal bubbles, although we have verified that our experiments involved only spherical objects, as confirmed by their degree of sphericity. Moreover, the corresponding evolutions of the droplet's aspect ratio χ, reported respectively in Fig. 3.7d 3.7e and 3.7f, are very close to one. For comparison, the bubble's wake instability reported by Magnaudet and Mougin [START_REF] Magnaudet | Wake instability of a fixed spheroidal bubble[END_REF] is observed for χ > 2.21. The observed zigzag trajectories are therefore probably the consequence of hydrodynamic instabilities giving rise to a 3D flow.

The amplitude of the zigzag trajectories was determined from the images. A maximum of three oscillations were observed in the camera field. It should be noted that in some cases, there appears to be no complete oscillation in the plane of vision, making it difficult to assess whether the bifurcation oscillates or is steady. In the explored operating conditions, the amplitude of the observed oscillations ranges between 6.59×10 -5 and, 3.84×10 -3 . It is difficult to unambiguously correlate the weakest oscillations with a 3D flow. On the other hand, the upper limit is closed to the case shown in Fig. 3.7c, where the droplet of diameter 2R = 2.91 mm undergoes oscillations with amplitude 2.56 mm, of the same order of magnitude of its size. More generally, in the experiments carried out, the amplitude of the oscillations observed ranged between 3.15 × 10 -2 R and 2.34R, confirming that it can be significant and measurable when processing the trajectory of a droplet.

The velocity magnitude U and velocity components U x and U y of the droplets shown in Fig. 3.7 are reported in Fig. 3.8. The droplets with the straight trajectories exhibit constant velocity, both in magnitude and in the plane of the camera (Fig. 3.8a and 3.8d, and Fig. 3.8b and 3.8e), whereas in the case of the droplet with a zigzag trajectory, velocity fluctuations are observed, that are consistent with the 3D flow features (Fig. 3.7c and 3.8f). The oscillation frequencies are different between x and y direction, this can be explained by the fact that the x-velocity is maximal when the y-velocity is equal to zero and minimal when the y-velocity is at an extremum, and therefore the x-velocity oscillation has a doubled frequency compared to the y-velocity oscillation.

Experimental flow map of rising droplets

All the experimental results are reported in the form of a flow regime diagram (µ * , Re) in Fig. 3.9. Droplets that achieved a straight trajectory are represented with gray symbols, while those evidencing oscillating trajectories are shown using black symbols. Results presented in Fig. 3.9, and literature on solid sphere unsteady bifurcation, suggest that a freely rising droplet could experience a similar bifurcation as the one observed for solid spheres. The origin of this transition will be investigated in the following, thanks to numerical simulations. In particular, the role of the flow inside the droplet will be emphasized in the following.

Drag coefficient

The drag coefficients deduced from experimental droplet trajectories already considered in Figs. 3.7 and 3.8 are compared in Fig. 3.10 with the evolution of C d deduced from axisymmetric simulations, at the same viscosity and density ratios. We can observe that at high Re, i.e. in Figs. [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF].10b and 3.10c, the experimental values deviate from the 2D predictions.

The average drag coefficients, C d deduced from experiments with Heptane and n-Pentane, are plotted in Fig. 3.11 and compared with the drag coefficient derived from the axisymmetric simulation. Only the droplets satisfying the sphericity criteria are represented (more than 50 cases for each fluid). We can clearly see that for Reynolds numbers above 93 (n-Pentane case) and 156 (Heptane case), the experimental drag coefficient are significantly higher than those predicted by assuming axisymmetric flow.

For example, at Re = 350, the deviation is +38.7% for n-pentane and +6.31% for Heptane.

These experimental observations provide clear evidences of the development of non- 84CHAPTER 3. 3D BIFURCATIONS STUDY axisymmetric flows resulting in oscillatory trajectories and a subsequent significant increase in the drag coefficient.

DNS study of the transition from 2D to 3D flow

The transition from 2D to 3D flow has been studied using DNS. Numerical simulation enables us to visualize and study the different flow fields reached outside and inside the drop, and to define the origin of the bifurcation.

Flow description

Two-dimensional axisymmetric flow

When

Steady 3D flows

For certain values of Reynolds number, viscosity and density ratios, a bifurcation from axisymmetric flow to stable non-axisymmetric flow is observed. Here, we analyze the velocity field and streamlines of three different configurations undergoing this stable bifurcation:

• The instability comes from the increase in Reynolds number that should reduce the pressure at the droplet center until a critical Reynolds number is reached, the main driving force for the flow becomes the inertia term (as viscous and pressure force get negligible), providing the condition for the development of the azimuthal velocity disturbance (O(10 -5 )U ) into a 3D flow.

When instability occurs, two phenomena are observed: the appearance of an azimuthal velocity distribution and the displacement of internal vortex cores. The latter resemble a twisted torus (see Fig. 3.15a) or even a completely deformed one, when increasing Re (see Fig. 3.15c. Indeed, in case 3, the vortex core moves towards the back of the drop, leading to the formation of four other vortex cores (one in each quarter), each coupled to the outside via the interface boundary conditions. Despite the absence of detachment (mu * < 2), the displacement of the vortex core creates a new stagnation point at the interface, leading to the formation of a vortex in the continuous phase. As a result, if the flow field is observed in different planes, the shape of the flow can be entirely different. This can be seen in Figs. 3.14a, 3.14d, 3.14c and 3.14f. In the first raw (plane xy), the streamlines are similar to the Hill vortex, while in the second raw (plane xz), the streamlines are radically different. The heart of the vortex is stuck to the rear part of the interface, and numerous streamlines follow a straight trajectory from the rear stagnation point of the droplet, towards the front of the droplet.

For a viscous droplet (Case 2, µ * = 5), the bifurcation is mainly visible in the outer flow (see Figs. [START_REF] Mei | A note on the history force on a spherical bubble at finite Reynolds number[END_REF].14b and 3.14b). The threshold Reynolds number is around 300. The axisymmetry of the flow breaks down into two opposing vortices rotating around the x axis (see Fig. 3.15b) and the wake flow undergoes competition between centrifugal acceleration and a local pressure minimum inside its core. When the wake vortex shifts from the dominance of viscous effects to that of pressure gradient effects, the second flow regime occurs. The structure of the wake vortex structure is best understood in this regime if we focus on vortex motions around local pressure minima, as proposed by Jinhee Jeong [START_REF] Hussain | On the identification of a vortex[END_REF]. The typical flow generated is clearly visible in the yz plane, where he streamlines reveal the symmetry breaking behind the drop, but we can hardly see any 3D structure inside the drop in this case.

3D unsteady flows

If the Reynolds number increases further, the stable 3D flow can evolve into an unstable regime. We propose here to analyze the velocity field and streamlines of two different cases undergoing this unstable bifurcation:

• µ * = 0.05, ρ * = 1 and Re = 50

• µ * = 10, ρ * = 0.2 and Re = 500

In the first case, corresponding to µ * < 2, the development of the 3D unsteady flow is mainly observed inside the droplet. We can see in Fig. 3.16a, that the two counterrotating vortices grow and shrink alternately due to an oscillating motion of the two vortices. The oscillation maintains a certain symmetry in the alternating movement of the vortices. If the viscosity ratio and Reynolds number are higher, the flow is unstable. As in the previous example, we can observe swirling movements around local pressure minima in the vortex core. The onset of azimuthal velocity in the axisymmetric bifurcating vortex is no longer manifested by two equal and opposite sections, but we observe instead one section growing faster behind the droplet and eventually dominating the other, which disappears in the vicinity of the droplet, while at a greater distance from the interface, the opposite phenomenon is observed. (Figs 3.16b 3.16d) This imbalance between the two opposing parts of the recirculation behind the drop is particularly visible in the xy plane (Fig. 3.16b). The vortex core is no longer contained within the yz axis, but stretches along x. The strongest vortex, close to the drop, takes up most of the space behind the interface, while the other is pushed further back. Eventually, this position returns to the initial equilibrium and the same process repeats itself. In the xz plane, the sequence of steps is similar.

In the yz plane at x = 2 (Fig. 3.16f), the flow oscillates around a mean axis (pointing to the left of the image). We can see that there are two main regions in the flow: one where the flow detaches from the droplet, the other where the vortex grows to take most of the place in the wake of the droplet. In this 2D projection, the flow is colored by the planar velocity M ag(U y , U z ), and we can observe that there is only a small local part of the flow that acquires most of the velocity.

Internal or external origin of instability

As indicated previously in Sec. 3.4.3, the simulations are first carried out in 2D until reaching a steady state, which is then used to initialize the 3D simulations, inside and outside the drop. In cases where 3D flow develops, 2 additional simulations were performed to determine the internal or external origin of the bifurcation. Each simulation was restarted, from the converged 2D solution, by freezing either the internal flow or the external flow.

In the example of Fig. 3.18 (with µ * = 5, ρ * = 1 and Re = 350), we observe that the flow remains axisymmetric when the external flow is frozen, while it becomes 3D when a 2D flow is imposed inside the drop. We can therefore conclude in this case that the origin of the bifurcation comes from the external phase. Note in Fig. 3.18h that the location of the plane of symmetry selected by the flow is slightly different from that observed in the case where no constraint is imposed on the internal flow. Fig. 3.18i.

In the example of Fig. 3.19 (with µ * = 0.2, ρ * = 1 and Re = 100), the flow remains axisymmetric when the internal flow is frozen, while it becomes 3D when the external This procedure was repeated for the more than 50 cases of simulations leading to 3D flow. In each case, we observed that the bifurcation originated either inside or outside the droplet, but never simultaneously inside and outside. We were also able to identify that, regardless of the value of the density ratio ρ * , the origin of the bifurcation is always internal when µ * < 2 and external otherwise (µ * >2). This information is reported in the flow phase diagram presented in the following section.

Flow phase diagram from numerical simulations

The determination of the transition threshold between stable and unstable flow regimes has been conducted by considering the amplification rate of a perturbation.

The amplification of computer round-off errors due to numerical discretization leading to the bifurcation of the flow can be a very time-consuming. To reduce the computational time of the 3D simulations and control the level of the disturbance, we introduce a 92CHAPTER 3. 3D BIFURCATIONS STUDY weak flow perturbation during the computation at a given time step, which ensures that the 2D regime is converged. The velocity perturbation consists in a random azimuthal velocity distribution, with maximum amplitude ±10 -5 U which is added either inside or outside the droplet when the outer or inner flow is frozen, respectively.

In some cases, a 3D bifurcation can develop with two orthogonal symmetry planes. Therefore, a validation test was made to verify that the same simulation configuration as explained in Sec 3.4.3 with a random perturbation yields the same flow features. As observed in Fig. 3.20, the same 3D flow should have symmetry with different angles to y and z to validate the mesh independence. In Fig. 3.20, we can verify this requirement with two identical meshes and simulation configurations, we obtain the same flow state but with a different orientation. The maximum azimuthal velocity W max (either inside or outside the droplet) is a metric that allows us to detect the bifurcation of the 3D flow. The linear amplification of the disturbance can be plotted in a semi-log graph (time, log(W max ) see Fig. 3.21).

We then determine α, the corresponding rate of damping (resp. amplification of the perturbation), as illustrated in Fig. 3.21. If α is negative, we observe a damping of the perturbation: the flow evolves towards the axisymmetric case, which corresponds to the Hill's vortex for these conditions of Reynolds number and viscosity ratio. If α is positive, the azimuthal velocity increases, until it saturates to a plateau value corresponding to the emergence of the 3D regime.

For all simulations, we gather the amplification/ damping factors. For each viscosity and density ratio, we can then determine by interpolation the transition threshold from positive to negative α as a critical Reynolds number (corresponding to the marginal state). Furthermore, for each Reynolds number, we can interpolate a critical viscosity or density ratio for which the amplification factor is equal to zero. We have also added simulations to refine the threshold determination. Three interpolation examples are presented in Fig. 3.22, demonstrating that the Reynolds number has an important effect on the amplification factor. However, it appears that the effect of density ratio is different for internal and external bifurcations, weak effect on external bifurcation and significant effect when the flow instability originates from inside the droplet, see Fig. For the specific case ρ * = 1, we have increased up to 174 simulations the initial 49 simulations defined by our parametric study, to accumulate enough data for a precise determination of the Reynolds number threshold. All the data from post-processing of flow and amplification factor at ρ * = 1, are reported in Fig. 3.23. Different colors are used for axisymmetric flow, 3D steady (internal or external origin) as well as 3D unsteady (internal or external origin) flows.

We observe, in Fig. 3.23, that all the simulations for µ * < 2 present an internal 3D bifurcation. On the other hand, the simulations for (µ * > 2) yield external 3D bifurcation. It is important to note that for all the simulations conducted, none of them have shown simultaneous internal and external 3D instability onset, and the marginal curve between internal and external flow bifurcation is located at µ * = 2. This frontier is similarly observed in the case of droplet axisymmetric steady wake in Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF], and Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF]. It is a peculiar outcome of this work which deserves further study the droplet flow stability in function of µ * , to understand why µ * = 2 appears to be stable whatever the Reynolds number. The critical Reynolds number between axisymmetric and 3D steady flow for a fixed solid sphere (corresponding to µ * → ∞) is Re c = 212. The external bifurcation marginal states we determined tend to this value when the viscosity ratio becomes very large. The critical Reynolds number for a solid sphere unsteady bifurcation is Re c,unsteady ≈ 270 and is verified for large viscosity ratios as well in our simulations. The flow structure before the onset of flow instability is composed of detached streamlines forming a steady axisymmetric wake. We can relate Eq. (3.29) for the evolution of the separation angle of the flow to the observed bifurcation marginal curve. Those two features of the flow (appearance of the recirculation and appearance of 3D bifurcation) seems to be closely related when the viscosity ratio and the Reynolds number are varied (and almost no effect of density ratio). From Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], the evolution of the separation angle depending on the Reynolds number is given by

θ d =
Re c (µ * ) = 212e 0.1 1+µ * coef 1µ * coef 2 (3.31)
We use the curve fit algorithm from the python library scipy.optimize to fit Eq. (3.31) on our external bifurcation critical points. The result is shown in Fig. 3.24 and Eq. (3.32), the correlation fits really well with our data with a model max error = 24.9, a model mean absolute error = 5.80, a model R2 score = 0.995 and a model mean absolute percentage error = 0.0272%.

Re c (µ * ) = 212e 0.1 1+µ * 0.132µ * 13.7 (3.32)
It is considerably more challenging to obtain meaningful data on the critical points between steady and unsteady bifurcation, as the interpolation between amplification factor cannot be used. But for what concerns external bifurcation, we know that at µ * → ∞, the steady bifurcation critical point is Re c = 212 for steady bifurcation and Re c ≈ 270 for the unsteady bifurcation. In our simulation, this ratio 270 212 = 1.27, is verified: every simulation 1.27 times over the critical bifurcation Reynolds number is unsteady for all µ * . Therefore, an approximation of the unsteady bifurcation critical point could be given by Eq. (3.33). In Fig. 3.24b, we can observe that for the external bifurcation, the impact of density ratio is negligible because the density ratio is not involved in the definition of the Reynolds number (external). Therefore, we can conclude that the key parameters to study the external bifurcation are only the Reynolds number and the viscosity ratio.

The left part of the Fig. 3.23 (internal 3D bifurcations), shows a similar trend to the left part when µ * → 2 but does not to converge to an asymptotic value when µ * → 0, see Fig. 3.25. This is very intriguing because at µ → 0, we tend towards the case of a bubble in a liquid, and we do not expect any 3D bifurcation for this specific case. Nevertheless, to correctly reproduce the case of bubble, we should also include the effect of the density ratio. When we study internal bifurcation, we should look at the quantity of interest of the concerned phase, being the internal Reynolds number. Fig. 3.26 shows the internal bifurcation map with this variable, and it appears that using Re i and µ * to characterize the internal bifurcation is coherent, accounting for the effect of ρ * in the internal Reynolds number Re i , see 3.26. Similarly to the external bifurcation, it is challenging to obtain meaningful data on the critical points between steady and unsteady bifurcation, as the interpolation on amplification factor cannot be used. We have seen for external bifurcations that the ratio between the unsteady and steady critical bifurcation point is constant equal to 1.27, so we observe that a similar approach with a ratio between steady and unsteady bifurcation also works in the case of internal bifurcations but with a constant around 2.1, as it is shown in Eq. (3.37) and (3.38). In Fig. 3.22b, we show that the density ratio as a strong impact on the amplification factor for the internal bifurcation (on the contrary of Fig. 3.22c where the impact of ρ * is weak for external bifurcation). It is therefore straightforward to use the internal Reynolds number (equal to Reρ * /µ * ) which characterizes the strength of inertial effects inside the droplet as the dimensionless number for the internal bifurcation.

Compiling those results for different ρ * , we can place the data points from our simulations in the previous internal bifurcation graph, see Fig.

Effect of 3D flow on drag

The overall flow dynamics is modified when the wake structure or the internal vortex flow experiences a hydrodynamic transition from steady 2D, to steady 3D and eventually unsteady 3D. We consider first the consequences on the drag coefficient. In Fig. 3.27 we can see the effect of every type of bifurcation previously studied on the drag coefficient.

DNS STUDY OF THE TRANSITION FROM 2D TO 3D FLOW 101

As observed, the bifurcation, regardless of its type, is responsible for a drag coefficient increase compared to the equivalent axisymmetric case.

This increase is quantified with a drag enhancement factor. For all µ * we have investigated, the results on 3D drag coefficients divided by the corresponding axisymmetric drag coefficient is calculated. The ratio is reported as a function of the external Reynolds number Re in Fig. 3.28. We can clearly distinguish the axisymmetric case ( Cd 3D Cd axi = 1), from the case with 3D bifurcation ( Cd 3D

Cd axi ̸ = 1). More interesting is that we can also observe a difference between the steady bifurcation and unsteady bifurcation. The steady bifurcation drag ratio align, until we reach the Reynolds number at which the droplet experiences a 3D bifurcation, and then the drag ratio value drops. For the steady bifurcation, the drag ratio can be fitted by a linear function, which intercept and slope are reported in Tab. From this table and some linear interpolations, we deduce the following drag coefficient enhancement factor (valid for ρ * = 1) for internal bifurcation, see Eq. (3.41) and for external bifurcation, see Eq. (3.42). Those correlations have been tested on some simulation data of steady bifurcation cases, we obtained the following error metrics:

• Max error: 0.181, the maximum difference between predicted and simulated C d (3D),

• Mean absolute error (MAE): 0.159, the average of the differences measured between predicted and simulated values,

• R2 score: 0.628, measures the proportion of the variance in the simulation data that can be explained by the prediction,

• Mean absolute percentage error (MAPE): 0.111%, similar to MAE but based on the relative differences, i.e., divided by the value from DNS.

Cd 3D

Cd axi (µ * < 2) = 0.235 -0.258 ln(µ * ) + (0.00411 + 0.000550 ln(µ * ))Re (3.41)

Cd 3D Cd axi (µ * > 2) = 1.06 -0.0757 ln(µ * ) + 0.000411Re (3.42)
If we aim at a comprehensive model to predict the drag coefficient, it is necessary to consider the effect of density ratio on the drag coefficient of 3D cases. For axisymmetric cases, we know that the impact of the density ratio stands within 5% order of magnitude from the axisymmetric case at ρ * = 1. This observation is verified when doing 3D simulations of axisymmetric cases. In Fig. 3.30b, we compare the evolutions in function of internal Reynolds number for internal bifurcation case and external Reynolds number for external bifurcation cases. From Fig. 3.30a, we can see that the density ratio as a weak impact on the external drag coefficient even after the bifurcation. To obtain the drag coefficient for a case with external bifurcation, it is sufficient to compute the corresponding axisymmetric one and utilize Eq. (3.42) to determine the corresponding 3D drag coefficient (for steady bifurcation).

From Fig. 3.30b, we can see that the bifurcation behavior is identical with an equivalent internal Reynolds number, but the external Reynolds number being different, the drag coefficient is also different. The change in drag coefficient with the 3D bifurcation respects the same logic as in (3.41) in regard to Re i : Cd 3D Cd axi (ρ * = 0.628) = 0.123 and We are comparing results from freely moving droplets (experiments) and fixed droplets in a uniform flow (numerical simulations).

Bifurcation diagram from simulations and experiments

Placing the bifurcation criteria from Eq. (3.36) (steady) and (3.38) (unsteady), in the bifurcation graph from Sec. 3.5.2, we obtain Fig. 3.32. In this figure, the experimental data are separated depending on their trajectory (straight or oscillatory). We have previously commented that this figure suggests a similar behaviour between a freely rising droplet and solid spheres. This finding is remarkable considering that we do not expect any external bifurcation, similar to that observed for solid spheres, for µ * < 2. The results shown in Fig. 3.32 confirm this particularity, the freely rising droplet trajectory doesn't oscillate only when the internal bifurcation appears, so the fact that the droplet is freely rising might introduce a phenomenon that we couldn't observe in the fixed droplet simulation. We use the value of drag coefficient as an identification criterion for the bifurcation, so we will identify both steady and unsteady bifurcations from experiments. We put the experimental data in the bifurcation graph together with the bifurcation 106CHAPTER 3. 3D BIFURCATIONS STUDY criterion from Eq. (3.36) (steady) and (3.38) (unsteady), in Fig. 3.33. To identify an experimental point as axisymmetric or 3D, we use the closest drag coefficient simulation data between axisymmetric simulation and 3D simulation (with bifurcation). For Reynolds number below the critical bifurcation point, the 3D simulation is equal to the axisymmetric simulation result, so we had to conclude on axisymmetry for those experimental conditions. However, it is interesting that as soon as the Reynolds number passes the bifurcation criterion, the axisymmetric result is no more the closest for most cases. Rare exceptions are identified in high Reynolds case where the experimental uncertainty and small deformation make possible for the axisymmetric and unsteady drag coefficients to be the actual solution. This figure also explains why we were unable to observe 3D bifurcations for Cyclohexane droplets in experiments. The drag coefficient appears more suited to identify the bifurcation than the actual trajectory, as it better correlates the simulation observation. So the phenomenon that makes a discrepancy between the simulation bifurcation criterion and the experimental path does not have a significant impact on the drag coefficient.

Drag coefficient

We also use experiments to discuss the relevance of the significant impact of the 3D flow on the drag coefficient observed in the numerical simulations. We have determined . This clearly validates the axisymmetric simulation as long as the flow remains axisymmetric. More interestingly, the experimental results follow the 3D simulation results after the critical Reynolds number, confirming the significant difference in the drag force when the flow becomes 3D. The transition from 2D to 3D is also confirmed when looking at the drag coefficient evolution. For applications, we recommend that the specific drag coefficient expression proposed in this work to be used to describe accurately droplet dynamics for Reynolds numbers larger than the critical Reynolds number determined in this work.

Note that Fig. 3.34a also indicates that for high Reynolds number cases, a drag coefficient larger than the 3D simulation is observed in experiments. This can be explained by a slight deformation of the droplet in the experiments, resulting in a larger drag force.

Strouhal number

The configuration of a fixed sphere for simulations is expected to have a significant impact on unsteady cases when compared to experimental data obtained from freely rising droplets at similar Reynolds number, viscosity, and density ratios. The motion of the droplet due to the lift force induced by the 3D flow should modify the external 109 sphere study. However, the results for Heptane are not sufficiently numerous to use for correlations. The results with n-Pentane droplets are more numerous, and we can see that, for small Reynolds number, the experimental outcomes align correctly and increase with Re before becoming more scattered.

Comparison with Kiya et al. [START_REF] Kiya | Near-wake instabilities and vortex structures of three-dimensional bluff bodies: A review[END_REF] solid sphere results in Fig. 3.35c suggests that our experimental results are in the same order of magnitude than results for solid sphere bifurcation. So, major droplet bifurcation mechanisms could be similar to those occurring for solid spheres, thus explaining Fig. 3.32 result's. For N-pentane, in Fig. 3.36a, the simulation results largely underestimate the experimental ones, for some cases with a factor of hundred. This indicates that the vortices in simulation are less frequent and weaker than in experiments. Nevertheless, for 200 < Re < 300, the slope is identical, depicting that despite the factor difference, the mechanisms behind the unsteady bifurcation remain the same. For simulation cases with Re > 300, the simulation outcomes appear to be chaotic. The data visualization indicates that there are two primary frequencies, and the post-processing Fourier analysis 110CHAPTER 3. 3D BIFURCATIONS STUDY yields a similar modulus for both of them. Consequently, the calculation of Strouhal number is performed with the frequency of the largest modulus, resulting in a certain scarcity in our results. Experimental data appear to be spared by this phenomenon, and only one frequency is clearly identified globally. This might be caused by the uncertainty in the post-processing, blurring the small displacement generated by the higher frequency, or the difference is due by the freely moving drop.

Results for Heptane, in Fig. 3.36b, do not have as much experimental data as N-Pentane. Nevertheless, the simulation data have the same order of magnitude as the experimental ones, with this experimental data being generally smaller than the ones for N-Pentane.

The comparison with Kiya et al. [START_REF] Kiya | Near-wake instabilities and vortex structures of three-dimensional bluff bodies: A review[END_REF]'s figure referenced in Fig. 3.36c indicates a similar evolution for transition between laminar and turbulent flow around solid spheres. Therefore, the bifurcation being internal rather than external demonstrates similar mechanisms. A plateau value for the Strouhal number evolution with Reynolds number is verified for experimental data. For simulations, we suppose that this transition regime is influenced by the configuration of a fixed sphere.

Statistical analysis

We have conducted a thorough investigation of the drag coefficient obtained through three distinct methods, namely experimental, simulation with axisymmetric assumption, and 3D simulation. Our results, show that the axisymmetric assumption is not verified in certain conditions. But the most important remain to quantify the degree of accuracy of the 3D simulation and the deducted bifurcation criterion with the experimental results, knowing that other phenomenon could corrupt the experimental results (deformation, contamination). We propose to perform a t-test with pooled standard error to consider the experimental and simulation uncertainty in the probability of a model representing the experimental results.

The pooled standard error, denoted as sp, see Eq. (3.43) is calculated by combining the standard deviations of the two groups (by example experimental results and corresponding axisymmetric simulation), weighted by their respective sample sizes and uncertainty. This pooled estimate of the standard error provides a more accurate measure of the overall variability. The p-value indicates that the two groups have an important probability of being accurate descriptions of each other, while the t-statistic quantifies the two group mean gap. We propose to verify three hypotheses:

• Hypothesis zero (H 0 ): All experimental drag coefficient results could be represented by corresponding axisymmetric simulations drag coefficient What we can conclude from this analysis is that the axisymmetric assumption couldn't be considered valid for our experiment. The axisymmetric assumption is verified for cases under critical Reynolds number, and leads to large underestimation overwise, therefore validating the critical bifurcation point Eq. (3.40). The 3D simulation explains the experimental data for case over the critical Reynolds number.

Conclusion

We have studied the 3D bifurcation inside and outside a spherical droplet under uniform flow condition, experimentally and numerically. The experiments indicate that 3D 112CHAPTER 3. 3D BIFURCATIONS STUDY bifurcation can be observed, while also verifying the spherical droplet assumption. Numerically, we intend to conduct a thorough parametric investigation to gain an in-depth understanding of this phenomenon. We also discuss the impact of the fixed droplet in simulations compared to the freely rising droplet in experiments.

We first identify that bifurcation could be experimentally observed through a zigzagging droplet path. From simulation, we identify two kinds of bifurcation that could be observed: the external and internal bifurcations. It means that the bifurcation origin is always situated only in one phase, despite that, the flow in the other phase can be largely impacted by the bifurcation due to boundary conditions at the interface. But, we couldn't find any case for which the analysis conducted by freezing one phase flow yields a bifurcation to appear on both sides of the interface. No matter where the bifurcation originated, we found that depending on the Reynolds number, there can be a steady or unsteady bifurcation. The mechanism can be compared, but the effect on drag and lift is significantly different.

Internal bifurcations appear in the center of Hill's vortex. As the Reynolds number increases, the pressure in the center of the droplet decreases. This is because the internal Reynolds number becomes larger, and the viscous force becomes negligible compared to the inertia force. The driving force for the bifurcation is the inertia force, and a small disturbance in the azimuthal velocity or pressure can initiate the bifurcation to a 3D flow.

If we increase the Reynolds number even more, the previous steady three-dimensional flow is no longer stable, resulting in an unsteady bifurcation.

External bifurcation similarly appears when the wake vortex experiences a competition between centrifugal acceleration and a local pressure minimum within its core. The flow breaks down, into two steady counter-rotating vortices creating an important lift on the droplet, but a smaller drag increase than the internal bifurcation. With an increase in the Reynolds number, it was revealed that the wake is composed of unsteady periodic hairpin vortices with constant orientation, as well as new oppositely oriented hairpin vortices induced by the interaction between the near-wake flow and the outer flow.

From the analysis of the amplification factor of some small disturbance inserted in the flow, we were able to precisely interpolate the critical bifurcation threshold of the steady bifurcation, both internal and external. By utilizing the correlation for axisymmetric flow separation, we were able to precisely fit this expression and accurately predict the critical steady bifurcation threshold. The unsteady bifurcation criterion was not correlated with the same accuracy because the method with the amplification factor is no longer possible. Thus, we observed that a ratio of 1.24 between unsteady and steady external bifurcation critical point and of 2.1 between unsteady and steady internal bifurcation critical point is a good approximation.

One important part of this investigation is the impact on the drag coefficient of the 3D bifurcation. Our findings indicate that for steady bifurcation, the ratio of 3D drag coefficient to axisymmetric drag coefficient increases linearly with the Reynolds number for every µ * and ρ * . Therefore, we can now predict the drag coefficient for steady bifurcation cases with only the drag coefficient corresponding to axisymmetric behavior, the bifurcation criterion and the correlation for this drag coefficient enhancement factor.

Finally, we had a discussion of experimental and simulation results together. The assumption of a fixed sphere in simulation creates a deviation for unsteady bifurcation, resulting in a different Strouhal number than the experimental one. Nevertheless, we obtained the same bifurcation diagram and similar drag coefficient results in experiments and simulations. So, we conducted a statistical analysis that gives a quantitative metric that invalidates the axisymmetric assumption and validates the 3D simulation results, all while verifying to be over the bifurcation criterion.

With all this information, we can predict for which parameters the bifurcation will appear, what the consequence on drag. We can easily find the correct drag coefficient value based on a simple axisymmetric simulation result. The experimental study evidences that those phenomena could be observed under real-world conditions. We believe that those 3D flows could be a real asset to liquid-liquid extraction process enhancement. And we identify an interesting outlook on the subject of flow stability around µ * = 2. This chapter is an under review article, see Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]. Based on the observation that most studies of heat or mass transfer in dispersed two phase flows are limited to simplified or specific operating conditions, we used DNS to address the needs for more robust correlations, arising from the development of new liquid-liquid solvent extraction processes, involving e.g., more viscous solvents or different hydrodynamic conditions.

A large parametric study of mass transfer from an isolated droplet was achieved to propose a physically based correlation of the apparent Sherwood number, Sh, that is sufficiently generic to apply to emerging and future applications. The updated correlations derive first from a sound analysis of the features, in terms of concentration fields inside and outside the droplet and interface properties, to challenge the validity of the assumptions of mass transfer correlations.

Then, we propose proper modifications in the formulation of existing models. In a second step, the accuracy and the robustness of the new models have been assessed by comparison with Sh obtained from DNS. The proposed correlations are applicable for mass transfer as well as for heat transfer problems, in most configurations encountered in chemical plants.

Introduction

Solvent extraction is a major separation process in the chemical, mining and pharmaceutical industries, for selectively extracting substances from complex mixtures. This operation is based on the ability of a solvent to specifically dissolve chosen target compounds, thus allowing their separation from the other constituents. In conditions typical of liquid-liquid extraction, one phase (in general the solvent) is dispersed in the other non-miscible one to enhance their contact area. Strong agitation prevails to generate small droplets, and therefrom to maximize the surface of contact between the two immiscible fluids. The mass flux of a solute transferred from one phase (in kg.s -1 ), with a concentration C i (in kg.m -3 ), to a second phase, with a concentration C e (in kg.m -3 ), is expressed by Eq. (4.1):

Φ = hS(C i -kC ∞ e ), (4.1)
with S the contact area (in m 2 ), and h the mass transfer rate (in m.s -1 ).

Beyond equilibrium principles, understanding the underlying mass transfer mechanisms and their kinetics is essential to predict and improve the process efficiency. The apparent (or global) mass transfer coefficient h is usually expressed by its dimensionless form, the Sherwood number Sh, for which numerous correlations, mainly empirical, exist. Most of them have been reviewed by Kumar and Hartland [81]. Sh compares the flux of solute Three mass transfer regimes are generally distinguished, according to a criterion proposed by Wylock et al. [START_REF] Wylock | Gas adsorption into a spherical liquid droplet: numerical and theoretical study[END_REF] and based on the diffusivity ratio of the solute between the inner and the external phases, D * = D i /D e , and on the concentration ratio achieved at thermodynamic equilibrium across the interface, k (see definition in Eq. (4.14)).

• The isolated internal problem, for k √ D * << 1, where the main mass transfer resistance is located inside the droplet (Sh = Sh I ). In this configuration, the solute concentration is uniform along the interface and fixed by the external conditions Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF], Grassia and Ubal [START_REF] Grassia | Streamline-averaged mass transfer in a circulating drop[END_REF].

• The isolated external problem, for k √ D * >> 1, where the main mass transfer resistance is located in the continuous phase (Sh = Sh E ), resulting in a uniform concentration inside the droplet and along the interface Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF], Juncu [START_REF] Gh | Conjugate heat and mass transfer from a solid sphere in the presence of a nonisothermal chemical reaction[END_REF].

• When none of the transfer resistances can be neglected, the problem is said "conjugated" or "coupled".

The flux conservation of the solute across the interface (with no accumulation) results in the following relation:

1 Sh = 1 Sh i + kD * Sh e (4. 3
)
where Sh i = h i d D i and Sh e = hed De are the internal and external Sherwood numbers (h i and h e being the internal and external mass transfer coefficients), to be distinguished from the values Sh I and Sh E when considering only the internal and external problems, respectively. In the coupled problem, the apparent mass transfer coefficient is classically deduced from the double-film model (Fig. 4.1), where the mass transfer process is seen as the series connection of the two (finite and non-zero) resistances located on each side of the interface, generally assuming Sh i = Sh I and Sh e = Sh E . The main correlations of mass-transfer from an isolate droplet have been reviewed and discussed recently by Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], based on a comprehensive hydrodynamic and mass-transfer study by Direct Numerical Simulation (DNS). The authors achieved a parametric study covering the three regimes, with particular emphasis on the conjugate mass-transfer regime, investigating the influence of hydrodynamic conditions on both the solute spatial distribution and the transfer evolution. They evidenced that in most cases, the internal (resp. external) mass-transfer coefficient, Sh i , (resp. Sh e ), deviates from the corresponding coefficient typical of the isolated internal (resp. external) problem, Sh I (resp. Sh E ), and that in most conjugated problems these differences between Sh i and Sh I (resp. Sh e and Sh E ) result from the non-uniform distribution of the solute prevailing in or/and around the droplet. Consequently, the use of the coefficients of the asymptotic cases in Eq. ( 4.3) to assess the global coefficient Sh typical of a coupled problem, violates the founding assumptions of the double-film model, and can lead to errors greater than 30% Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. However, a current problem encountered in designing emerging liquid-liquid extraction processes is that the available correlations of Sh are mostly empirical and valid only within the range of hydrodynamic/thermodynamic conditions in which they were determined, sometimes far from the one considered. Moreover, except for very specific conditions, correlations are mainly available for Sh I and Sh E , that rigorously cannot be used in the double-film model to calculate the Sherwood number of conjugated problems.

In this study, we aim at deriving general correlations for the internal and external mass transfer coefficients Sh i and Sh e valid for conjugated problems, and hence convenient to calculate Sh from Eq. (4.3). As in our previous work Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], the study relies on extensive use of DNS of the mass transfer from a spherical droplet in uniform and steady flow (Sec. 4.2). The parametric study has been extended, and adapted, to explore a wider range of hydrodynamic and thermodynamic conditions, and to better represent liquidliquid systems encountered in solvent extraction processes. Typical solute distributions encountered in Oil-in-Water (O/W) and Water-in-Oil (W/O) configurations, for both cases where the thermodynamic equilibrium is favorable or opposed to solute transfer, are analyzed regarding state-of-art correlations (Sec. 4.3), allowing a better understanding of the coupling between local conditions and mass transfer resistances. Finally, based on physical principles applying to both asymptotic and coupled cases, new correlations are proposed for Sh i and Sh e (Sec. 4.4), which when used in Eq. ( 4.3) leads to an excellent agreement between the estimate of the apparent Sherwood number and DNS results.

Simulation strategy

Governing equations and mesh

The study is based on the direct numerical simulation of hydrodynamics and mass transfer across a rigid spherical droplet embedded in a uniform flow. The problem is limited to 2D axisymmetric conditions.

The dimensionless Navier-Stokes equations for an incompressible Newtonian fluid and the solute mass balance equation are solved both in the dispersed, δ = i, and the continuous, δ = e, phases.

∇ • U δ = 0, (4.4) ∂U δ ∂t + U δ • ∇U δ = -∇P δ + 1 Re δ ∇ 2 U δ , ( 4.5
)

∂C δ ∂t + ∇ • (U δ C δ ) = 1 P e δ ∇ 2 C δ . (4.6)
These equations have been made dimensionless considering the droplet diameter d = 2R, the external flow velocity in the reference frame of the droplet U ∞ (the slip velocity), while the concentration is normalized by the difference between the initial concentration inside the droplet and the initial concentration of the external fluid. Specific boundary conditions are set at the liquid-liquid interface to ensure the flux continuity for momentum and mass. They express in particular the continuity of the tangential speed (U θ ) and the shear stress (µe rθ ) at the fluid-fluid interface, while the normal component of the speed (U r ) is set to zero. The normal stress balance at the interface is not necessary since there is no deformation of the spherical droplet. Consequently, the hydrodynamic jump conditions at the interface are expressed as follows:

U r,e = U r,i = 0, U θ,e = U θ,i , µ e e rθ,e = µ i e rθ,i , (4.7)

Additionally, the mass flux local density of the transferred species is assumed to be continuous, and thermodynamic equilibrium prevails at the interface. The following 120CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER equations are used to calculate the interface concentration and concentration gradients normal to the interface.

C s i = kC s e , -D i ∇C i = -D e ∇C e , (4.8)
The equations are solved with the IMFT code JADIM Legendre [START_REF] Legendre | Quelques aspects des forces hydrodynamiques et des transferts de chaleur sur une bulle sphérique -Partie 3[END_REF], which implementation for viscous particles has recently been implemented and validated in Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF].

The spatial discretization is based on a second-order centered scheme. A three-step Runge-Kutta explicit scheme (implementing a Crank-Nicolson scheme to calculate the diffusive terms) is used to compute the time-advancement, with the time-step chosen to satisfy the CFL stability criterion.

The internal and external Sherwood numbers can be deduced from the simulation results (see Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF] for the details of the calculation):

Sh i = 1 S(C i -C i s ) S ∂C ∂ξ 2 dS, ( 4.9
)

Sh e = 1 S(C e s -C ∞ ) S ∂C ∂ξ 2 dS, (4.10)
with C i the average concentration in the droplet, and C i s and C e s the average concentrations on both sides of the interface. ξ 2 is the radial coordinate.

An orthogonal curvilinear coordinates system is used, where the droplet is discretized with a polar mesh while the cell frontiers follow the equipotential and streamlines of the potential flow around a cylinder in the continuous phase. An expansion factor is used to ensure sufficient refinement at the interface, i.e. at least four grid points are located inside the hydrodynamic and/or the mass transfer boundary layers, the thickness of which varies with Re -1 2 and P e -1 2 respectively (see the mesh convergence analysis in Tab. 4.1 and in Fig. 4.2). The computational domain is centered on the droplet and extends from -50R to +50R in the axial direction, and from 0 to 50R in the radial direction, with R the droplet radius. Similar meshes and domain sizes have been used previously Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], Godé et al. [START_REF] Godé | Bassetboussinesq history force acting on a drop in an oscillatory flow[END_REF] without noticeable effect on the flow dynamics. It is important to note that the singularity at the droplet center causes some numerical scheme instability. We had to create a specific boundary condition that caused a small deviation in the flow at this point visible in the post-treatment, we believe from our hydrodynamic study validation that this has a negligible impact. To test mesh independence, we have done some simulations at the refinement level defined for our study (normal), double refinement (refined) and half refinement (coarse). The refinement level is defined by the number of points in the mass transfer or hydrodynamic boundary layer, this represents 2, 4 or 8 points for respectively coarse, normal and refined mesh. The first mesh cell close to the interface has a typical width of 10 -3 R in each case, and the total grid points evolves from 3.3 10 3 (coarse), 8.4 10 3 (normal) to 31.2 10 3 (refined). The goal is to use Richardson extrapolation Roache [START_REF] Patrick | Verification and validation in computational science and engineering[END_REF] to deduce the result (Sherwood) at infinite refinement, and therefore compute the error at each refinement level. The difference between our default mesh strategy and Richardson's extrapolation is less than 1%, so we consider our meshing strategy validated.

Case

Definition of dimensionless numbers

By convention, the problem is described considering the external Reynolds and Péclet numbers, Re and P e, defined with the external fluid properties: All ratios are defined as the property of the fluid in the droplet reported to the same property in the external fluid, e.g.:

Re = Re e = 2ρ e U ∞ R µ e , ( 4
µ * = µ i µ e . ( 4.13)
Similarly and although it is common to express the relative affinity of the solute for the solvent with Henry's coefficient He (i.e., the ratio of the solute concentration achieved at thermodynamic equilibrium between the solvent phase and the aqueous phase), and for consistency with the other fluids properties, we consider in this study the ratio between the solute concentration, at equilibrium, between the droplet and the external phase:

k = C eq i C eq e , ( 4.14)
meaning that k = He in the case of an organic droplet in water (µ * > 1), which is the most common configuration encountered in solvent extraction columns, and k = 1/He in the reversed situation (µ * < 1).

In a solvent extraction process, extraction designates the process of transferring a solute from the aqueous phase to the organic phase, for which it has a greater affinity. By opposition, stripping consists in transferring the solute back from the organic to the aqueous phase. As we can disperse either the organic phase or the aqueous phase (although this configuration is less common because it generates more pressure droplets along the extraction column), this results in four possible combinations, and therefore four different types of configuration to simulate (Fig. 4.3). At last, for a given chemical system, the molecular diffusivity of the solute cannot be varied independently of the viscosity of the carrier phase, both properties being related by the Stokes-Einstein law Eq. (4.16) Einstein [START_REF] Einstein | Über die von der molekularkinetischen theorie der wärme geforderte bewegung von in ruhenden flüssigkeiten suspendierten teilchen[END_REF]:

D = kT 6πRµ , ( 4.16)
with r the characteristic size of the solute's molecule. Hence, the ratio of the internal and external diffusivity, D * , is related to the viscosity ratio through D * = 1/µ * , which reduces the number of cases to consider, compared to our previous study Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF].

Methodology

The methodology is the same as in Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. The hydrodynamic problems are considered first, by solving Eq. (4.4) and (4.5) for the desired sets of hydrodynamic parameters Re, and µ * . The density ratio ρ * is fixed at 1. This is indeed a good approximation in liquid-liquid systems, and it was shown that this parameter has little impact on the flow Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF].

Then the mass transfer problems are considered, by solving the transient solute balance Eq. (4.6) on the frozen flow fields corresponding to the selected Re and µ * and for the chosen "chemical" parameters P e and k. We remind that by convention, Re and P e are related to the external phase. As indicated in Sec. 4.2.1, the refinement of the mesh is adjusted to the hydrodynamic boundary layer thickness, which depends on Re value. When necessary, i.e. for cases with the highest P e, further refinement of the mesh in the boundary layer is performed prior to solving the mass transfer problem. In this case, the flow-fields obtained on the "coarse" mesh are interpolated on the "refined" and additional hydrodynamic simulations are performed, to achieve steady-flow on the new mesh.

The simulations are initialized as follows. For the hydrodynamic problem, velocity is set to zero in the whole domain, while the slip velocity U ∞ is applied at the domain inlet. The pressure reference is fixed at the domain outlet. For the mass-transfer simulations, an initial normalized uniform concentration equal to C i = 1 is imposed in the droplet, while it is set to C e = 0 in the continuous phase. As depicted in Fig. 4.4, the different Sherwood numbers quickly reach their steady values, although the mean concentration in the droplet decreases. They remain unchanged until the concentration in the droplet reduces to very small values. The calculations are stopped when the mean concentration in the droplet reaches C i = 10 -5 . In each case, the corresponding values of the steady-state Sherwood numbers are such that the Sherwood number variation between two consecutive time steps satisfies ∆Sh/Sh < 10 -9 .

Parametric study

The effects of the slip velocity U ∞ , the inner, and outer phases viscosity (and their diffusivity according to Eq. (4.16)), the solute affinity and the convection are studied by varying Re, µ * , k and P e, respectively. We investigated 42 hydrodynamic configurations, and from them 1050 mass-transfer cases, by considering every possible combinations between the dimensionless parameters in the following sets:

• Re ∈ [1, 10 , 20, 50, 100, 200, 500]
• µ * ∈ [0.05, 0.1, 0.2, 5, 10, 20], and hence D * ∈ [20, 10, 5, 0.2, 0.1, 0.05]

• k ∈ [0.1, 0.25, 1 , 4, 10]
• P e ∈ [10, 50, 100, 500, 1000]

For each parameter, the range was selected to cover for solvent properties already used in hydrometallurgical processes (as e.g. the PUREX process for nuclear fuel treatment) and for possible new solvents considered in the emerging rare earth or lithium recycling sectors. The physical and chemical properties were estimated according to the PubChem database Kim S [START_REF] Cheng T Kim | Pubchem in 2021: new data content and improved web interfaces[END_REF]. 44 additional simulations have been carried out at P e = 0.01 and P e = 10, 000 to explore further the asymptotic behavior of the external Sherwood number, Sh e , in Sec. 4.4.

To verify that the non-deformable spherical droplet assumption of the model remains 126CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER valid in each case, the maximum Re was bounded based on the Clift's diagram Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF] and the Tadaki number Myint et al. [START_REF] Myint | Shapes of Single Drops Rising Through Stagnant Liquids[END_REF]. The range of Re explored is however large enough regarding the slip velocities typical of usual liquid-liquid contactors (pulsed column, stirred column, mixer-settler).

Analysis of mass-transfer regimes

In the following, the effect of Re and P e on the solute spatial distribution at steady-state is discussed for both O/W (µ * > 1) and W/O (µ * < 1) configurations. The solute has either a higher (k < 1) or a smaller (k > 1) affinity for the surrounding phase. Cases where thermodynamic equilibrium yields continuity of concentration across the interface (k = 1), as in heat transfer problems, are also considered.

The section is divided in two parts, depending on the criterion k √ D * proposed by Wylock et al. [START_REF] Wylock | Gas adsorption into a spherical liquid droplet: numerical and theoretical study[END_REF] . The steady-state contours of concentration obtained by DNS, and the concentration profiles along and perpendicular to the flow direction (respectively corresponding to the green and red lines in Fig. 4.2) are given in each case. The concentrations are normalized by the mean concentration in the droplet at the considered time

(C = C C i
). The same color-scale is used for all contour plots, ranging from 10 -5 (dark blue) to 10 (dark red).

For the 24 case-study used for illustration (corresponding to 6 different liquid-liquid systems explored for the same combinations of Re and P e values), the time evolutions of the internal, external and apparent Sherwood numbers are reported in 7.4. The Sherwood values achieved at steady-state are gathered in Tab. 7.1, where they are moreover compared with the corresponding Sh I and Sh E values predicted by the stateof-art correlations for the isolated internal and external problems, from Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] and Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF] respectively.

Internal problem and problems with expected high inter

- nal resistance, k √ D * < 1

State of the art regarding the isolated internal problem

If k √ D * << 1, mass transfer is expected to be governed by the internal resistance. In this "internal regime" the Sherwood number is noted Sh I .

The first theoretical prediction for the problem of pure diffusion inside a sphere was proposed by Newman [START_REF] Albert B Newman | The drying of porous solids: diffusion and surface emission equations[END_REF]. In this Re = 0 and P e = 0 configuration, the Sherwood number was found to tend to the asymptotic value Sh N ewman given by Eq. (4.17):

Sh I(Re=0,P e=0) = Sh N ewman = 2π 2 3 ≈ 6.58. (4.17)
The study was extended by Kronig et al. [START_REF] Kronig | On the theory of extraction from falling droplets II[END_REF] who derived an analytical solution for the two extreme cases of zero and infinite Péclet number using the solution of the

ANALYSIS OF MASS-TRANSFER REGIMES 127

Stokes equations inside a spherical droplet. At P e = 0, they obtained the same value as Newman [START_REF] Albert B Newman | The drying of porous solids: diffusion and surface emission equations[END_REF], whereas for P e → ∞ another asymptotic value was found:

Sh I(Re=0,P e→∞) = Sh Kronig = 17.9. (4.18)
More recently, Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] considered the case of a moving bubble in the range 0 ≤ Re ≤ 100, for finite Péclet values, between 0 and ∞. Based on numerical simulations, they proposed a more general correlation for Sh I :

Sh I = Sh I(Re,P e→∞) -Sh I(Re,P e=0) 1 + exp(-α 1 (ln(P e max ) -β 1 )) + Sh I(Re,P e=0) , (4.19)
with Sh I(Re,P e=0) = Sh I(Re=0,P e=0) = 6.58, with no Re correction of Eq. (4.17), and Sh I(Re,P e→∞) = 17.9(1 + Re) 0.0044 , with slight effect of the Reynolds number compared to Eq. (4.18).

Eq. ( 4. [START_REF] Collins | A second approximation for the velocity of a large gas bubble rising in an infinite liquid[END_REF]) uses an activation function to deduce Sh I from the asymptotic values at zero and infinite Péclet numbers. To improve the agreement with experimental observations, the authors additionally accounted for the interfacial advection by introducing P e max , based on the maximum tangential velocity U max at the interface as:

P e max = U max d D e = U max U ∞ P e, ( 4.20)
with for the case of bubble (validity will be discussed for droplet in Sec. 4.4.2): At moderate Re and P e (Fig. 4.5a) inner concentration contours are typical of a diffusion process, from the center of the droplet, where the concentration is maximal, toward the interface, which is consistent with the high viscosity of the inner fluid (µ * = 5). Diffusion also prevails in the external fluid (small P e), resulting in the appearance of a halo of solute around the droplet, and a quite uniform external concentration, superior to the internal one, and therefore opposed to diffusion mass-transfer. This behavior is also described by the Fig. 4.6a, where profiles in the axial (green line) and perpendicular (red line) directions are almost superimposed. The jump in concentration at the interface is consistent with k = 0.1. In this case, the Sh i and Sh values are close, assessing the impact of the internal process. However, the mass-transfer kinetics appears to be half slower than predicted by correlation related to the internal problem (Sh ≈ 1/2Sh I ).

U max U ∞ = 1
When the value of the external P e is increased (Fig. 4.5b), convective transport appears to be more effective in both phases. Sh e increases by an order of magnitude (see Tab. 7.1) and the halo surrounding the droplet is significantly reduced, in particular in front of the droplet, compared to the case at P e = 10. The effect of the internal recirculation is also reinforced, yielding an increase of Sh i and qualitatively by the displacement of the maximum concentration in the center of the Hill vortex. Oscillations are observed in the time evolution of Sh i (Fig. 7.4). This results in a drastic modification of the concentration profiles (Fig. 4.6b) with in particular the appearance of a bump on the red profile, corresponding to the cores of the Hill vortex, where the concentration is maximum and more than twice higher than the average relative concentration. The green At Re = 200, an external recirculation appears behind the viscous droplet (Fig. 4.5c and 4.5d), in agreement with the previous results of Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. The latter may act as a reservoir for the solute, as evidenced in the Fig. 4.6b. This results in a concentration distribution which decreases from the back to the front of the droplet, acting as if a reverse transfer were occurring. At fixed P e, increasing the Re has little effect on the internal, external and apparent Sherwood numbers (Tab. 7.1), and the same trends are observed at Re = 200 and Re = 20 when varying convection strength from P e = 10 to P e = 1000. As for Re = 20, at high P e, equivalent resistance to mass transfer is opposed by the two phases, and the solute spatial distribution deviates strongly from that expected for an internal regime.

Hence, besides the low value of k √ D * , strong deviations, up to 50% are evidenced at low P e between Sh and Sh I (see Tab. 7.1), due to non-negligible external resistance, and at high P e, the mass-transfer process shares similarities to coupled regime. At low P e, the contour plots (Fig. 4.7a and 4.7b) and the concentration profiles (Fig. 4.8a) evidence more significant concentration heterogeneity than in the previous case. A concentration gradient develops along the flow direction inside the droplet, leading to a noticeably steep profile (green curves). The local large concentration at the rear of the droplet, already observed at µ * = 5, is more pronounced at µ * = 0.2, due to the higher viscosity (hence lower diffusivity) of the external phase. This leads to a massive accumulation of the solute extracted in the first instants of transfer, responsible for the appearance of the axial concentration gradient highlighted in the droplet, corresponding to solute counter-transferred from the back to the front. This has a great impact on the external concentration, and therefore on the global Sherwood number, as highlighted by its time-evolution, similarly to the external mass transfer coefficient, h e (see Fig. As in the previous O/W case, increasing P e has a significant impact on the solute spatial distribution and hence on mass-transfer: appearance of a bump due to solute accumulation in the core of the Hill vortex is observed (Fig. 4.8b) as well as a drastic increase, of two orders of magnitude, of Sh e and consequently on Sh (Tab. 7.1). However, increasing P e also reduces the area experiencing the counter-transfer, therefore limiting its contribution to the Sherwood numbers. This is moreover illustrated by the comparison of the concentration profiles between P e = 10 and P e = 100 (7.3) and the subsequent tenfold value of Sh e (hence Sh).

At last, the effect of Re remains moderated. Note that in this case (µ * < 0.2), there is no external recirculation at Re = 200, in agreement with previous observations of Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF].

In this configuration typical of a conjugated problem, important discrepancies are observed, especially between the actual Sherwood and the value for an isolated internal problem, at P e = 10, due to the concentration gradient in the droplet tuning the problem to an external limited one, whereas at elevated P e the problem is clearly a coupled one. At much higher P e the problem is clearly in the coupled regime (resistance to transfer is present on both sides of the interface). It is interesting to notice that only small differences are found between Sh i and Sh I , and between Sh e and Sh E (Tab.

External problem and problems with expected high ex-

ternal resistance, k √ D * > 1

State of the art regarding the isolated external problem

The transfer process is expected to be governed by the external transfer only when k √ D * ≫ 1. The Sherwood number is noted Sh E . As for the internal asymptotic case, correlations exist to estimate Sh E from the flow and fluid properties. For external problems, Sh E grows with the Péclet number, and different relations depending on Re and P e are collected in Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF].

In the pure diffusion problem, an analytic solution can easily be derived from the governing equation, Sh E(Re=0,P e=0) = 2. From expansion series analysis, Acrivos and Taylor [START_REF] Acrivos | Heat and mass transfer from single spheres in stokes flow[END_REF] have extended this relation to small but finite Péclet and Reynolds numbers:

Sh E(Re≪1,P e≪1) = 2 + 1 2 P e + 1 6
2 + 3µ * 2 + 2µ * P e 2 ln(P e) + O(P e 2 ). (

Eq. (4.23) is valid for a droplet of any viscosity but remains based on very limiting assumptions for both Re and P e.

For larger Re, Boussinesq [START_REF] Boussinesq | Calcul du pouvoir refroidissant des courants fluides[END_REF] derived an analytical solution for bubble valid for P e ≫ 1:

134CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER Sh E(Re≫1,P e≫1) = 2 √ π P e 1 2 , (4.24)
and Clift et al. [START_REF] Clift | Bubbles, Drops, and Particles[END_REF] proposed a relation to relate low Re and high Re.

Sh E(Re,P e) = 1 + (1 + 0.564P e 2 3 ) 3 4 , (4.25)
For solid particle, Ranz and Marshall [START_REF] We Ranz | Evaporation from droplets[END_REF] proposed a correlation to connect small Re and large Re as:

Sh E(Re,P e ) = 2 + 0.60Re 0.5 Sc 0.33 , (

with Sc = ν e /D e the Schmidt number, and ν = µ/ρ the kinematic viscosity.

The most general correlation available for droplet is that of Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF], who correlated the growth rate of Sh E with P e to the viscosity ratio (i.e. the diffusivity ratio) through three expressions relating respectively to µ * = 0 (spherical bubble), µ * = 2 (fluid sphere) and µ * → ∞ (solid sphere), in Eq. (

:

Sh E = 2 -µ * 2 Sh E(µ * =0,P e,Re) + 4µ * 6 + µ * Sh E(µ * =2,P e,Re) 0 ≤ µ * ≤ 2, Sh E = 1 µ * + 2 Sh E(µ * =2,P e,Re) + µ * -2 µ * + 2 Sh E(µ * →∞,P e,Re) µ * > 2. ( 4.27)
With: The concentration discontinuity at the interface is consistent with k = 10 assessing the higher affinity of the solute for the fluid in the droplet. The inner concentration profiles are not affected by the hydrodynamics (no impact of Re). However, increasing the P e again results in a steeper decrease in the outer concentration perpendicular to the flow direction, and a flatter profile at the rear (Fig. 4.12b) compared to the case at moderate P e (Fig. 4.12a). The effect of Re is hardly noticeable for a given P e. Besides the analogy of the solute spatial distribution with the features of the asymptotic external regime, and the very close values of Sh e and Sh E , the approximation Sh = Sh E is not valid because using Sh E and Sh I in Eq. (4.3) overestimates the apparent Sherwood number. The time evolution of the Sherwood numbers in this case is shown in Fig. 7.7. The concentration contours are similar to cases in Fig. 4.5 for the same viscosity (resp. diffusivity) ratio, leading to very similar internal and external Sherwood numbers, although of course the apparent Sh is significantly lower here due to the much less favorable thermodynamic coefficient k = 10. As in the equivalent case for k = 0.1, the oscillations due to the Hill vortex can be observed in the time-evolution of Sh i at high P e (see Fig. The conservation of mass flux at the interface implies that all resistances to mass transfer are added, reading Eq. (4.3). However, in the double-film model (Fig. 4.1), the Eq. ( 4.3) at the interface is expressed by Eq. (4.29), where it is assumed that Sh i = Sh I and Sh e = Sh E , the Sherwood numbers considering the internal and external problems being calculated by one of the abovementioned correlations.

Sh E(µ * =0,P e,Re) = 0.651P e 1/2

Second case-study

1 Sh = 1 Sh I + kD * Sh E . ( 4.29)
As evidenced by Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF] Eq. (4.29) may imply important errors in actual coupled mass transfer problems. Indeed, as it was pointed in the previous DNS examples, the solute distributions inside, outside and/or at the droplet interface generally deviate from assumptions based on a fixed and uniform concentration at the interface for internal and external problems. The differences between Sh i and Sh I , between Sh e and Sh E , and between Sh Eq.( 4.3) and Sh Eq. (4.29) for the illustration case-study are given in Tab.7.1 and reported as error metrics in Tab. 4.2 for the 1050 configurations studied by DNS, in terms of:

• Max error, the maximum difference between predicted (Sh I , resp. Sh E ) and simulated (Sh i , resp. Sh e ) among the 1050 cases;

• Mean Absolute Error (MAE), the average of the differences measured between predicted and simulated values for all cases;

• Mean Absolute Percentage Error (MAPE), similar to MAE but based on the relative differences, i.e. divided by the value from DNS;

• Determination coefficient R 2 , that measures the proportion of the variance in the simulation data that can be explained by the prediction. This is moreover illustrated by the parity plots between the usual correlations and the DNS results in Fig. 4.17. One can see that the Sh I correlation, Eq. (4.22), by constraining the value between two asymptotes induces a systematic source of error, while the correlation for Sh E , Eq. (4.27), better describes the trends but generally underestimates the actual external Sherwood number Sh e .

Using these correlations to deduce the apparent Sherwood number Sh by Eq. (4.29) leads to significant errors, often exceeding 50% of absolute percentage error, hence confirming the observations of Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF] for configurations typical of liquid-liquid extraction, and the need for specific correlations of Sh i and Sh e to be proposed in coupled problems.

Correlation of the internal mass transfer coefficient for coupled problems, Sh i

Maximum interfacial velocity

As mentioned in Sec. 4.3.1, Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] proposed to use the maximum velocity U max at the interface to correct for the effects of interfacial flow on mass-transfer. Indeed, compared to solid particles, for which no-slip condition at the interface imposes a zero-tangential velocity, it is well-known that non-zero interfacial velocity prevails in the case of "fluid particles". In the case of droplets, Hadamard [START_REF] Hadamard | Mouvement permanent lent d'une sphère liquide et visqueuse dans un liquide visqueux[END_REF] determined that the maximum tangential velocity U max tends toward 1 2 U ∞ when the Reynolds number tends to zero, and 3 2 U ∞ when the Reynolds number tends to infinity. The link between these two asymptotic values was expressed using Eq. (4.21) by Legendre [START_REF] Legendre | On the relation between the drag and the vorticity produced on a clean bubble[END_REF] in the case of a bubble.

In the present case of droplet, we seek to relate the asymptotic values of U max to the viscosity ratio µ * , which tends to 0 for a bubble and to infinity for a solid sphere. By plotting the data from our simulations in a semi-log graph, Fig. 4.18, a sigmoid-like trend is observed. Using the same approach as proposed by Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF], we fitted Eq. (4.30) to our data by least squares minimization and obtained:

U max,µ * = 1 2
16 + 3.315Re We suppose here that similar asymptotic values exist for the internal Sherwood number Sh i , of which Sh I is expected to be the limiting case at k √ D * ≪ 1.

We have also discussed in Sec. 4.3, that the assumption used to define this "internal mass-transfer" problem and to derive the corresponding asymptotic values Sh I , is incorrect in most coupled transfer cases:

• Concentration is generally not uniform along the interface and around the droplet in the continuous phase. The deviation from this assumption is particularly highlighted in Water-in-Oil configurations (as in Fig. 4.15), but it can also be observed for more viscous droplets (Fig. 4.9) even for expected internally limited regimes (see Fig. 4.5);

• Although little effect of the Reynolds number is observed in most cases, the simulations clearly show an effect of the flow on transfer, with a significant impact of the Hill vortex (at high P e) on the solute distribution in the droplet and hence at the interface (Fig. 4.6);

• An impact of the solute affinity k was revealed in particular at low P e, where significant differences are observed in the distribution of the solute, and therefore on the value of Sh i (Tab. 7.1), depending on whether the solute has a higher (Fig. 4.7) or a lower (Fig. 4.11) affinity for the external phase, with an enhanced (k = 0.1) or a damped (k = 10) effect on the "counter-transfer" likely to occur at the rear of the droplet.

Neglecting the impact of the actual solute distribution and of the equilibrium constant k on the internal mass-transfer process is responsible for significant discrepancies between the asymptotic values of Sh i at P e → 0 and P e → ∞ and the corresponding Sh I values predicted respectively by Newman [START_REF] Albert B Newman | The drying of porous solids: diffusion and surface emission equations[END_REF] and Kronig et al. [START_REF] Kronig | On the theory of extraction from falling droplets II[END_REF] for the internal problem, and given by Eq. (4.17 As the effect rapidly decreases with P e, we therefore correlated it to the value of P e max defined in Eq. (4.20) and based on the maximum velocity at the droplet surface, U max,µ * , given by Eq. (4.30), by analogy with the expression proposed by Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] for the internal problem.

Correlation of the external mass transfer coefficient for coupled problems, Sh e

As mentioned in Sec. 4.3.2 there is no finite asymptotic value for the Sherwood number of the external problem, Sh E . It is a growing function of the Péclet number proposed by Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF]. Similarly, all the DNS results presented in Sec. 4.3 confirm the strong impact of an increase of P e = 10 to P e = 1000 on the value of Sh e .

Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF] distinguished several cases in their correlation for Sh E , Eq. (4.27), depending on the viscosity ratio, hence to express the case of droplets (represented by µ * = 2) to the simplified cases of bubbles (µ * → 0) and solid particles (µ * → ∞). We compare, in Fig. 4.20, the value of Sh e obtained from DNS to the predictions of Eq. (4.27) for Sh E . The equilibrium constant, k, has also a non-negligible influence on Sh e , especially at low viscosity ratio (Fig. 4.20a), where the internal resistance is supposed to be low. It is interesting to observe that in this case, the higher the value of k, the closer Sh e from Sh E , consistently with the unfavorable thermodynamics, while it seems to be the opposite trend in the case of droplets (Fig. 4.20b), at least at moderate P e. At last, the effect of k appears to be negligible for solid particles (Fig. 4.20c). Further analyzing the evolution of Sh e with k in the case of bubbles (see Fig. Accordingly, we defined new numerical parameters for Eq. (4.28), by proceeding to a least squares optimization on the entire set of DNS results. From there, we deduced an expression of the external Sherwood number, Sh e , based on the correlation of Michaelides [START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF] for the external problem:

Sh e = 2 -

Relevance of the new correlations

The proposed correlations were tested on the whole set of 1050 simulations, resulting in a significant improvement of the prediction of both the internal, Sh i , the external, Sh e , and the apparent, Sh, Sherwood numbers. This is illustrated by the parity plots in Fig. 4.22, where most predictions fall under 10% error margin, which is a great improvement compared to the previous correlation results (Fig. 4.17).

Conclusion

Mass transfer rate from a droplet has been mostly studied analytically under very restrictive assumptions, and empirically, either in dedicated devices, leading to hydrodynamic conditions usually different from real processes, or from global measurements in pilot devices. for the development of liquid-liquid extraction processes. Today, for recycling purposes, we need to develop some new processes with the state-of-the-art solvents that are, for the most cases, more viscous and therefore not in the range of parameters considered for the development of correlations currently in use. In this study, we have considered the mass transfer from a single isolated spherical droplet thanks to DNS simulations, to cover a wide range of parameters, such as the Reynolds and Péclet numbers, the Henry coefficient and the viscosity ratio, for the development needs of future processes. This study holds significance for heat transfer, as it has also been simulated cases with a Henry coefficient to 1. In the case of liquid-liquid extraction, we need to master the expression of the transfer resistance inside and outside the droplet. Until now, most of the studies have used correlations developed for internal or external transfer regimes only: one phase transfer is quicker than the other, so we can consider that one phase resistance is zero. Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF], and our simulations have demonstrated that using those correlations in the transfer resistance additivity law can generate significant discrepancies. So, our objective was to understand the origin of these discrepancies and to provide more accurate expressions for the description of the mass transfer across the droplet surface.

Thanks to DNS simulations, we have observed that the key mechanism to consider is the presence of a solute reservoir at the rear of the droplet due to the combined effect of the flow both inside and outside the droplet. As a consequence, the lowest resistance phase does not always have a constant and uniform concentration at the interface, as stated when combining the internal and external problems.

Because the internal Sherwood correlation could be described as an activation function between two asymptotic values, we needed to revise the asymptotic values, and for 152CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER the activation dynamics to add the impact of the reservoir, creating, in some cases, a counter-transfer in the correlation. Using those observations in the form of a new formulation of the correlation in an optimization code, we were able to fit our simulation data with a R 2 score: R 2 = 0.971 (instead of R 2 = 0.623 previously with the original correlations).

For the external Sherwood correlations, the state of the art was based on some least squares optimizations and did not have as much physical meaning as for the internal Sherwood correlation, so we used the same method as previous authors to fit our simulation data, and we were able to achieve a R 2 score: R 2 = 0.989 (instead of R 2 = 0.859 previously with the original correlation). Using those new correlations in the mass transfer resistance additivity law, we obtained a R 2 score: R 2 = 0.992 for the global Sherwood (instead of R 2 = -0.680 previously with the original correlations).

We have limited our study to the restrictive case of axisymmetric and uniform flow around an isolated spherical droplet. Some studies have indicated that the flow inside a spherical droplet in a uniform flow could be 3D for some specific Reynolds numbers and viscosity ratios, see Edelmann et al. [START_REF] Edelmann | Numerical investigation of different modes of internal circulation in spherical drops: Fluid dynamics and mass/heat transfer[END_REF], Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF]. In those cases, the developed mass transfer correlations are invalidated and there is a need to consider the 3D flow influence on mass transfer. It would also be important to study the shape oscillation effect as found in Charin et al. [START_REF] Amanda | On the dynamic behavior of rising droplets[END_REF] on the mass transfer. We also need to understand the collective effect on mass transfer, for which a promising approach seems to be the use of neural networks, see Chouai et al. [START_REF] Chouai | Use of neural networks for liquid-liquid extraction column modelling: An experimental study[END_REF].

SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

In most industrial applications, mass transfer is accompanied by a chemical reaction, likely to strongly modify the apparent transfer reaction. In this chapter, we compare different methods to predict the mass transfer enhancement induced in the configuration of Chap. 4, when a first order chemical reaction occurs in the external phase.

The first method is the classical film theory methods, which is an analytical solution to the problem with the assumption of flat interface and diffusive transfer only. The second one relies on the same assumptions, but considers the effect of sphericity. The third one is also based on the film theory method, but considers the convective part in the solute transfer rate, thanks to a more general definition of the Hatta number. In all cases, the application of the film theory is based on the correlations introduced in Chap. 4 for the estimation of the internal and external resistances without reaction/ It is found that the third method precisely predicts the enhancement factor due to the chemical reaction on mass transfer. This chapter will be submitted shortly as a research article.

Introduction

Mass transfer in reactive fluid-fluid systems is of most importance for many applications: liquid-liquid extraction as in Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] or Brunson and Wellek [START_REF] Brunson | Mass transfer inside liquid droplets and gas bubbles accompanied by a second-order chemical reaction[END_REF], gas-liquid extraction Lebrun et al. [START_REF] Gaelle Lebrun | Gas-liquid mass transfer around a rising bubble: Combined effect of rheology and surfactant[END_REF] or Sa'adiyah et al. [START_REF] Setiorini Sa'adiyah | Effects of chemical absorption on mass transfer from single carbon dioxide bubbles in aqueous sodium hydroxide solution in a vertical pipe[END_REF] for CO 2 absorption. Because it is a complex problem, most studies are empirical, aimed e.g. at measuring the overall apparent coefficient k L a in bubble columns Sa'adiyah et al. [START_REF] Setiorini Sa'adiyah | Effects of chemical absorption on mass transfer from single carbon dioxide bubbles in aqueous sodium hydroxide solution in a vertical pipe[END_REF], or consider simplified flow configurations that enable analytical solutions (film theory) or simple numerical resolution, i.e. without coupling with the full Navier Stokes equations. Such configurations generally deviate from typical conditions prevailing in industrial applications and leading to coupled mass transfer situations studied numerically in Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] and Kleinman and Reed [START_REF] Kleinman | Unsteady conjugate mass transfer between a single droplet and an ambient flow with external chemical reaction[END_REF].

In the framework of the double film theory Schaschke [101], as shown in Fig. 5.1, the system is seen as a one-dimensional problem corresponding to planar interface, where steady transfer of the solute is represented by linear concentration gradients through the thin layers (films) located on either side of the interface separating the two immiscible fluids. Thermodynamic equilibrium is assumed at the interface. For convenience, the mass-transfer rate Φ (in kg.s -1 ) is expressed from the difference between the bulk concentrations (i.e. at "infinity", in kg.m -3 ), introducing a global mass transfer coefficient or conductance, h (in m.s -1 ), expressed by Eq. (5.1):

Φ = h S (C ∞ i -kC ∞ e ), (5.1)
with S the total contact area (in m 2 ).

In this model, each film is considered as a stagnant fluid layer of width δ where the solute transport is achieved by molecular diffusion only, hence leading for the left side:

h i = D i δ i . (5.2)
with D i the molecular diffusivity of the transferred solute in the fluid, and similarly for the other side. The apparent resistance 1/h results from the addition of the two resistances, 1/h i and 1/h e , placed is series.

Regardless of the transfer mechanisms on each side of the interface (convection and/or diffusion), the continuity of the mass transfer rate enables to consider the additivity of the resistances. In this framework, the dimensionless form of the conductance h, the Sherwood number, Sh = hd D i is considered, with d a characteristic length for the transfer process. The additivity of series resistances reads: .3) where D * = D i /D e is the ratio of diffusion coefficients of the solute in each fluid, and k is the equilibrium constant related to respective solubility of the solute in the liquids.

1 Sh = 1 Sh i + kD * Sh e . ( 5
When the solute is consumed by a chemical reaction in the receiving phase, the concentration profiles may be modified, as depicted in Fig. 5.1b.

156CHAPTER 5.

SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

• If the reaction is slow, it takes place in the bulk without impacting the mass transfer rate. The concentration's profiles are unchanged compared to Fig. 5.1a.

• If the reaction is fast, a stronger gradient prevails in the film, where the reaction takes place in a thin zone close to the interface, materialized by the dashed line in Fig. 5.1b, and this yields acceleration of the solute transfer.

The ratio between the diffusive time across the film and the characteristic time of the chemical reaction, called the Damköhler number, Da, is the relevant dimensionless number for reactive mass-transfer systems (with Da = 0 for non reactive problems).

The acceleration of transfer due to the chemical reaction can be accounted for using an enhancement factor E, to correct the conductance of the phase where the reaction occurs, phase "e" in the case of Fig. 5.1:

1 Sh Da = 1 Sh i,Da=0 + kD * E Sh e,Da=0
(5.4)

with Sh Da the apparent Sherwood number in the reactive system, Sh i,Da=0 and Sh e,Da=0 the Sherwood numbers of phase "i" and "e" without reaction. The enhancement factor E is defined as:

E = Φ Da Φ Da=0 = Sh e,Da Sh e,Da=0 (5.5)
For the studies of the opposite scenario with a reaction inside the droplet, it is possible to apply the enhancement factor to the internal Sherwood number to express the global (or apparent) Sherwood number as:

Sh(Re, P e, D * , k, Da) = 1 E Sh i,Da=0
+ kD * Sh e,Da=0 -1 .

(5.6)

E is generally correlated to the Hatta modulus, Ha, which compares the maximum flux of solute likely to be consumed by the chemical reaction within the film, Φ R,max , to the maximum flux likely to be transported through the film, Φ D,max . Its expression depends on the type of chemical reaction considered.

Ha = Φ R,max Φ D,max (5.7)
In the case of a first order chemical reaction, the film theory reads in Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] and Villermaux [START_REF] Villermaux | Génie de la réaction chimique[END_REF]: The most comprehensive study of the effect of a chemical reaction on the rate of coupled mass transfer (i.e. in cases where convection is also involved in the solute transport) was carried out by Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] in the case of a spherical fluid particle. The author considered either first order or second order reaction, taking place inside or outside the sphere.

E = E F T = Ha tanh(Ha) (5.
The study was limited to creeping flow, to simplify the numerical resolution of the governing equations, and a single intermediate value of P e = 100, and 3 values of Da, Da = 1, 100 and 1000. Thanks to this model, the relevant Sherwood numbers for both the non-reactive, Da = 0, and the reactive, Da ̸ = 0, cases were calculated and compared with the apparent Sherwood number resulting from the numerical simulations. Acceleration of transfer was deduced from Eq. (5.4), considering the enhancement factor E F T predicted by the film theory and a general definition of Ha based on the Sherwood number, Eq. (5.10) in the case of an external reaction:

Ha = Ha J = 2 √ Da Sh e,Da=0
(5.10)

Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] shows the drawbacks of this approach, by demonstrating that the additivity relations in Eq. ( 5.3) do not offer accurate predictions for all parameters values. Recently, direct numerical simulations (DNS) were used to get a better understanding of the specific features of coupled mass transfer, without reaction, and to highlight the main differences with the film theory Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. Based on these findings, accurate correlations for Sh i and Sh e , and therefrom for Sh, were proposed and validated in the case of the non-reactive coupled mass-transfer problem (Da = 0) from a spherical droplet in a uniform flow, for conditions up to Re = 500 and P e = 1000 Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], Godé et al. [START_REF] Godé | Corrigendum to Towards more predictive correlations for conjugated mass transfer across a droplet interface[END_REF].

The case study is described in Sec. 5.2 where the governing equations and numerical techniques, allowing to solve the fully coupled problem, are rapidly introduced. The simplified solution for pure diffusion transfer in an infinitely long domain (e.g. same assumption as in the film theory) is given, to derive a definition of the enhancement factor that accounts for finite spherical geometry (E SF T ).

Based on the DNS results, the effect of a first order chemical reaction on the mass transfer process is analyzed in Sec. 5.3, for increasing values of Da and various typical coupled regimes.

The performance of the modified additivity of transfer resistances model, Eq. (5.4), is discussed in Sec. 5.4, based on the non-reactive correlations of Sh i,Da=0 and Sh e,Da=0 introduced in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], Godé et al. [START_REF] Godé | Corrigendum to Towards more predictive correlations for conjugated mass transfer across a droplet interface[END_REF]. The different possible combinations of the enhancement factor and Hatta criterion given by the film theory (E F T and Ha F T Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF]), as well as its extension to spherical geometry (E SF T and Ha SF T ), are considered.

Problem statement and governing equations

The configuration under investigation is the same as in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]. We consider a single non-deformable sphere of fluid in a uniform flow of an immiscible liquid phase.

The density and viscosity of the droplet and continuous phase are respectively ρ i and ρ e , and µ i and µ e . The far-field velocity of the external fluid is U ∞ . The solute A, with initial concentration C ini i , is transferred from the droplet to the continuous phase where its initial concentration is C ini In this study, we additionally assume that the solute is consumed in the external phase by a first order chemical reaction: A → B, which rate R is given by: R = -νC e (5.12)

By convention, the thermodynamic equilibrium at the interface is expressed at any time, and any location on the surface, by:

C s i = kC s e (5.13)
with k the equilibrium constant.

∇ • U δ = 0, (5.14)
ρ δ ∂U δ ∂t + ∇ • (U δ U δ ) = ∇ • T δ + ρ δ g, ( 5.15)
with T δ the stress tensor and g the gravity. The solute mass balance in each phase is given below for the two fluids:

∂C i ∂t + ∇ • (U i C i ) = D i ∇ 2 C i , ( 5.16
)

∂C e ∂t + ∇ • (U e C e ) = D e ∇ 2 C e -νC e
(5.17)

with -νC e the rate of disappearance of A due to the chemical reaction.

At last, the droplet interface is mobile, leading to the following interfacial conditions, expressing the continuity of the tangential velocity and stress:

U r,e = U r,i = 0, U θ,e = U θ,i , µ e e rθ,e = µ i e rθ,i , , (

Continuity of the local flux density is imposed across the interface, i.e. no accumulation of the solute:

-D i ∇C i .n = -D e ∇C e .n (5.19)
where n is the local vector normal to the interface. The dimensionless formulations of Eqs. (5.14)-(5.17) are:

∇ • U i = 0, (5.20)
∇ • U e = 0, (5.21
)

∂U i ∂t + ∇ • (U i U i ) = -∇P i + 1 Re i ∇ 2 U i ( 5
∂C i ∂t + ∇ • (U i C i ) = 1 P e i ∇ 2 C i ( 5

Numerical model

The set of dimensionless Eqs. )), are solved simultaneously in the two phases using the IMFT in-house code JADIM Legendre [START_REF] Legendre | On the relation between the drag and the vorticity produced on a clean bubble[END_REF]. We use the same methodology, i.e. same computational domain, same mesh and same numerical strategy, as in the non-reactive case, which is fully described in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]. The 2D axisymmetric computational domain, the central part of which is shown in Fig. 5.3a, extends from -50R to 50R in the flow direction and from 0 to 50R in the perpendicular direction. The droplet is discretized by a polar mesh, whereas the mesh used for the external fluid is built on the equipotential and streamlines of the potential flow around a cylinder (orthogonal mesh lines). This configuration was found better than a fully polar mesh in terms of stability and computation cost Rachih et al. [START_REF] Rachih | Numerical study of conjugate mass transfer from a spherical droplet at moderate Reynolds number[END_REF]. The internal and external Sherwood numbers, defined respectively by Eq. (5.26) and (5.27), are directly obtained from the simulation results, and the global Sherwood number, Sh, is determined using the mass transfer continuity Eq. (5.3). Note that these definitions are valid regardless of the value of the Damköhler number, i.e. for both non reacting (as in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]) and reacting mass transfer cases.

Sh i = 1 S C i -C s i drop ∂C ∂ξ 2 I - dS, ( 5.26
)

Sh e = 1 S C s e -C ∞ drop ∂C ∂ξ 2 I + dS, ( 5.27)
with C i the average concentration in the droplet, C s i the average concentration on the inner side of the interface (I -), C s e the average concentration on the outer side of the interface (I + ), and ξ 2 the radial coordinate.

Mesh refinement validation

In Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF] the mesh refinement at the interface was adjusted to ensure that at least 4 grid points are located within the external mass transfer boundary layer, which width roughly scales with P e e -1 . To verify that this criterion is still consistent in the case of reactive mass transfer (where, as depicted in Fig. 5.1b, the solute concentration is likely to go to zero at a reaction plane located closer to the interface), we compared the DNS results obtained in a diffusion-dominated problem, with the corresponding analytical solution (7.5). As illustrated in Fig. 5.3b, excellent agreement is observed between the results of the numerical simulations and the analytical solution, hence validating the chosen mesh for the reactive case up to Da = 100.

Simulation strategy and explored domain

The range of parameters is the same as in our previous study dedicated to non-reactive transfer (Da = 0) Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], and the steady state of hydrodynamic simulations were used, by imposing the converged steady velocity field, to conduct additional mass transfer simulations for different values of Da ̸ = 0. This allows to obtain the steady-state Sherwood numbers, as shown in Fig. 5.4. • P e e ∈ [10, 50, 100, 500, 1000]

• Da ∈ [0, 1, 10, 100]

The pure diffusion case

Assuming by analogy with the film theory that the transport of the solute is achieved only by diffusion (U ∞ = 0), the problem reduces to a 1D spherical geometry, fully described by the following radial mass-balance equations (in spherical coordinates):

D i 1 r 2 ∂ ∂r r 2 ∂C ∂r = 0, for 0 < r < R (5.28) D e 1 r 2 ∂ ∂r r
2 ∂C ∂r -νC = 0, for r > R (5.29)

RESULTS AND DISCUSSION

163

The analytical expression (see 7.5) of the solute distribution in the spherical domain is obtained for both the reactive (ν ̸ = 0) and the non-reactive cases (ν = 0). From the concentration profiles, the fluxes of mass transfer Φ Da̸ =0 and Φ Da=0 can be derived. It is shown that the expression of the enhancement factor, Eq. (5.8), writes simply:

E SF T = 1 + Ha (5.30)
with Ha = Ha F T defined by Eq. (5.9). We notice that the internal Sherwood number is only slightly impacted by the chemical reaction (Fig. 5.5a), compared to the external Sherwood number (Fig. 5.5b), as the reaction takes place in the external/continuous phase. The apparent effect of the chemical 164CHAPTER 5.

Results and discussion

SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

reaction hence mainly comes from the increase of the Sh e , which is consistent with the proposed definition of Sh Da for coupled problems with external reaction, Eq. (5.3). Accordingly, a faster decrease of the mean droplet concentration is observed over time (Fig. 5.6). The concentration fields corresponding to the four cases of Fig. 5.5 are shown in Fig. 5.7, highlighting that, although the hydrodynamic conditions are unchanged, the extent of the external concentration gradient is significantly reduced by the chemical reaction, while both the internal and external concentration distributions tend towards spherical symmetry. Consequently, the case that could have been considered as an isolated external problem at Da = 0, progressively evolves toward an almost internal problem at Da = 100 (the internal resistance being mainly related to the low diffusivity, D * = 0.2), due to the sole effect of the solute consumption by the chemical reaction. The increasing gap with increasing Da between the mean internal concentration and the average surface concentration reported in Fig. 5.6 is also consistent with this finding. Complementing Fig. 5.7, the reaction-induced modifications of the concentration profiles, along and perpendicular to the flow direction (materialized by respectively the green and red lines in Fig. 5.3b), are illustrated in Fig. 5.8. We observe that the chemical reaction steepens the outer concentration gradient, which explains the significant increase of Sh e (defined by Eq. (5.27)). In the fastest reactive case Da = 100, an effect is also noticeable for the inner concentration gradient, as the interface concentration becomes smaller, in agreement with the slight increase in Sh i (defined by Eq. (5.26)). Despite the presence of closed streamlines inside the droplet, that are typical of a Hill's vortex shape, it is difficult to assess the coupled effect of reaction and convection on the concentration field in this case, due to the small Peclet number considered, P e e = 10.

The Peclet number has a major effect, more significant than the Reynolds number, on the (non-reactive) mass transfer rate from an isolated droplet Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]. On the other hand, no dependence of the enhancement factor on the external Peclet number can be drawn from the film theory model, since molecular diffusion is the only transport mechanism considered in this simplified model (P e = 0). To quantify the effect of both the external flow and the chemical reaction, the simulations discussed in Fig. 5.7 and 5.8 were reproduced at P e e = 1000 (Fig. 5.9 and 5.10). The time evolutions of the mean droplet concentrations for 0 ≤ Da ≤ 100 are compared in Fig. 7.10 Both the preferential accumulation of the solute at the center of the internal vortex close to the surface, rather than at the center of the droplet, and the stretching of the concentration wake observed behind the droplet (Fig. 5.9) are typical of the high P e regime. As in the case of small P e, the chemical reaction mainly modifies the external concentration profile. However, at high P e, this effect is not uniform around the droplet and is mostly observed in the wake of the droplet, where the transferred solute is likely to accumulate. (Fig. 5.10a). Note that this difference in the axial and perpendicular directions cannot be accounted for by the film theory.

Enhancement effect of the chemical reaction with expected high internal resistance according to the film theory (

√ kD * < 1)
This case with µ * < 1 is typical of water droplet in oil. In this other coupled problem, the thermodynamic equilibrium is favorable to mass transfer (k = 0.1) and the solute diffusivity in the external phase is small compared to the internal one (D * = 5), meaning that at small P e, i.e. poorly convective flow, the solute transport in the outer phase is weak. In the non-reactive case, these conditions result in a massive accumulation of the transferred solute at the rear of the droplet, compared to the previous case discussed in Sec. 5.3.1. This accumulation is also materialized by an increasing concentration in the profile along the axial direction (Fig. 5.12a). This features of the non-reactive case (Da = 0, Fig.

PREDICTING THE GLOBAL SHERWOOD NUMBER FOR CONJUGATE TRANSFER ENHAN

As in the case of the viscous droplet (Sec. 5.3.1), at such high Péclet number, a thin concentration boundary layer develops on the droplet front. The impact of the internal Hill's vortex is also well pronounced with a larger concentration at its center. In this case, however, the chemical reaction appears to have little impact on the solute transport, when comparing the concentration fields at Da = 0 and Da = 100 both inside and outside the droplet, suggesting a little enhancement effect of the chemical reaction. This is also highlighted in the concentration profiles reported in Fig. 7.11 and by the similar slopes of the time evolution of droplet average concentration shown in Fig. 7.12. The concentration in the droplet wake at Da = 100 is however reduced (see Fig. 7.11) and we can expect larger Da to efficiently consume the concentration that acumulates in the wake.

Summary

In coupled mass transfer configurations without reaction, the distribution of the solute is impacted by advection transport, with a maximum of concentration in the core of the internal Hill's vortex and in the wake of the droplet. As expected, the impact of advection is clearly increased when increasing the P e number. In the moderate Da range considered, the DNS results show that the effect of the chemical reaction (taking place in the external phase) is mainly visible in regions where the solute concentration is low, and therefore rapidely consumed. In more concentrated areas, such as the halo surrounding the droplet at moderate P e and its wake at high P e, a higher Da would probably have resulted in a more perceptible effect. The effect is more pronounced when the chemical reaction is faster. In cases where Da is high and P e low, the effect of the chemical reaction is also noticeable on the spatial distribution of the solute inside the droplet.

The chemical reaction tends to smooth out the concentration heterogeneity induced by convection in the coupled systems, which as a result approach the configuration described by the film theory i.e. where the concentration gradients are concentrated in the vicinity of the interface and almost uniform.

Predicting the global Sherwood number for conjugate transfer enhanced by chemical reaction

Predicting Sh for a fluid-fluid system is not straightforward and requires some careful considerations, depending on both the flow and chemical configurations. This section is intended to clarify the reasoning.

As for non-reactive systems Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF], DNS allows to highlight the deviation of coupled problems from the simplified asymptotic assumption of the film theory, and to relate them to the typical solute distributions achieved by the flow transport inside and outside the droplet. Although a wide range of parameters was 172CHAPTER 5.

SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

studied, only systems where the droplet remains spherical and the flow axisymmetric were investigated. It was demonstrated that using the external and internal Sherwood number correlations from isolated transfer problems (i.e. where either the internal or the external resistance is negligible, respectively Sh I and Sh E ) in the resistance additivity law, Eq. ( 5.3), could lead to important discrepancies in conjugate problems.

To go beyond, the authors proposed in a previous study Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF] physically based correlations of the internal, Sh i , and external, Sh e , Sherwood numbers prevailing in coupled problems, where the solute distribution can be far from uniform inside and outside the droplet. They also validated the accuracy of the global Sherwood number prediction, Sh, by Eq. ( 5.3).

In this last section, we propose to complement the description of coupled transfer without reaction following our approach in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF] together with the enhancement factor used in film theory, in order to develop a better prediction of the general behavior of coupled and reactive systems.

Non-reactive mass transfer coefficients

The correlation for Sh i at Da = 0 is based on the work of Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] for bubbles (i.e. µ * → 0). It is extended here to a wider range of µ * by introducing an activation function that corrects the maximal interfacial velocity (expressed by P e max ), which varies strongly between the two limiting cases of the bubble and the solid sphere (i.e. µ * → 0 and µ * → ∞, respectively):

Sh i = Sh i,P

Enhancement factor

Several definitions of Ha and E have been introduced in the previous sections, and can possibly be combined in Eq. (5.4).

• Regarding the Hatta number, Eq. (5.9), defining Ha F T , is an expression deduced from the film theory, for which the solute transport in the two phases is only driven by molecular diffusion. Alternately, Eq. (5.10), defining Ha J , is a more general expression, which accounts for convection transport. It was used in Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] for creeping flows at moderate P e, achieving good results.

• Regarding the enhancement factor, Eq. (5.8), defining E F T , is the original expression derived from the film theory Schaschke [101], i.e. based on planar geometry; whereas Eq. (5.30), defining E SF T , results from the transposition of E F T to a spherical geometry, more appropriate for a droplet, still in the framework of the "diffusion" film theory.

The evolution of the enhancement factor with the Hatta criteria for the different possible definition combinations, is illustrated in Fig. 5.14 and 5.15, where the enhancement factors obtained from the 1260 DNS results (using Eq. (5.5)), have been reported, respectively distinguished by their Da or their external Péclet P e e values, as both parameters are likely to affect mass transfer. As expected, the enhancement factors for the coupled mass-transfer configurations obtained by DNS agree much better with the prediction of the film theory, E F T and E SF T when the general definition of the Hatta number Ha J is used (Fig. 5.14b) and (Fig. 5.15b). It is however difficult to assess which of the two definitions is the most appropriate.

PREDICTING THE GLOBAL SHERWOOD NUMBER FOR CONJUGATE TRANSFER ENHAN

The deviation from the film theory with Ha F T increases significantly with Da and is up to one order of magnitude at Da = 100 (Fig. 5.14a). No trend emerges regarding the effect of the Péclet value on the discrepancy between DNS results and the film theory (Fig. 5.14a). Conversely, when the general definition of the Hatta criterion Ha J is considered, a better correlation is observed between the DNS results and the film theory models. The agreement is better at larger values of Da (Fig. 5.14b), i.e. when the chemical reaction is fast.

While no trend could be drawn about the Péclet numbers using Ha F T , Fig. 5.14b indicates that at moderate chemical reaction rate (Da < 100), the deviation between the DNS results and the enhancement factor predicted by the film theory is higher at moderate Peclet number. More precisely, the largest errors correspond to cases where 1 ≤ P e e ≤ 10 and 0.2 < k √ D * < 0.5. These configurations, an example of which is illustrated in Fig. 5.11, have been observed to undergo significant changes in the mass transfer regime when Da is increased. They evolve from an initial coupled transfer problem, where the concentration gradient is not uniform along the interface,due to the effect of the prevailing interfacial velocity in liquid-liquid systems as in Colombet et al. [START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] and Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], to a mostly internal problem where the surface concentration is uniform as in the film theory. It is likely that using information from a coupled problem at Da = 0 in the Hatta criterion Ha J to describe a mostly isolated internal case isolated at Da > 10 causes this error.

Accuracy of the model prediction

Based on the previous results, we chose to use Eq. (5.4) with the same definition as Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] for the enhancement factor, i.e. E F T = f (Ha J ). The global Sherwood number predicted in this way is compared to the value obtained from DNS in Fig. 5 The straight lines reprensent the limits of 0% error (black), 10% error (blue) and 30% error (red).

The possible effect of the external chemical reaction on the internal Sherwood number is noticeable on Fig. 5. [START_REF] Joseph H Gibbons | Effect of a surface active agent on the velocity of rise of benzene drops in water[END_REF], where a poor agreement between the model and the DNS results is observed at small Sh e , whereas excellent agreement was achieved in the nonreactive case Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF]. Improving the agreement between DNS and correlation would deserve further analysis. However, for the majority of configurations, the relative difference between the modelled and simulated Sh and Sh e is less than 30%. The agreement was significantly improved by accounting for the specific features induced by advective mechanisms, thanks to using the proper Sh i and Sh e correlations. The relevance of the proposed methodology is moreover confirmed quantitatively by the quite satisfactory values of the metric differences as reported in Tab. 5.1 in terms of:

• Mean Absolute Error (MAE), the average of the differences measured between predicted and simulated values for all cases;

• Mean Absolute Percentage Error (MAPE), similar to MAE but based on the relative differences, i.e. divided by the value from DNS;

• Determination coefficient R2, that measures the proportion of the variance in the simulation data that can be explained by the prediction.

Max Difference MAE MAPE R2 Sh i DNS vs correlation 6.28 2.09 0.207% 0.815 Sh e DNS vs correlation with Eq. (5.9) 25.9 8.14 0.698% 0.105 Sh DNS vs correlation with Eq. (5.9)

3.59 0.674 0.144% 0.964

Conclusion

The influence of an external first-order chemical reaction on mass transfer around a spherical droplet in a uniform flow was studied. The simulation strategy, based on recent work by Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], relies on the direct numerical simulation (DNS) of flow and reactive transport using the IMFT inhouse code JADIM.

In the non-reacting system, DNS results highlighted the complexity of the coupling between transport phenomena (hydrodynamics) and chemistry (thermodynamic equilibrium), which in the majority of real cases (i.e. coupled transfer) leads to configurations very far from the implicit hypotheses of the mass transfer correlations widely used in chemical engineering. Generalized correlations for the internal Sh i and external Sh e Sherwood numbers, to be used in the double film model, were proposed and validated Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF] ,Godé et al. [START_REF] Godé | Corrigendum to Towards more predictive correlations for conjugated mass transfer across a droplet interface[END_REF].

178CHAPTER 5.

SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

In this work, the extension to cases where the transferred solute is consumed by a chemical reaction confirms the acceleration of the mass transfer rate provided by the chemical reaction. However, this additional coupling between chemistry and hydrodynamics further complicates the prediction of the apparent transfer rate, especially since, as the DNS results underline, when convection is moderate (e.g. at P e e = O( 10)) the stationary concentration profiles, very different from those assumed in the classical film model at Da = 0, ultimately tend towards this asymptotic behavior when the reaction is fast (Da = 100).

The film theory, Eq. (5.4), was tested on more than 1,000 reacting cases, considering i) the improved Sh i and Sh e correlations proposed in Godé et al. [START_REF] Godé | Towards more predictive correlations for conjugated mass-transfer problems encountered in current and future for solvent extraction processes[END_REF], Godé et al. [START_REF] Godé | Corrigendum to Towards more predictive correlations for conjugated mass transfer across a droplet interface[END_REF], and ii) different definitions of the enhancement factor and the Hatta modulus. The DNS results, in terms of Sherwood numbers, were used as validation data.

The generalized correlations for the internal and external Sherwood numbers, Eq. (5.31) and (5.35), associated with the enhancement factor of Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF] were shown to accurately predict the overall Sherwood number. This method definitely brings a significant improvement of the global Sherwood number compared to the original film theory. For this, the mass transfer resistance additivity law is based on the non-reactive internal and external correlations, respectively Eq. (5.31) and (5.35), with the enhancement factor, Eq. ( 5.3), computed with the "general" definition of the Hatta modulus, Ha J , Eq. (5.9), as proposed by Juncu [START_REF] Juncu | The influence of the Henry number on the conjugate mass transfer from a sphere: II -Mass transfer accompanied by a first-order chemical reaction[END_REF]. This allows to provide satisfying, and in most cases accurate, predictions of the global Sherwood number for the coupled mass-transfer problem with a first order external chemical reaction, based only on the physical and chemical properties. Liquid-liquid extraction has seen renewed interest over the past ten years with the rise of recycling processes, whether for critical materials or new energy technologies or for rejuvenating nuclear fuel treatment. Recent advances in the field of extractive molecules, as well as new environmental constraints, require revisiting this separation process, although it is already mature on an industrial scale. In this context, numerical simulation is a valuable tool to adapt existing installations and/or accelerate the deployment of new processes. However, this requires the use of robust, validated and predictive models.

In this work, we have harnessed the power of direct numerical simulation (DNS) to develop new models for simulating liquid-liquid processes. The case study considered is that of a spherical drop, flowing in an immiscible liquid phase and possibly undergoing mass (or heat) transfer.

The study focused on four main topics for which knowledge was still lacking in the literature: two purely hydrodynamic topics, respectively the modeling of history force and the study of the appearance of 3D flow bifurcations; and two topics relating to chemistry-transport couplings, respectively dedicated to coupled mass transfer without and with chemical reaction.

On droplet hydrodynamics

The history force is a transient force often neglected in two-phase flow simulations due to the difficulty of integrating it in time, moreover in the case of (viscous) droplets due to the lack of analytical expression. However, as highlighted in Chap. 2 the contribution of this force to the slip velocity, and therefrom the residence time of the droplets in liquid-liquid devices, can be significant. In this work, we extended the work of Legendre et al. [START_REF] Legendre | Basset-Boussinesq history force of a fluid sphere[END_REF] to the case of an oscillatory flow, simulating for example the oscillating flow component in pulsed columns used in liquid-liquid extraction. The contribution of the history force to the total force exerted on the droplet was studied as a function of the Stokes number (here related to the frequency of oscillations) and the viscosity ratio, enabling us to identify the flow domain in which its contribution is not negligible. In particular, the history force contribution found to be significant (typically larger than 10% of the total RMS) for a wide range of Stokes number (between 0.5 and 5) and viscosity ratio (between 0.5 and 10), where all important dynamics of turbulence induced clustering is expected. At last, the far-field oscillatory velocity conditions were replaced by a slip velocity signal predicted by a DNS study of particles in a turbulent flow, which demonstrated the relevance of the proposed history force kernel under conditions more representative of those encountered in industrial contactors, operating in turbulent regime 3D bifurcations and their impact on the drag coefficient remain an open and important problem to correctly represent the slip velocity of droplets in an immiscible liquid. In Chap. 3, we observed that two types of bifurcation could occur by increasing the Re: external or internal bifurcation, depending on the phase in which it originates.

ON COUPLED MASS TRANSFER FROM A DROPLET

Regardless of where the bifurcation begins, both phases can be significantly impacted and we also found that depending on the Reynolds number, the bifurcation could be stable or unstable, with a different effect on drag and lift. A correlation was proposed to calculate the real drag coefficient from that predicted by the axisymmetric hypothesis, and a criterion was given for the appearance of 3D flow. At last, an original experimental study was carried out with carefully chosen fluids to verify the sphericity hypothesis.

The results confirmed the appearance of 3D bifurcation, with good DNS/experiment agreement both on the bifurcation criterion and on the drag coefficient obtained.

On coupled mass transfer from a droplet

The mass transfer of an undeformable spherical droplet in a uniform flow was studied on the basis of a detailed understanding of the characteristics of internal and external flows encountered in a wide range of Reynolds and Péclet number, with the axisymmetric flow assumption.

In Chap. 4, state-of-the-art correlations for internal and external Sherwood numbers have been revisited through an extensive DNS study. Typical coupled transfer phenomena, such as the accumulation of solutes at the core of Hill vortices and/or in the wake of droplets, that these correlations cannot predict have been identified. The original correlations have been corrected accordingly, yielding new correlations which, when used in the additivity law of mass transfer resistance, significantly improve the quality of the transferred flux prediction, which is valuable for the development of new liquid-liquid extraction processes.

The study was extended in Chap. 5 to cases where a first order chemical reaction occurs in the continuous phase. DNS results confirmed the enhancement effect provided by the chemical reaction and highlighted the mechanisms involved. The available expressions of the enhancement factor provided by the film theory, assuming either plane or spherical interface, and accounting or not for the convection contribution to mass transfer in the Hatta criteria, have been compared. We selected the best combination, which, based on the non-reactive Sherwood number from Chap. 4, made it possible to predict in a large majority of configurations, the reactive coupled transfer flow with an accuracy greater than 30%.

Towards typical conditions of industrial implementation

In this work, and using DNS, we provide new insights into droplet hydrodynamics in liquid-liquid systems and on the validity range of typical assumptions made to predict the slip velocity and mass transfer rates. More predictive correlations are also given for the prediction of mass transfer in the general case where convection plays a significant 182CHAPTER 6. CONCLUSION AND OUTLOOK role in the transport of species, and where the solute transferred is consumed by a first order reaction.

Although this thesis provides a new building block in the modeling of multiphase processes, and in particular liquid-liquid extraction, many assumptions still need to be overcome to be fully representative of industrial process implementation case conditions. Among the open questions concerning liquid-liquid extraction column simulation and modeling, the effect of the droplet deformations and surface oscillations, observed e.g. consecutive to drop's bouncing on packing elements, and of the high-volume fractions of droplets, up to 20 or even 30%, prevailing in solvent extraction contactors still need to be quantified.

The collective behaviour of droplets in liquid-liquid extraction columns can have a significant impact on the efficiency of the extraction process: neighbouring droplets can collide, hence changing their trajectory and velocity, or even coalesce, hence modifying the exchange surface ; additional turbulence and mixing may be induced by the local restrictions of the fluid vein between the particles, in their wake, etc. The literature on this topic proposes two approaches to tackle the collective effects: the first one being based on physical modeling, see Yu et al. [START_REF] Yu | Mass transfer and droplet behaviors in liquid-liquid extraction process based on multiscale perspective: A review[END_REF] and the second one by taking advantage of machine learning, see Chouai et al. [START_REF] Chouai | Use of neural networks for liquid-liquid extraction column modelling: An experimental study[END_REF].

Droplet's deformations may impact both hydrodynamics and the mass transfer rate. They are usually addressed by VOF, see Asghar et al. [START_REF] Asghar | Numerical wetting benchmarks -advancing the plicrdf-isoadvector unstructured volume-offluid (vof) method[END_REF], or other Front Tracking methods Bois [START_REF] Bois | offluid/level-set algorithm of triocfd[END_REF]. However, these methods, although they can be convenient for hydrodynamics alone, may suffer from lack of resolution for the case of mass transfer, where a steep concentration gradient prevails at the deforming interface. Rachih [START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF] simulate the mass transfer from an ellipsoidal droplet with JADIM. These preliminary developments, relating to non-spherical but still non-deformable droplets, should make it possible to understand the impact of the shape of the droplets on their motion and Sherwood number, and to 1st order, to cover a wider range of droplet shapes encountered. in an extraction column. Deformations through shape oscillations, as studied in Piedfert et al. [START_REF] Piedfert | Numerical simulations of a rising drop with shape oscillations in the presence of surfactants[END_REF], have a significant impact on internal recirculation and therefore should also have an impact on mass transfer. A process enhancement method to be studied is to set the extraction column pulsation at a precise droplet resonance frequency.

3D bifurcations 3D bifurcations significantly modify the flow around and/or in the drop, which results in a notable change in drag forces (Chap. 3). The 3D flow that takes place at high Re is also likely to impact the kinetics of mass transfer. The internal bifurcation, which results in an opening of the streamlines of the 2D Hill vortex, at the heart of which an accumulation of solute was observed, should favor the internal transfer towards the interface of the drop. Conversely, the External bifurcation can reinforce the "reservoir" effect observed at the rear of the droplet in 2D cases. Additional simulations were carried out to understand the effect of 3D bifurcations on the Sherwood number. In this aim, we used the same strategy as in Chap. 4, i.e. by considering the transport of a solute in a frozen flow, however using the 3D mesh. We therefore limited the analysis to the case and P e = 1000. The external bifurcation turns out to have a negligible effect on Sh i and a slight effect on Sh e , and hence on Sh. In the case considered, where the solute transfer is favored by thermodynamics (k = 10), we observe the same trends as with the 2D assumptions, with a more pronounced "reservoir" effect as P e e increases.

These simulations, while confirming the effect of bifurcations on mass transfer, show that their interpretation, and a fortiori their correlation to Sh, is not straightforward and merits further study, for which the model developed in this thesis, with JADIM, is well suited. However, and even if only 2 examples have been treated, which it would be inappropriate to generalize, it appears that the effect is less spectacular than expected, and that the 2D correlations proposed in Chap. 4 remains relevant as a first approximation.

Approximation of the kernel to avoid overflow error

History force has often been neglected due to the computational cost for the integration of the memory kernel in equation (2.5). For droplets, the problem is even worse because it involves multiplications between an exponential and a complementary error function. Indeed, each term can easily overpass the maximum number of digits a computer can store, even if they compensate each over (one tends to infinity when the other tends to zero). A solution is to approximate the complementary error function, with an exponential function with a negative power term, therefore compensating the exponential one: Abramowitz et al. [START_REF] Abramowitz | Handbook of mathematical functions with formulas, graphs, and mathematical tables[END_REF]. Considering the following function f : We show in Fig. 7.1, that the approximation, precisely, respects the exact solution.

f = e y × erf c(x), ( 7
The expression of the history force with this simplification is then If the analytic expression of dW dt ′ is known, it is possible, in some cases, to integrate this expression exactly and get rid of the numerical integration. But calculating the history force in general cases needs to compute dW dt ′ along the trajectory. We cannot solve the integral exactly, but the reader can refer to the work of Daitche [START_REF] Daitche | Advection of inertial particles in the presence of the history force: Higher order numerical schemes[END_REF] for more details.

F h (t) = 6πµ e R

Derivation of the analytical solution of the reactiondiffusion problem in a spherical domain

For the non-reactive transfer, we use the fact that in spherical coordinate △C = 0, thus we obtain Eq. (7.5). For the reactive transfer, we make the following change of variables:

C * = C Cs = f (r)
r with C * a normalized concentration, C s the interfacial concentration and f (r) an introduced function, with some calculation steps we obtain Eq. (7.6). If we go back to the definition of the enhancement factor, but this time defined with the transfer rate, see equation (7.11). With the precedent results, we can deduce a solution of the enhancement factor in pure diffusive spherical configuration, see Eq. (7.12).

1. 1

 1 Mass transfer according to film theory scheme. . . . . . . . . . . . . . . 1.2 Small scale pulsed column used for R&D study in nuclear environment at CEA Marcoule, from left to right: column, droplet swarm, droplet interaction and isolated droplet mass transfer schematic. . . . . . . . . 1.3 Clift diagram, relating droplet shape to physical properties. . . . . . . . 1.4 Hill's vortex visualization. . . . . . . . . . . . . . . . . . . . . . . . . . 2.1 Time evolution of the forces acting on a droplet over one period of a uniform oscillatory flow (t * = tf ): Total force, Steady drag force F SD , Inertial or pressure gradient force & . . . . . . . . . . . . . . . 2.2 Visualization of the streamlines inside and around a viscous droplet at different times of the oscillation period, for Re = 0.1, µ * = 0.2 and f * = 1.

8 LIST OF FIGURES 2 . 8 61 3. 1

 828611 2.2f). . . . . . . . . . . . . . . . . 2.4 Visualization of the streamlines over the oscillation period for Re = 0.1, µ * = 5 and f * = 1. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.5 Visualization of the streamlines over one period for Re = 0.1, µ * = 0.2 and f * = 0.01. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.6 Visualization of the streamlines over one period for Re = 0.1, µ * = 0.2 and f * = 10 3 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.7 (a) Evolution of the normalized tangential velocity U θ /W 0 as a function of the polar angle θ at the interface. (b) Radial evolution of the normalized axial velocity U x (x = 0, r)/W 0 . The times considered are the same as in Fig. 2.6: Fig. 2.6a, Fig. 2.6b, Fig. 2.6c, Fig. 2.6d, Fig. 2.6e, Fig. 2.6f. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Evolution of all RMS contributions with f * . Total force RMS, Steady force RMS, Inertial (Pressure gradient) and . . . . . . . . . . . . . . 49 2.9 Evolution of each force RMS relative to the total force RMS with f * . Steady drag force RMS, Inertial (Pressure gradient) and . . . . . 50 2.10 Normalized transfer function F = Fh /6πµ e R W (f * ) as a function of f * for different viscosity ratios: (a) µ * = 0.05, (b) µ * = 0.2 and (c) µ * = 20. Numerical simulations at Re = 0.1. Solid sphere (relation (2.14), Spherical bubble (relation (2.15)), Fluid particle of viscosity µ * (relation (2.13)), relation (2.27) with slip length λ = R/3µ * (relation 2.4), relation( 2.27) with unsteady slip length given by relation (2.24), and dashed line asymptotic value of our kernel when f * → ∞ . . . . . 51 2.11 Relative contribution of the history force RMS to the total force in the phase diagram Stµ * . Color contours are created using Delaunay triangulation between our simulation data (black dots). . . . . . . . . . 53 2.12 (a) Interfacial slip length λ/R as a function of θ at t=0.25T, for: Re = 0.1, µ * = 0.2 and f * = 0.01, Re = 0.01, µ * = 0.05 and f * = 1, Re = 0.1, µ * = 5 and f * = 10., (b) Interfacial slip length λ/R as a function of time over a period, for Re = 0.1, µ * = 0.2 and f * = 0.01 . . 54 2.13 Normalized history force as a function of the normalized time t/T . Comparison between direct simulations ( ) and equation (2.5) ( ) with the kernel given by Eq. (2.5), the slip length by relation( 2.4) and ( ) with the kernel given by Eq. (2.5) and the slip length by relation( 2.24) 56 2.14 Stokes and viscosity ratio map of (left) History force model M AE * , (right) Total force model M AE * , with slip length λ = R/3µ * . . . . . . . . . . 57 2.15 Stokes and viscosity ratio map of (left) History force model M AE * , (right) Total force model M AE * , with slip length expression( 2.24). . . . . . . 57 2.16 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 59 2.17 Comparison of the time evolution of the history force from simulations with the model for different viscosity ratios. . . . . . . . . . . . . 60 2.18 Stokes number and viscosity ratio map of (left) History force model M AE * , (right) Total force model M AE * , with slip length expression( 2.24), with result of simulation in a "turbulent" flow. . . . . . . . . . . . . . . . . . 60 2.19 Comparison of the time evolution of the total force from simulations with the model for different viscosity ratios. . . . . . . . . . . . . Experimental signature of a nonlinear trajectory. ⃝: droplet, : trajectory, →: instantaneous velocity vector and →: instantaneous acceleration vector. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68 3.2 Picture (left) and scheme (right) of the experimental setup. . . . . . . . 71 3.3 Clift diagram with maximum Reynolds and radius deduced for n-Pentane , Heptane and Cyclohexane . . . . . . . . . . . . . . . . . . . . 72 3.4 Example of raw image (left) and corresponding post-processed image with droplet center and radius detection (right) . . . . . . . . . . . . . 73 LIST OF FIGURES 3.5 Description of the meshing procedure and zoom on the mesh in the vicinity of the droplet (the actual simulation domain is 50 times larger than the droplet radius). . . . . . . . . . . . . . . . . . . . . . . . . . . 3.6 Evolution of the drag coefficient as a function of the Reynolds number.

  (line) color corresponds to µ * : µ * = 20, µ * = 10, µ * = 5, µ * = 3, µ * = 2, µ * = 0.5, µ * = 0.33, µ * = 0.2, µ * = 0.1, µ * = 0.05. . . . 3.7 Trajectory (top figures) and time evolution of the aspect ratio in experiments (bottom figures) of Heptane or n-Pentane droplets in water. . . . 3.8 Time evolution of the velocity for experimental cases, velocity modulus raw data , Y-velocity raw data, X-Velocity raw data, and rolling mean data. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.9 Flow map based on experimental observations of droplets rising in water, straight vertical trajectory, oscillating trajectory. From left to right: n-Pentane µ * = 0.231, Heptane µ * = 0.416 and Cyclohexane µ * = 1.051. 3.10 comparison of drag coefficient Reynolds number evolution of experimental cases , with axisymmetric simulation data . . . . . . . . . . . . . . . 3.11 Experimental mean drag coefficient C d and comparison with corresponding axisymmetric simulations . Experimental results are reported with an error bar estimated from Eq. (3.20). . . . . . . . . . . . . . . . . . . 3.12 Examples of axisymmetric flows. Left: low-viscosity droplet with low Re; right: viscous droplet with moderate Re, showing detachment of the outer boundary layer. Top: xy plane at z = 0; Middle: xz plane at y = 0; Bottom: yz plane at x = 0. . . . . . . . . . . . . . . . . . . . . . . . . 3.13 Vortex core inside and outside the droplet, for the cases shown in Fig. 3.13 3.14 Comparison of internal and external steady 3D flow cases. Top: xy plane at z = 0; Middle: xz plane at y = 0; Bottom: yz plane at x = 0. . . . . 3.15 Comparison of internal and external bifurcation consequences for the 3D steady flows of Fig. 3.14. . . . . . . . . . . . . . . . . . . . . . . . . . 3.16 Examples of unsteady 3D flow cases. . . . . . . . . . . . . . . . . . . . 3.17 Comparison of unsteady 3D flow cases. . . . . . . . . . . . . . . . . . . 3.18 Determination of the bifurcation origin from single phase flow simulations and visualization for µ * = 5, ρ * = 1 and Re = 350 . . . . . . . . . . . . 3.19 Determination of the bifurcation origin from single phase flow simulations and visualization for µ * = 02, ρ * = 1 and Re = 200 . . . . . . . . . . . 3.20 Two identical meshes and simulation settings (µ * = 0.1, ρ * = 1, and Re = 50), performed with two initial random perturbations. The streamlines only inside the droplet are presented to highlight the orientation of symmetry planes. The flow direction follows the x-axis which is orthogonal to the plane of visualization. . . . . . . . . . . . . . . . . . . . . . . . .LIST OF FIGURES 3.21 Evolution of maximum azimuthal velocity w max and damping or amplification factor α as a function of time. The grey region represents the plateau where the damping or amplification factor is determined. . 93 3.22 Evolution of amplification factor in function of Reynolds number or density ratio. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94 3.23 Flow phase diagram for ρ * = 1, axisymmetric cases, steady internal bifurcation marginal state, steady 3D flow state with internal bifurcation, 3D unsteady internal bifurcation, unsteady internal bifurcation marginal state, steady external bifurcation marginal state, steady 3D flow state with external bifurcation, 3D unsteady flow with external bifurcation, unsteady external bifurcation marginal state. . . . . . . 95 3.24 Phase diagram of external bifurcation with and without the effect of density ratio, axisymmetric cases, ρ * ̸ = 1 axisymmetric cases, external bifurcation marginal state, external bifurcation marginal curve, external steady bifurcation cases, external bifurcation unsteady cases, external bifurcation unsteady marginal curve, ρ * ̸ = 1 external bifurcation cases. 97 3.25 Internal bifurcation map at ρ * = 1, axisymmetric cases, internal steady bifurcation critical point, internal steady bifurcation critical line, steady state internal bifurcation cases, internal unsteady bifurcation cases, internal unsteady critical bifurcation cases. . . . . . . . . . . . 98 3.26 Internal bifurcation map with density ratio parametric study simulations, axisymmetric cases, ρ * ̸ = 1 axisymmetric cases, internal steady bifurcation critical point, internal steady bifurcation critical line, internal steady bifurcation cases, internal unsteady bifurcation cases, internal unsteady bifurcation critical line, ρ * ̸ = 1 internal bifurcation cases. . . 99 3.27 Time evolution of drag coefficient in function of time for simulations presented in Sec. 3.6.1. . . . . . . . . . . . . . . . . . . . . . . . . . . . 100 3.28 3D drag coefficient divided by the axisymmetric drag coefficient in function of the Reynolds number, simulation results and linear fit . . . . . . 101 3.29 3D drag coefficient divided by the axisymmetric drag coefficient in function of the viscosity ratio, simulation results and regression . . . . . . . 102 3.30 Comparison of density ratio effect on drag coefficient for internal and external bifurcations. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 103 3.31 Strouhal number evolution with Reynolds number for simulation results. 104 3.32 Bifurcation map from experimental data, straight trajectory, oscillating trajectory, internal steady bifurcation critical line, and internal steady bifurcation critical line. . . . . . . . . . . . . . . . . . . . . . . . 105 3.33 Internal bifurcation map with experimental data, axisymmetric cases, internal steady bifurcation critical line, steady state internal bifurcation cases, internal unsteady bifurcation cases, internal unsteady critical bifurcation cases. . . . . . . . . . . . . . . . . . . . . . . . . . . 106 3.34 Drag coefficient from the simulations ( (axisymetric) and 3D steady and for unsteady) and experimental results with error bar. . . . . . . . . 107 LIST OF FIGURES 3.35 Evolution of Strouhal number with Reynolds number, experimental results for a spherical droplet freely rising in water. . . . . . . . . . . . . . . . 3.36 Strouhal number comparison of simulations , and experimental results .4.1 Illustration of possible mass-transfer regimes in the film theory. . . . . 4.2 Detail of the domain discretization: The droplet, with a polar mesh, is in the center. The green and red lines, extending respectively on -1 < x/R < 2 and 0 < y/R < 2 in the parallel and perpendicular directions of the external flow, materialize the lines used to plot the concentration profiles in Sec. 4.3. . . . . . . . . . . . . . . . . . . . . . 4.3 Possible configurations encountered in solvent extraction processes. The arrows indicate the direction of the transfer: red for extraction, green for stripping. Blue (resp. yellow) stands for the aqueous phase (resp. organic): (a) and (c) represent O/W configurations (for which µ * > 1), (b) and (d) W/O configurations (µ * < 1). . . . . . . . . . . . . . . . . . 4.4 Time evolution of the global (•), internal ( * ), and external (△) Sherwood numbers in the case of a solute transferred from the continuous phase to the droplet (a) (initial concentrations C i = 0, C e = 1, black lines) and from the droplet to the continuous phase (c) (initial concentrations C i = 1, C e = 0, red symbols). In both cases: Re = 20, µ * = 5, k = 10, P e e = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.5 Streamlines and concentration contours in the case "Oil in Water" (µ * = 5). Effect of the Reynolds and external Péclet numbers for k = 0.1 (k √ D * = 0.04). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.6 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.5. "Oil-in-Water", k √ D * = 0.04 (Symbols: Re = 20, lines Re = 200). . . . . . . . . . . . . . . . . . 4.7 Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the Reynolds and external Péclet numbers for k = 0.1 (k √ D * = 0.22). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.8 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.7. "Water-in-Oil", k √ D * = 0.22. (Symbols: Re = 20, Lines: Re = 200). . . . . . . . . . . . . . . . 4.9 Streamlines and concentration contours in the case "Oil in Water" (µ * = 20). Effect of the Reynolds and external Péclet numbers for k = 1 (k √ D * = 0.22). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.10 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.9. "Oil-in-Water", k √ D * = 0.22 with k = 1. (Symbols: Re = 20, Lines: Re = 200). . . . . . . . . 4.11 Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the Reynolds and external Péclet numbers for k = 10 (k √ D * = 22). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . LIST OF FIGURES 4.12 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.11. "Water-in-Oil", k √ D * = 22.36. (Symbols: Re = 20, Lines: Re = 200). . . . . . . . . . . 136 4.13 Streamlines and concentration contours in the case "Oil in Water" (µ * = 5). Effect of the Reynolds and external Péclet numbers for k = 10 (k √ D * = 4.47). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 137 4.14 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.13. "Water-in-Oil", k √ D * = 4. (Symbols: Re = 20, Lines: Re = 200). . . . . . . . . . . . 138 4.15 Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.05). Effect of the Reynolds and external Péclet numbers for k = 1 (k √ D * = 4.47). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 139 4.16 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.15. "Water-in-Oil", k √ D * = 4.47 with k = 1 (Symbols: Re = 20, Lines: Re = 200). . . . . 139 4.17 Parity plots of isolated problem correlation results compared to coupled problem simulations results. . . . . . . . . . . . . . . . . . . . . . . . . 142 4.18 Comparison between maximum interfacial velocity U max correlation fitted for droplet in function of µ * (line) and DNS results (points). . . . . . . 143 4.19 Asymptotic values of the internal Sherwood number, Sh i , at low (left) and high (right) Péclet. . . . . . . . . . . . . . . . . . . . . . . . . . . . 145 4.20 Comparison between Sh E Eq. (4.27) (line) and DNS results Sh e (points) for coupled problems at Re = 10. k = 0.1, k = 0.25, k = 1, k = 4 and k = 10. . . . . . . . . . . . . . . . . . . . . . . . . . . 147 4.21 Evolution of simulation Sh e with k for bubbles (µ * = 0) at Re = 10 and P e = 100 , and comparison with value from Eq. (4.28) for bubble, Sh E (black dashed line). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 148 4.22 Parity plots of new correlation results compared to coupled problem simulations results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 150 5.1 Film theory scheme for conjugate mass transfer. . . . . . . . . . . . . . 155 5.2 Simulation configuration at initial time with streamlines. . . . . . . . . 158 5.3 Simulation domain and mesh (left) and validation test (right). (a) Domain discretization, zoom on the droplet. The red and green lines are used in Sec. 5.3 to plot the axial and perpendicular concentration profiles, respectively. (b) Comparison between the enhancement factor E deduced from DNS (symbols) and the analytical solution Eq. (7.12) (line) for different Da values. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 160 5.4 Simulation strategy for the mass transfer investigation based on avalaible hydrodynamic simulations from Godé et al. [5] . . . . . . . . . . . . . 162 LIST OF FIGURES 5.5 Time evolution of the Sherwood numbers obtained by DNS for Re = 20, D * = 0.2, k = 10 and P e e = 10 (with F o = Det R 2 the dimensionless time). Effect of an external first order chemical reaction with Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). Steady external Sherwood values are reported in Tab. 7.2. . . . . . . . . . . . . 5.6 Time evolution of the mean concentration in the droplet (line) at the droplet interface (point): effect of an external first order chemical reaction with Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue) (Re = 20, D * = 0.2, k = 10 and P e e = 10). Concentrations are normalized by the mean concentration in the droplet. . . . . . . . 5.7 Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 0.2, k = 10 and P e e = 10 . . . . . . . . . . . . . . . . . 5.8 Concentration profiles for Re = 20, D * = 0.2, k = 10 and P e e = 10 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet. . . . . . . . . 5.9 Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 0.2, k = 10 and P e e = 1000 . . . . . . . . . . . . . . . . 5.10 Concentration profiles for Re = 20, D * = 0.2, k = 10 and P e e = 1000 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet. . . . . . . . . 5.11 Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 5, k = 0.1 and P e e = 10 . . . . . . . . . . . . . . . . . . 5.12 Concentration profiles for Re = 20, D * = 5, k = 0.1 and P e e = 10 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet. . . . . . . . 5.13 Concentration spatial distribution and streamlines for different Da at Re = 200, D * = 5, k = 0.1 and P e e = 1000. . . . . . . . . . . . . . . . 5.14 Enhancement factor: comparison of simulation results (symbols, colored from light blue for Da = 1 to dark blue for Da = 100) with predictions of the simplified models (solid line: E F T , dash line: E SF T ). (a) Hatta from film theory, (b) General expression of Hatta from Juncu [6]. . . . 5.15 Enhancement factor: comparison of simulation results (symbols, colored from light blue for P e e = 10 to dark blue for P e e = 1000) with predictions of the simplified models (continuous line: E F T , dashed line: E SF T . (a) Hatta from classic film theory, (b) General expression of Hatta from Juncu [6]. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 5.16 Parity plots of external, internal and global Sherwood numbers obtained from DNS simulation and model prediction. In each case, the darker the symbol, the greater the value of the parameter, i.e. D * in (a), k in (b), Sh i in (c) and Sh e in (d). The straight lines reprensent the limits of 0% error (black), 10% error (blue) and 30% error (red). . . . . . . . . . . . 176 6.1 Evolution of Sh e , Sh i , and Sh with P e e at steady-state, in the case of an internal bifurcation (Re = 50, µ * = 0.1, and k = 0.1), and comparison with the corresponding axisymmetric predictions (dashed line)183 6.2 Streamlines and concentration field visualization in the case of an internal bifurcation at Re = 50, µ * = 0.1, k = 0.1, at small (left) and elevated (right) external Péclet. . . . . . . . . . . . . . . . . . . . . . . . . . . . 184 6.3 Evolution of Sh e , Sh i , and Sh with P e e at steady-state, in the case of an external bifurcation (Re = 500, µ * = 5, and k = 10), and comparison with the corresponding axisymmetric predictions (dashed line)185 6.4 Streamlines and concentration field visualization in the external bifurcation at Re = 500, µ * = 5, k = 10, comparison at different Peclet number. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 185 7.1 Comparison of kernel computation with or without the approximation shown in Eq. 7.3. The curve without approximation is computed using the python library mpmath to ensure accuracy and avoid overflow errors. 189 7.2 Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the external Reynolds and Péclet numbers for k = 0.1 (k √ D * = 0.22). . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 7.3 Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.7a and 7.2 "Water-in-Oil", Re = 20, µ * = 0.2, k √ D * = 0.22 with k = 0.1 (Lines: P e = 10, Symbols: P e = 100. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 191 7.4 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.1. "Oil-in-Water", k √ D * = 0.04. . . . . . . . . . . . . . . . . . . . . . . . . . . . 192 7.5 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.1. "Waterin-oil", k √ D * = 0.22. . . . . . . . . . . . . . . . . . . . . . . . . . . . 193 7.6 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.1. "Oil-in-Water", k √ D * = 0.22 with k = 1. . . . . . . . . . . . . . . . . . . . . 194 7.7 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Waterin-Oil", k √ D * = 22.36. . . . . . . . . . . . . . . . . . . . . . . . . . . 195 LIST OF FIGURES 7.8 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Oil-in-Water", k √ D * = 4.47. . . . . . . . . . . . . . . . . . . . . . . . . . . . 7.9 Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Waterin-Oil", k √ D * = 4.47 with k = 1. . . . . . . . . . . . . . . . . . . . . . 7.10 Time evolution of the mean internal concentration in the droplet (line) and interfacial (point), for Re = 20, µ * = 5, k = 10 and P e e = 1000 at Da : Da = 0, Da = 1, Da = 10, Da = 100 . . . . . . . . . . . . . . 7.12 Time evolution of the mean internal concentration for Re = 200, µ * = 0.2, k = 0.1 and P e e = 1000 at Da : Da = 0, Da = 1, Da = 10, Da = 7.11 Concentration evolution for Re = 200, µ * = 0.2, k = 0.1 and P e e = 1000 at Da : Da = 0, Da = 1, Da = 10, Da = 100 . . . . . . . . . . . .

1 . 1 1 . 2 Implementation and modelling of liquid-liquid extraction 1 . 3 1 . 4

 11121314 An introduction to PUREX process. . . . . . . . . . . . . . Focus on isolated spherical droplet . . . . . . . . . . . . . . Thesis plan . . . . . . . . . . . . . . . . . . . . . . . . . . . .

28CHAPTER 1 .Figure 1 . 1 :

 111 Figure 1.1: Mass transfer according to film theory scheme.

Figure 1 . 2 :

 12 Figure 1.2: Small scale pulsed column used for R&D study in nuclear environment at CEA Marcoule, from left to right: column, droplet swarm, droplet interaction and isolated droplet mass transfer schematic.

4 )

 4 With g the gravity acceleration, µ the dynamic viscosity (with index e stating if we use the exterior viscosity, i the interior), ρ the density, σ the surface tension, U the velocity, and r the droplet radius. This graph is called the Clift diagram and is shown in Fig.1.3.

Figure 1 . 3 :

 13 Figure 1.3: Clift diagram, relating droplet shape to physical properties.

Fig. 2 .

 2 Fig. 2.1 compares the magnitude of these forces over a period of the oscillating flow for two specific cases: Re = 0.01, µ * = 5, f * = 0.1 (f * i = 0.02), and Re = 0.1, µ * = 5, f * = 100 (f * i = 20). All forces are scaled by the mean steady drag force over the period:

Figure 2 . 1 :

 21 Figure 2.1: Time evolution of the forces acting on a droplet over one period of a uniform oscillatory flow (t * = tf ): Total force, Steady drag force F SD , Inertial or pressure gradient force & Added-mass force F I&AM , History force F H . (Left): Re = 0.01, µ * = 5 and f * = 0.1 (f * i = 0.02). (Right): Re = 0.1, µ * = 5 and f * = 100 (f * i = 20). All forces are scaled by the mean steady drag force over the period F 0 given by Eq. 2.20.

44CHAPTER 2 .

 2 BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW

(

Figure 2 . 2 :

 22 Figure 2.2: Visualization of the streamlines inside and around a viscous droplet at different times of the oscillation period, for Re = 0.1, µ * = 0.2 and f * = 1. More instants are shown during the second half of the period to provide a detailed description of the two recirculation vortices that develop both inside and outside the droplet during the deceleration phase.

Figure 2 . 3 :

 23 Figure 2.3: (a) Evolution of the normalized tangential velocity U θ /W 0 with the polar angle θ at the interface. (b) Radial evolution of the normalized axial velocity U x (x = 0, r)/W 0 . The times considered are the same as in Fig. 2.2: t = 0T , W (t) = W 0 (Fig. 2.2a), t = 0.554T , W (t) = 0.667W 0 (Fig. 2.2b), t = 0.582T , W (t) = 0.507W 0 (Fig. 2.2c), t = 0.670T , W (t) = 0.124W 0 (Fig. 2.2d), t = 0.800T , W (t) = 0.0489W 0 (Fig. 2.2e), t = 0.883T , W (t) = 0.329W 0 (Fig. 2.2f). .

46CHAPTER 2 .Figure 2 . 4 :

 224 Figure 2.4: Visualization of the streamlines over the oscillation period for Re = 0.1, µ * = 5 and f * = 1.

  2.5). For f * = 10 3 the flow structure around and inside the droplet cannot follow the fast dynamics of the external flow. This is observed in Fig. 2.7 where the radial profile of the external flow shows no change inside the drop. The transition from external to internal flow occurs across a very thin diffusion layer. However, the development of two vortices is still observed in the outer region, as shown in Fig. 2.6e.

Figure 2 . 5 :Figure 2 . 6 :

 2526 Figure 2.5: Visualization of the streamlines over one period for Re = 0.1, µ * = 0.2 and f * = 0.01.

Figure 2 . 7 :

 27 Figure 2.7: (a) Evolution of the normalized tangential velocity U θ /W 0 as a function of the polar angle θ at the interface. (b) Radial evolution of the normalized axial velocity U x (x = 0, r)/W 0 . The times considered are the same as in Fig. 2.6: Fig. 2.6a, Fig. 2.6b, Fig. 2.6c, Fig. 2.6d, Fig. 2.6e, Fig. 2.6f.

Figure 2 . 8 :

 28 Figure 2.8: Evolution of all RMS contributions with f * . Total force RMS, Steady force RMS, Inertial (Pressure gradient) and Added-mass forces RMS, History force RMS. (Left): simulation results for Re = 0.01 and µ * = 20. (Right): simulation results for Re = 0.1 and µ * = 0.05.

Figure 2 . 9 :

 29 Figure 2.9: Evolution of each force RMS relative to the total force RMS with f * . Steady drag force RMS, Inertial (Pressure gradient) and added-mass forces RMS, History force RMS, ---Limit 10% contribution. (Left): simulation results with Re = 0.01 and µ * = 5. (Right): simulation results with Re = 0.1 and µ * = 0.05.

2. 6 .

 6 HISTORY FORCE IN THE FREQUENCY DOMAIN 51 (a) µ * = 0.05 (b) µ * = 0.2 (c) µ * = 20

Figure 2 . 10 :

 210 Figure 2.10: Normalized transfer function F = Fh /6πµ e R W (f * ) as a function of f * for different viscosity ratios: (a) µ * = 0.05, (b) µ * = 0.2 and (c) µ * = 20. Numerical simulations at Re = 0.1. Solid sphere (relation (2.14), Spherical bubble (relation (2.15)), Fluid particle of viscosity µ * (relation (2.13)), relation (2.27) with slip length λ = R/3µ * (relation 2.4), relation( 2.27) with unsteady slip length given by relation (2.24), and dashed line asymptotic value of our kernel when f * → ∞ .

.26) 2 . 7 . 53 Fig. 2 .

 27532 Fig.2.11 shows the history force RMS reported in the phase diagram St-µ * for St varying between 10 -3 to 10 +3 and µ * between 0.05 and 20. The history force RMS contribution is found to be more than 10 %, for Stokes number in the range 0.01 < St < 5. Similar results have been observed for radial migration of particles in a vortex for St < 0.01 with a significant history force contribution Candelier et al.[START_REF] Candelier | On the effect of the boussinesq-basset force on the radial migration of a stokes particle in a vortex[END_REF]. At lower Stokes numbers, the history force RMS is almost equal to the steady drag RMS contribution, while at high frequency the pressure gradient and added-mass force are the dominant forces. However, the range of parameters (Stokes number and viscosity ratio) over which the history force must be considered is wide and corresponds to the regime where preferential accumulation is expected (0.05 < St < 5) for droplets (0.5 < µ * < 10). Therefore, an accurate expression for the history kernel able to calculate the unsteady force for Lagrangian tracking is highly desirable.

10 ¡2 10 ¡1 10 Figure 2 . 11 :

 1010211 Figure 2.11: Relative contribution of the history force RMS to the total force in the phase diagram Stµ * . Color contours are created using Delaunay triangulation between our simulation data (black dots).

Figure 2 .

 2 Figure 2.12: (a) Interfacial slip length λ/R as a function of θ at t=0.25T, for: Re = 0.1, µ * = 0.2 and f * = 0.01, Re = 0.01, µ * = 0.05 and f * = 1, Re = 0.1, µ * = 5 and f * = 10., (b) Interfacial slip length λ/R as a function of time over a period, for Re = 0.1, µ * = 0.2 and f * = 0.01

56CHAPTER 2 .

 2 BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW Re = 0.01, µ * = 5, f * = 0.1 Re = 0.1, µ * = 0.2, f * = 1 Re = 0.01, µ * = 0.05, f * = 1000

Figure 2 . 13 :

 213 Figure 2.13: Normalized history force as a function of the normalized time t/T . Comparison between direct simulations ( ) and equation (2.5) ( ) with the kernel given by Eq. (2.5), the slip length by relation( 2.4) and ( ) with the kernel given by Eq. (2.5) and the slip length by relation( 2.24) . (a) Re = 0.01, µ * = 5, f * = 0.1 (f * i = 0.02), r 2 = 0.9619, (b) Re = 0.1, µ * = 0.2, f * = 1 (f * i = 5), r 2 = 0.9452, (c) Re = 0.01, µ * = 0.05, f * = 100 (f * i = 2000), r 2 = 0.1970.

10 ¡2 10 ¡1 10 Figure 2 . 14 :Figure 2 . 15 :

 1010214215 Figure 2.14: Stokes and viscosity ratio map of (left) History force model M AE * , (right) Total force model M AE * , with slip length λ = R/3µ * (expression 2.4). Contour is created using (Delaunay et al. [51]) triangulation between data points.

60CHAPTER 2 .

 2 BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW (a) µ * = 0.2 (b) µ * = 0.5 (c) µ * = 2 (d) µ * = 5

Figure 2 . 17 :Figure 2 . 18 :

 217218 Figure 2.17: Comparison of the time evolution of the history force from simulations with the model for different viscosity ratios.

2. 9 .

 9 OUTLOOK, ACCURACY OF HISTORY FORCE KERNEL IN A TURBULENT FLOW61 (a) µ * = 0.2 (b) µ * = 0.5 (c) µ * = 2 (d) µ * = 5

Figure 2 . 19 :

 219 Figure 2.19: Comparison of the time evolution of the total force from simulations with the model for different viscosity ratios.

3 70 3. 3 . 1 70 3. 3 . 2 72 3. 4 74 3. 4 . 1 76 3. 4 . 3 79 3. 5 79 3. 5 . 1 79 3. 5 . 2 81 3. 5 . 3 82 3. 6 84 3. 6 . 1 63 64CHAPTER 3 . 3D BIFURCATIONS STUDY 3 . 7 3 . 8

 3703170327247441764379579517952815382684616333738 Experimental methods . . . . . . . . . . . . . . . . . . . . . Experimental setup . . . . . . . . . . . . . . . . . . . . . . . . Fluids selection . . . . . . . . . . . . . . . . . . . . . . . . . . 71 3.3.3 Droplet detection and post-processing . . . . . . . . . . . . . Numerical simulations . . . . . . . . . . . . . . . . . . . . . Numerical methods . . . . . . . . . . . . . . . . . . . . . . . . 74 3.4.2 Validation tests . . . . . . . . . . . . . . . . . . . . . . . . . . Outline of the parametric study . . . . . . . . . . . . . . . . . Experimental observations . . . . . . . . . . . . . . . . . . . Droplet trajectories . . . . . . . . . . . . . . . . . . . . . . . Experimental flow map of rising droplets . . . . . . . . . . . Drag coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . DNS study of the transition from 2D to 3D flow . . . . . Flow description . . . . . . . . . . . . . . . . . . . . . . . . . 84 3.6.2 Internal or external origin of instability . . . . . . . . . . . . 90 3.6.3 Flow phase diagram from numerical simulations . . . . . . . 91 3.6.4 Effect of 3D flow on drag . . . . . . . . . . . . . . . . . . . . 100 Discussion of experimental and numerical outcomes. . . . 3.7.1 Bifurcation diagram from simulations and experiments . . . . 3.7.2 Drag coefficient . . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.3 Strouhal number . . . . . . . . . . . . . . . . . . . . . . . . . 3.7.4 Statistical analysis . . . . . . . . . . . . . . . . . . . . . . . . Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 3 . 1 :

 31 Figure 3.1: Experimental signature of a nonlinear trajectory. ⃝: droplet, : trajectory, →: instantaneous velocity vector and →: instantaneous acceleration vector.

70CHAPTER 3 .

 3 3D BIFURCATIONS STUDYΦ =A e A = Surface area of volume equivalent sphere Measured surface area(3.14)

Figure 3 . 2 :

 32 Figure 3.2: Picture (left) and scheme (right) of the experimental setup.

72CHAPTER 3 .Figure 3 . 3 :

 333 Figure 3.3: Clift diagram with maximum Reynolds and radius deduced for n-Pentane , Heptane and Cyclohexane .

Figure 3 . 4 :

 34 Figure 3.4: Example of raw image (left) and corresponding post-processed image with droplet center and radius detection (right)

3. 4 2 8. 8 , 2 8. 8 )

 42828 . NUMERICAL SIMULATIONS 75 ap close = 1.1345 Expansion factor of streamline inside the droplet and outside near the droplet ac close = 1.30 Expansion factor of equipotentials close to the droplet np f ar = 10 Number of streamlines, far from the droplet nc f ar = 10 Number of equipotentials for from the droplet n ϕ = 64 Number of planes around the x-axis j max = 30 Radial discretization of the droplet nc arc = 30 Angular discretization of the quarter of a droplet lx = 50 Size of the domain in x relative to the droplet radius ly = 50 Size of the domain in y relative to the droplet radius dp min = smallest( Re -1 P e -1 Distance between the interface of the droplet and the first pressure point np close -= ln ly dp min -1 np far -1 1-ap close ln(ap close ) + 1 Number of streamlines close to the droplet nc close = ln 2 * lx * ncarc πr -1 nc far -1 1-ac close ln(ac close ) Number of equipotentials close to the droplet Table 3.2: Definition of the parameters of the mesh.particular, we guaranteed that the boundary layer, of size estimated by 2R Re -1/2 , had four minimum grid points.

  Description of mesh settings (see Tab.

3 . 2 )Figure 3 . 5 :

 3235 Figure 3.5: Description of the meshing procedure and zoom on the mesh in the vicinity of the droplet (the actual simulation domain is 50 times larger than the droplet radius).

  µ e e e Rθ = µ i e i Rθ (3.24) U e ψ = U i ψ , µ e e e θψ = µ i e i θψ (3.25)

Figure 3 . 6 :

 36 Figure 3.6: Evolution of the drag coefficient as a function of the Reynolds number. Comparison of simulation results (point) and correlations, Hadamard [1] and Rybczynski et al. [2] with Mei et al. [3] and Schiller and Naumann [4] (line) color corresponds to µ * : µ * = 20, µ * = 10, µ * = 5, µ * = 3, µ * = 2, µ * = 0.5, µ * = 0.33, µ * = 0.2, µ * = 0.1, µ * = 0.05.

80CHAPTER 3 .Figure 3 . 7 :

 337 Figure 3.7: Trajectory (top figures) and time evolution of the aspect ratio in experiments (bottom figures) of Heptane or n-Pentane droplets in water.

Figure 3 . 8 :

 38 Figure 3.8: Time evolution of the velocity for experimental cases, velocity modulus raw data , Y-velocity raw data, X-Velocity raw data, and rolling mean data.

Despite a small area

  of overlap in Reynolds between conditions giving rise to linear and oscillating trajectories (around 8% of the threshold value), we can identify the critical Reynolds numbers for the flow bifurcation of n-Pentane µ * = 0.231, Heptane µ * = 0.416 and Cyclohexane µ * = 1.051 droplets, respectively: Re = 249, Re = 238 and Re = 256. For comparison, the corresponding internal Reynolds, Eq. 3.2, are respectively equal to Re i = 668, Re i = 391 and Re i = 189. It's interesting to note that for the considered ascending droplets, oscillations are observed from external Re close to the critical bifurcation Re of a fixed solid sphere, equal to Re = 270.

82CHAPTER 3 .Figure 3 . 9 :

 339 Figure 3.9: Flow map based on experimental observations of droplets rising in water, straight vertical trajectory, oscillating trajectory. From left to right: n-Pentane µ * = 0.231, Heptane µ * = 0.416 and Cyclohexane µ * = 1.051.

  Figure 3.10: comparison of drag coefficient Reynolds number evolution of experimental cases , with axisymmetric simulation data .

Figure 3 . 11 :

 311 Figure 3.11: Experimental mean drag coefficient C d and comparison with corresponding axisymmetric simulations . Experimental results are reported with an error bar estimated from Eq. (3.20).

  the Reynolds number remains low, a Hill vortex takes place inside the droplet. Examples are shown in Fig.3.12, for 2 cases: µ * = 0.1, ρ * = 1 and Re = 20. (left column), and µ * = 10, ρ * = 0.2 and Re = 200 (right column). In the second case, where Re is higher and µ * > 2, an additional vortex develops at the rear of the droplet Rachih[START_REF] Rachih | Étude numérique du transfert de matière à travers l'interface d'une goutte sphérique en mouvement : mise en évidence des effets 3D[END_REF].In each case, since the flow is x-axis axisymmetric, the figures representing the xy and xz planes are identical. Figures 3.12e and 3.12f, representing the yz plane, further illustrate the absence of azimuthal velocity, with straight streamlines going from the center towards the surface of the drop. The corresponding vortex trajectories are shown in Fig.3.13, indicating a perfectly round circular vortex core.

Case 1 : 100 •

 1100 µ * = 0.2, ρ * = 1 and Re = Case 2: µ * = 5, ρ * = 1 and Re = 350 • Case 3: µ * = 1.05, ρ * = 0.782 and Re = 400 The 3 types of stable 3D flows obtained are illustrated in Fig. 3.14. Compared to the previous cases, very different characteristics can be observed. The non-symmetry of the flow is particularly visible in the yz plane (see Figs. 3.14g and 3.14i) where 4 distinct regions with non-zero azimuthal velocity can be observed inside the droplet for cases 1 and 3.

  (a) xy plane. µ * = 0.1, ρ * = 1, Re = 20 (b) xy plane. µ * = 10, ρ * = 0.2, Re = 200 (c) xz plane. µ * = 0.1, ρ * = 1 Re = 20 (d) xz plane. µ * = 10, ρ * = 0.2, Re = 200 (e) yz plane. µ * = 0.1, ρ * = 1, Re = 20 (f) yz plane. µ * = 10, ρ * = 0.2 and Re = 200

Figure 3 . 12 :

 312 Figure 3.12: Examples of axisymmetric flows. Left: low-viscosity droplet with low Re; right: viscous droplet with moderate Re, showing detachment of the outer boundary layer. Top: xy plane at z = 0; Middle: xz plane at y = 0; Bottom: yz plane at x = 0.

Figure 3 . 13 :Figure 3 . 14 :

 313314 Figure 3.13: Vortex core inside and outside the droplet, for the cases shown in Fig. 3.13

( a )

 a Vortex cores, for µ * = 0.2, ρ * = 1 and Re = 100 (b) Longitudinal vorticity isovalue, blue and red have opposite signs, for µ * = 5, ρ * = 1 and Re = 350 (c) Vortex cores, for µ * = 1.05, ρ * = 0.782 and Re = 400

Figure 3 . 15 :

 315 Figure 3.15: Comparison of internal and external bifurcation consequences for the 3D steady flows of Fig. 3.14.

88CHAPTER 3 .

 3 3D BIFURCATIONS STUDY (a) xy slice at z = 0, for µ * = 0.05, ρ * = 1 and Re = 50 (b) xy slice at z = 0, for µ * = 10, ρ * = 0.2 and Re = 500 (c) xz slice at y = 0, for µ * = 0.05, ρ * = 1 and Re = 50 (d) xz slice at y = 0, for µ * = 10, ρ * = 0.2 and Re = 500 (e) yz slice at x = 0, for µ * = 0.05, ρ * = 1 and Re = 50 (f) yz slice at x = 2, for µ * = 10, ρ * = 0.2 and Re = 500

Figure 3 . 16 :

 316 Figure 3.16: Examples of unsteady 3D flow cases.

Figure 3 . 17 :

 317 Figure 3.17: Comparison of unsteady 3D flow cases.

90CHAPTER 3 .Figure 3 . 18 :

 3318 Figure 3.18: Determination of the bifurcation origin from single phase flow simulations and visualization for µ * = 5, ρ * = 1 and Re = 350

Figure 3 . 19 :

 319 Figure 3.19: Determination of the bifurcation origin from single phase flow simulations and visualization for µ * = 02, ρ * = 1 and Re = 200

Figure 3 . 20 :

 320 Figure 3.20: Two identical meshes and simulation settings (µ * = 0.1, ρ * = 1, and Re = 50), performed with two initial random perturbations. The streamlines only inside the droplet are presented to highlight the orientation of symmetry planes. The flow direction follows the x-axis which is orthogonal to the plane of visualization.

( a ) 1 Figure 3 . 21 :

 a1321 Figure 3.21: Evolution of maximum azimuthal velocity w max and damping or amplification factor α as a function of time. The grey region represents the plateau where the damping or amplification factor is determined.

  3.22b, and 3.22c. 94CHAPTER 3. 3D BIFURCATIONS STUDY (a) Effect of Reynolds number on amplification factor at: ρ * = 0.2, µ * = 0.05 (b) Effect of density ratio on amplification factor at: Re = 50, µ * = 0.05 (c) Effect of density ratio on amplification factor at: Re = 500, µ * = 20

Figure 3 . 22 :

 322 Figure 3.22: Evolution of amplification factor in function of Reynolds number or density ratio.

Figure 3 . 23 :

 323 Figure 3.23: Flow phase diagram for ρ * = 1, axisymmetric cases, steady internal bifurcation marginal state, steady 3D flow state with internal bifurcation, 3D unsteady internal bifurcation, unsteady internal bifurcation marginal state, steady external bifurcation marginal state, steady 3D flow state with external bifurcation, 3D unsteady flow with external bifurcation, unsteady external bifurcation marginal state.

Figure 3 . 24 :

 324 Figure 3.24: Phase diagram of external bifurcation with and without the effect of density ratio, axisymmetric cases, ρ * ̸ = 1 axisymmetric cases, external bifurcation marginal state, external bifurcation marginal curve, external steady bifurcation cases, external bifurcation unsteady cases, external bifurcation unsteady marginal curve, ρ * ̸ = 1 external bifurcation cases.

98CHAPTER 3 .Figure 3 . 25 : 3 ( 3 . 34 )

 33253334 Figure 3.25: Internal bifurcation map at ρ * = 1, axisymmetric cases, internal steady bifurcation critical point, internal steady bifurcation critical line, steady state internal bifurcation cases, internal unsteady bifurcation cases, internal unsteady critical bifurcation cases.

3. 6 . 15 ( 3 . 37 )

 615337 DNS STUDY OF THE TRANSITION FROM 2D TO 3D FLOW 99 Re c,unsteady (µ * ) = 2.1e [0.832 ln(µ * )+5.59] e Re i,c,unsteady (µ * ) = 2.1e [0.832 ln(µ * )+5.59] e

  3.26a. If one still wants to use the external Reynolds number, it is sufficient to divide the x-axis by ρ * to compare the data, see Fig.3.26b. We have a correct agreement with previous results at ρ * = 1, the density ratio has only a small impact independently of Re i (like for the drag coefficient), and therefore we can model the bifurcation criterion for the internal bifurcation with the impact of the density ratio, see Eq. (3.39) and equivalently expressed in terms of external Reynolds number(3.40).

( a )Figure 3 . 26 : 39 )

 a32639 Figure 3.26: Internal bifurcation map with density ratio parametric study simulations, axisymmetric cases, ρ * ̸ = 1 axisymmetric cases, internal steady bifurcation critical point, internal steady bifurcation critical line, internal steady bifurcation cases, internal unsteady bifurcation cases, internal unsteady bifurcation critical line, ρ * ̸ = 1 internal bifurcation cases.

Figure 3 . 27 :

 327 Figure 3.27: Time evolution of drag coefficient in function of time for simulations presented in Sec. 3.6.1.

Figure 3 . 28 :

 328 Figure 3.28: 3D drag coefficient divided by the axisymmetric drag coefficient in function of the Reynolds number, simulation results and linear fit .

Figure 3 . 29 : 5 :

 3295 Figure 3.29: 3D drag coefficient divided by the axisymmetric drag coefficient in function of the viscosity ratio, simulation results and regression

(a) µ * = 10 ,Figure 3 . 30 :

 10330 Figure 3.30: Comparison of density ratio effect on drag coefficient for internal and external bifurcations.

  Cd 3DCd axi (ρ * = 1) = 0.120. So for internal bifurcation, we could estimate the 3D drag coefficient at different ρ * by computing the corresponding axisymmetric drag coefficient and use Eq.(3.41) with the corresponding internal Reynolds number. 104CHAPTER 3. 3D BIFURCATIONS STUDY The evolution of the drag coefficient in Fig. 3.27 reveals that the 3D unsteady flow results in an oscillation of the drag coefficient. From this plot, the oscillation of larger amplitude is characterized by the Strouhal number St. Fig. 3.31 reports the evolution of the Strouhal number as a function of Reynolds number for different viscosity ratios. The density ratio here is fixed to ρ * = 1. The external Reynolds number is used for the external bifurcation (µ * > 2) while the internal Reynolds number is used for the internal bifurcation (µ * < 2). For internal bifurcation (left figure column) the Strouhal increases with the internal Reynolds number until a critical point located at Re i ≈ 800, where it starts to flatten or decrease. For external bifurcation (right figure column), the Strouhal increases with the internal Reynolds number also, but the evolution is not as clear, and we cannot establish a critical point where it begins to flatten or decrease.

  (a) µ * = 0.1 (b) µ * = 5 (c) µ * = 0.2 (d) µ * = 20

Figure 3 . 31 :

 331 Figure 3.31: Strouhal number evolution with Reynolds number for simulation results.

3. 7 . 105 3. 7

 71057 DISCUSSION OF EXPERIMENTAL AND NUMERICAL OUTCOMES. Discussion of experimental and numerical outcomes.

Figure 3 . 32 :

 332 Figure 3.32: Bifurcation map from experimental data, straight trajectory, oscillating trajectory, internal steady bifurcation critical line, and internal steady bifurcation critical line.

Figure 3 . 33 :

 333 Figure 3.33: Internal bifurcation map with experimental data, axisymmetric cases, internal steady bifurcation critical line, steady state internal bifurcation cases, internal unsteady bifurcation cases, internal unsteady critical bifurcation cases.

  the drag coefficient in both experiments and simulations, and we compare experimental and numerical results in a figure where the drag coefficient is plotted as a function of Reynolds number, Fig.3.34. To make the comparison, specific simulations were carried out with the specific viscosity and density of the fluids used in the experiments. The range of Reynolds number observed in the experiments is discretized with a step of 10.

Figure 3 . 34 :

 334 Figure 3.34: Drag coefficient from the simulations ( (axisymetric) and 3D steady and for unsteady) and experimental results with error bar.

108CHAPTER 3 .

 3 3D BIFURCATIONS STUDY flow. So we expect the simulation and experimental Strouhal numbers to be different. Experimental strouhal data are presented in Fig. 3.35, we add simulation Strouhal data are presented in Fig. 3.36.The Strouhal number St is used to characterize unsteady wakes, as observed here for n-Pentane at Re = 364. From Ern et al.[START_REF] Ern | Wakeinduced oscillatory paths of bodies freely rising or falling in fluids[END_REF] and Kiya et al.[START_REF] Kiya | Near-wake instabilities and vortex structures of three-dimensional bluff bodies: A review[END_REF], we know that for fixed solid sphere simulations, the Hopf bifurcation yields Strouhal St ≈ 0.127. From Ern et al.[START_REF] Ern | Wakeinduced oscillatory paths of bodies freely rising or falling in fluids[END_REF], it appears that a free-rising sphere instead of a fixed sphere reduces the Strouhal number and changes the threshold Reynolds number. But they were able to explain the discrepancy between Fig.3.35a and 3.35c. In our experimental study, unsteady bifurcation (Fig.3.7c) corresponds to Strouhal number St = 0.458 which indicates that the vortices are moderately frequent and strong. An important consideration would be to observe the evolution of Strouhal number with the Reynolds number to compare with the solid particle case, see Ern et al.[START_REF] Ern | Wakeinduced oscillatory paths of bodies freely rising or falling in fluids[END_REF]. All our experimental post-processed cases are gathered in Fig.3.35c.

  From Kiya et al.[START_REF] Kiya | Near-wake instabilities and vortex structures of three-dimensional bluff bodies: A review[END_REF] article (solid sphere) , and our experimental data .

Figure 3 . 35 :

 335 Figure 3.35: Evolution of Strouhal number with Reynolds number, experimental results for a spherical droplet freely rising in water.

  Data from Kiya et al.[START_REF] Kiya | Near-wake instabilities and vortex structures of three-dimensional bluff bodies: A review[END_REF] , simulation , and experimental results (both N-pentane and Heptane).

Figure 3 . 36 :

 336 Figure 3.36: Strouhal number comparison of simulations , and experimental results .

  sp(a, b, e) = (len(a) -1)σ(a) 2 + (len(b) -1)σ(b) 2 + len(a)len(b)e 2 len(a) + len(b) -2 , (3.43) with a and b the two groups, σ the standard deviation and e the uncertainty, see Eq. (3.20). The t-statistic is computed as the difference between the means of the two groups divided by the pooled standard error. It quantifies the standardized difference between the sample means, see Eq. (3.44). The p-value is calculated using the t-distribution and the cumulative distribution function (CDF) and represents the probability of obtaining the observed t-statistic under the assumption that there is no true difference between the group means, see Eq. (3.45). t-stat(a, b, sp) = a -(a, b, t-stat) = 2(1 -CDF (|t-stat|, df = len(a) + len(b) -2)). (3.45)

114CHAPTER 3 2 4 . 3 4 . 4

 324344 Simulation strategy . . . . . . . . . . . . . . . . . . . . . . . 119 4.2.1 Governing equations and mesh . . . . . . . . . . . . . . . . . 119 4.2.2 Definition of dimensionless numbers . . . . . . . . . . . . . . 122 4.2.3 Methodology . . . . . . . . . . . . . . . . . . . . . . . . . . . 124 4.2.4 Parametric study . . . . . . . . . . . . . . . . . . . . . . . . . 125 Analysis of mass-transfer regimes . . . . . . . . . . . . . . . 126 4.3.1 Internal problem and problems with expected high internal resistance, k √ D * < 1 . . . . . . . . . . . . . . . . . . . . . . 126 4.3.2 External problem and problems with expected high external resistance, k √ D * > 1 . . . . . . . . . . . . . . . . . . . . . . 133 New correlations for Sh i and Sh e valid in conjugated problems140 4.4.1 Conjugated mass-transfer problem: need for new correlations 140 4.4.2 Correlation of the internal mass transfer coefficient for coupled problems, Sh i . . . . . . . . . . . . . . . . . . . . . . . . . . . 143 4.4.3 Correlation of the external mass transfer coefficient for coupled problems, Sh e . . . . . . . . . . . . . . . . . . . . . . . . . . . 146 4.4.4 Relevance of the new correlations . . . . . . . . . . . . . . . . 149 4.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 151 116CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER

  interface to that transported by diffusion (Eq. (4.2)): d a typical dimension of the interface (by convention we will use the diameter of the droplet), and D i the diffusivity of the internal phase. The use of D i in the definition of the global Sherwood results from the use of internal and external sherwood number expression, Sh i = h i d D i and Sh e = hed De (4.3), with h i and h e the internal and external mass transfer coefficients. The external Sherwood number definition simplifies with the diffusivity ratio to the internal diffusivity.

118CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER (a) Isolated internal problem (b) Isolated external problem (c) Coupled or conjugated problem

Figure 4 . 1 :

 41 Figure 4.1: Illustration of possible mass-transfer regimes in the film theory.

Figure 4 . 2 :

 42 Figure 4.2: Detail of the domain discretization: The droplet, with a polar mesh, is in the center. The green and red lines, extending respectively on -1 < x/R < 2 and 0 < y/R < 2 in the parallel and perpendicular directions of the external flow, materialize the lines used to plot the concentration profiles in Sec. 4.3.

Figure 4 . 3 :

 43 Figure 4.3: Possible configurations encountered in solvent extraction processes. The arrows indicate the direction of the transfer: red for extraction, green for stripping. Blue (resp. yellow) stands for the aqueous phase (resp. organic): (a) and (c) represent O/W configurations (for which µ * > 1), (b) and (d) W/O configurations (µ * < 1).

124CHAPTER 4 .Figure 4 . 4 :

 444 Figure 4.4: Time evolution of the global (•), internal ( * ), and external (△) Sherwood numbers in the case of a solute transferred from the continuous phase to the droplet (a) (initial concentrations C i = 0, C e = 1, black lines) and from the droplet to the continuous phase (c) (initial concentrations C i = 1, C e = 0, red symbols). In both cases: Re = 20, µ * = 5, k = 10, P e e = 100.

FirstFig. 4 .

 4 Fig. 4.5 is typical of a viscous droplet (µ * = 5, O/W) from which transfer is favoured by thermodynamics (k = 0.1, stripping). The corresponding concentration profiles are shown in Fig. 4.6. The time evolution of the transfer process can be assessed from Fig. 7.4. According to the small value of k √ D * = 0.04, this configuration is likely to be governed by the internal mass-transfer resistance (Sh → Sh I ).

128CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER (a) Re = 20, P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 4 . 5 :

 45 Figure 4.5: Streamlines and concentration contours in the case "Oil in Water" (µ * = 5). Effect of the Reynolds and external Péclet numbers for k = 0.1 (k √ D * = 0.04).

4. 3 .

 3 ANALYSIS OF MASS-TRANSFER REGIMES 129 and red profiles are no longer superimposed. The position of the maximum concentration point has a large impact on the concentration gradient at the interface, and a more pronounced decrease in the external concentration is observed in the perpendicular direction. (a) P e = 10 (b) P e = 1000

Figure 4 . 6 :

 46 Figure 4.6: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.5. "Oil-in-Water", k √ D * = 0.04 (Symbols: Re = 20, lines Re = 200).

22 Fig. 4 .

 224 Fig.4.7 is typical of a W/O configuration (µ * = 0.2) for which, as in Sec. 4.3.1, the transfer of the solute is favored by the thermodynamic equilibrium that prevails at the interface, with the same constant k = 0.1. The two configurations hence differ only by their viscosity (resp. diffusivity) ratio, which is smaller (resp. higher) in the present case,

Figure 4 . 7 :

 47 Figure 4.7: Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the Reynolds and external Péclet numbers for k = 0.1 (k √ D * = 0.22).

  7.5).

  (a) P e = 10 (b) P e = 1000

Figure 4 . 8 :

 48 Figure 4.8: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.7. "Water-in-Oil", k √ D * = 0.22. (Symbols: Re = 20, Lines: Re = 200).

7 . 1 )

 71 Additional case-study at k√ D * = 0.22 with k = 1Although k √ D * is the same as in the previous case, the O/W configuration of Fig.4.9 (µ * = 20, k = 1) is similar to a heat transfer problem where no concentration (temperature) discontinuity prevails at the interface and only transport mechanisms (i.e.

132CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER convection and diffusion) can affect the rate of mass-transfer. The corresponding concentration profiles are shown in Fig. 4.10. The steady-state Sherwood numbers are gathered in Tab. 7.1, and their time evolutions are plotted in Fig. 7.6.

  (a) Re = 20, P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 4 . 9 :

 49 Figure 4.9: Streamlines and concentration contours in the case "Oil in Water" (µ * = 20). Effect of the Reynolds and external Péclet numbers for k = 1 (k √ D * = 0.22).

  (a) P e = 10 (b) P e = 1000

Figure 4 . 10 :

 410 Figure 4.10: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.9. "Oil-in-Water", k √ D * = 0.22 with k = 1. (Symbols: Re = 20, Lines: Re = 200).

Fig. 4 .

 4 Fig.4.11 is typical of a W/O configuration (µ * = 0.2). Given the high value of the Wylock et al.[START_REF] Wylock | Gas adsorption into a spherical liquid droplet: numerical and theoretical study[END_REF] criterion, transfer is expected to be governed by external transfer. This is indeed consistent with the uniform distribution of concentration inside the droplet and the rather constant concentration profiles shown in Fig.4.11 and 4.12.

(a) Re = 20 ,

 20 P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 4 . 11 :

 411 Figure 4.11: Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the Reynolds and external Péclet numbers for k = 10 (k √ D * = 22).

136CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER (a) P e = 10 (b) P e = 1000

Figure 4 . 12 :

 412 Figure 4.12: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.11. "Water-in-Oil", k √ D * = 22.36. (Symbols: Re = 20, Lines: Re = 200).

: k √ D * = 4 . 47 We

 447 consider here an O/W configuration (µ * = 5), for which the concentration contours are shown in Fig. 4.13, the corresponding concentration profiles in Fig. 4.14, and the time evolution of the Sherwood numbers in Fig 7.8. The corresponding steady-state values are gathered in Tab. 7.1.

(a) Re = 20 ,

 20 P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 4 . 13 :

 413 Figure 4.13: Streamlines and concentration contours in the case "Oil in Water" (µ * = 5). Effect of the Reynolds and external Péclet numbers for k = 10 (k √ D * = 4.47).

  7.8).

138CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER (a) P e = 10 (b) P e = 1000

Figure 4 . 14 : 1 Fig. 4 .

 41414 Figure 4.14: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.13. "Water-in-Oil", k √ D * = 4. (Symbols: Re = 20, Lines: Re = 200).

(a) Re = 20 ,

 20 P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 4 . 15 :

 415 Figure 4.15: Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.05). Effect of the Reynolds and external Péclet numbers for k = 1 (k √ D * = 4.47).

  (a) P e = 10 (b) P e = 1000

Figure 4 . 16 : 4 . 4 4 . 4 . 1

 41644441 Figure 4.16: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.15. "Water-in-Oil", k √ D * = 4.47 with k = 1 (Symbols: Re = 20, Lines: Re = 200).

142CHAPTER 4

 4

Figure 4 . 17 :

 417 Figure 4.17: Parity plots of isolated problem correlation results compared to coupled problem simulations results.

4. 4 .

 4 NEW CORRELATIONS FOR SH I AND SH E VALID IN CONJUGATED PROBLEMS143

Figure 4 . 18 :

 418 Figure 4.18: Comparison between maximum interfacial velocity U max correlation fitted for droplet in function of µ * (line) and DNS results (points).

  ) and (4.18) respectively. Therefore, similarly to Colombet et al.[START_REF] Colombet | Mass or heat transfer inside a spherical gas bubble at low to moderate Reynolds number[END_REF] for Sh I , we propose here to define asymptotic values of the internal Sherwood numbers Sh i at P e → 0 and at P e → ∞, that account for the effect of µ * and k. The results are shown in Fig.4.19.

4. 4 .

 4 NEW CORRELATIONS FOR SH I AND SH E VALID IN CONJUGATED PROBLEMS145 Comparison between Eq. (4.17) (Newman, dash line) and Eq. (4.31) (red line). Symbols stand for DNS results.

  Comparison between Eq. (4.18) (Kronig, dash line) and Eq. (4.32) (blue: k = 0.1, red: k = 1, green: k = 10). Symbols stand for DNS results.

Figure 4 . 19 :

 419 Figure 4.19: Asymptotic values of the internal Sherwood number, Sh i , at low (left) and high (right) Péclet.

4. 4 .

 4 NEW CORRELATIONS FOR SH I AND SH E VALID IN CONJUGATED PROBLEMS147 (a) Case of bubbles (µ * = 0) (b) Case of droplets (µ * = 2) (c) Case of solid spheres (µ * → ∞).

Figure 4 . 20 :

 420 Figure 4.20: Comparison between Sh E Eq. (4.27) (line) and DNS results Sh e (points) for coupled problems at Re = 10. k = 0.1, k = 0.25, k = 1, k = 4 and k = 10.

148CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER

  4.21) lets us consider a correction for coupled problems, in the case of bubbles (µ * = 0 in Eq. (4.28)).

Figure 4 . 21 :

 421 Figure 4.21: Evolution of simulation Sh e with k for bubbles (µ * = 0) at Re = 10 and P e = 100 , and comparison with value from Eq. (4.28) for bubble, Sh E (black dashed line).

150CHAPTER 4 .

 4 CORRELATIONS FOR CONJUGATED MASS TRANSFER

  Global Sherwood from Eq.(4.3)

Figure 4 . 22 :

 422 Figure 4.22: Parity plots of new correlation results compared to coupled problem simulations results.

2 5 . 3 5 . 4

 25354 Problem statement and governing equations . . . . . . . . 5.2.1 Numerical model . . . . . . . . . . . . . . . . . . . . . . . . . 5.2.2 Mesh refinement validation . . . . . . . . . . . . . . . . . . . 5.2.3 Simulation strategy and explored domain . . . . . . . . . . . 5.2.4 The pure diffusion case . . . . . . . . . . . . . . . . . . . . . Results and discussion . . . . . . . . . . . . . . . . . . . . . 5.3.1 Enhancement effect of a chemical reaction with expected high external resistance according to the film theory ( √ kD * > 1) . 5.3.2 Enhancement effect of the chemical reaction with expected high internal resistance according to the film theory ( √ kD * < 1)168 5.3.3 Summary . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Predicting the global Sherwood number for conjugate transfer enhanced by chemical reaction . . . . . . . . . . . 5.4.1 Non-reactive mass transfer coefficients . . . . . . . . . . . . . 5.4.2 Enhancement factor . . . . . . . . . . . . . . . . . . . . . . . 5.4.3 Accuracy of the model prediction . . . . . . . . . . . . . . . . 5.5 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 154CHAPTER 5.

Figure 5 . 1 :

 51 Figure 5.1: Film theory scheme for conjugate mass transfer.

158CHAPTER 5 .

 5 SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION

e= 0 .

 0 Fig. 5.2 gives an illustration of flow and concentration at the initial time in a general case. These DNS results highlight the internal flow occurring within the moving droplet, that is likely to influence the concentration patterns in the droplet, which should deviates from the diffusion profiles expected in solid particles.

Figure 5 . 2 :

 52 Figure 5.2: Simulation configuration at initial time with streamlines.

. 24 )

 24 ∂C e ∂t + ∇ • (U e C e ) = 1 P e e ∇ 2 C e -DaC e (5.25) 160CHAPTER 5. SOLUTE EXTRACTION FROM A SPHERICAL DROPLET ENHANCED BY AN EXTERNAL FIRST ORDER CHEMICAL REACTION where U δ = U δ U∞ and C δ = C δ C ini i are the dimensionless velocity and concentration, while the droplet diameter d = 2R is chosen as the characteristic length scale. The dimensionless momentum equations for the internal and external phases show the Reynold numbers, defined as Re i = 2ρ i U∞R µ i and Re e = 2ρeU∞R µe . Similarly, the dimensionless solute transport equations show the internal and external Peclet numbers, respectively defined as P e i = 2U∞R D i and P e e = 2U∞R De , and finally the Damköhler number Da = νR 2 De .

  (5.21)-(5.25), and the corresponding interface boundary conditions (dimensionless form of Eqs. (5.18)-(5.19

Figure 5 . 3 :

 53 Figure 5.3: Simulation domain and mesh (left) and validation test (right). (a) Domain discretization, zoom on the droplet. The red and green lines are used in Sec. 5.3 to plot the axial and perpendicular concentration profiles, respectively. (b) Comparison between the enhancement factor E deduced from DNS (symbols) and the analytical solution Eq. (7.12) (line) for different Da values.

162CHAPTER 5 .

 5 Figure 5.4: Simulation strategy for the mass transfer investigation based on avalaible hydrodynamic simulations from Godé et al. [5]

5. 3 . 1

 31 Enhancement effect of a chemical reaction with expected high external resistance according to the film theory ( √ kD * > 1) This case with µ * > 1 is typical of an oil droplet in water. Fig. 5.5 compares the time-evolutions of the Sh, Sh i and Sh e obtained by DNS without reaction (Da = 0), and for increasing kinetics of the chemical reaction (1 ≤ Da ≤ 100). For this coupled problem, unfavorable to transfer due to both a moderate advection effect with P e = 10 and the thermodynamic equilibrium condition with k = 10, mass transfer is observed to be globally enhanced by the chemical reaction: the steady-state Sh value being higher when Da increases (Fig. 5.5c) with steady values of 3.59, 4.67, 10.2 and 29.3 for Da = 0, 1, 10 and 100, respectively.

Figure 5 . 5 :

 55 Figure 5.5: Time evolution of the Sherwood numbers obtained by DNS for Re = 20, D * = 0.2, k = 10 and P e e = 10 (with F o = Det R 2 the dimensionless time). Effect of an external first order chemical reaction with Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). Steady external Sherwood values are reported in Tab. 7.2.

Figure 5 . 6 :

 56 Figure 5.6: Time evolution of the mean concentration in the droplet (line) at the droplet interface (point): effect of an external first order chemical reaction with Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue) (Re = 20, D * = 0.2, k = 10 and P e e = 10). Concentrations are normalized by the mean concentration in the droplet.

100 Figure 5 . 7 :

 10057 Figure 5.7: Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 0.2, k = 10 and P e e = 10

166CHAPTER 5 .Figure 5 . 8 :

 558 Figure 5.8: Concentration profiles for Re = 20, D * = 0.2, k = 10 and P e e = 10 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet.

( a ) 100 Figure 5 . 9 :

 a10059 Figure 5.9: Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 0.2, k = 10 and P e e = 1000

168CHAPTER 5 .Figure 5 . 10 :

 5510 Figure 5.10: Concentration profiles for Re = 20, D * = 0.2, k = 10 and P e e = 1000 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet.

( a ) 100 Figure 5 . 11 :

 a100511 Figure 5.11: Concentration spatial distribution and streamlines for different Da at Re = 20, D * = 5, k = 0.1 and P e e = 10

Figure 5 . 12 :

 512 Figure 5.12: Concentration profiles for Re = 20, D * = 5, k = 0.1 and P e e = 10 at Da = 0 (green), Da = 1 (lightblue), Da = 10 (blue), Da = 100 (navyblue). (a) Along the flow direction (green line in Fig. 5.3a), (b) Perpendicular to flow direction (red line in Fig. 5.3a). Concentrations are normalized by the mean concentration in the droplet.

= 100 Figure 5 . 13 :

 100513 Figure 5.13: Concentration spatial distribution and streamlines for different Da at Re = 200, D * = 5, k = 0.1 and P e e = 1000.

174CHAPTER 5 .Figure 5 . 14 :

 5514 Figure 5.14: Enhancement factor: comparison of simulation results (symbols, colored from light blue for Da = 1 to dark blue for Da = 100) with predictions of the simplified models (solid line: E F T , dash line: E SF T ). (a) Hatta from film theory, (b) General expression of Hatta from Juncu [6].

Figure 5 . 15 :

 515 Figure 5.15: Enhancement factor: comparison of simulation results (symbols, colored from light blue for P e e = 10 to dark blue for P e e = 1000) with predictions of the simplified models (continuous line: E F T , dashed line: E SF T . (a) Hatta from classic film theory, (b) General expression of Hatta from Juncu [6].

  .16d and colored by their Sh e(Da=0) values. Figs 5.16a, 5.16b and 5.16c illustrate for comparison the parity plots for internal and external Sherwood numbers with correlations from respectively Eq. (5.31) and (5.35) with Eq. (5.36) (i.e combining Sh e with enhancement factor based on Ha J ).

176CHAPTER 5 .Figure 5 . 16 :

 5516 Figure 5.16: Parity plots of external, internal and global Sherwood numbers obtained from DNS simulation and model prediction. In each case, the darker the symbol, the greater the value of the parameter, i.e. D * in (a), k in (b), Sh i in (c) and Sh e in (d).The straight lines reprensent the limits of 0% error (black), 10% error (blue) and 30% error (red).

6 . 1 6 . 2 6 . 3

 616263 On droplet hydrodynamics . . . . . . . . . . . . . . . . . . . On coupled mass transfer from a droplet . . . . . . . . . . Towards typical conditions of industrial implementation . 180CHAPTER 6. CONCLUSION AND OUTLOOK

6. 3 .

 3 Figure 6.1: Evolution of Sh e , Sh i , and Sh with P e e at steady-state, in the case of an internal bifurcation (Re = 50, µ * = 0.1, and k = 0.1), and comparison with the corresponding axisymmetric predictions (dashed line)

184CHAPTER 6 .

 6 CONCLUSION AND OUTLOOK(a) P e e = 10 (b) P e e = 1000

Figure 6 . 2 :

 62 Figure 6.2: Streamlines and concentration field visualization in the case of an internal bifurcation at Re = 50, µ * = 0.1, k = 0.1, at small (left) and elevated (right) external Péclet.

. 1 )a 1

 11 the numerical approximation of f consists in considering the following simplification of the complementary error function for x > 1 asf = (a 1 b + a 2 b 2 + a 3 b 3 + a 4 b 4 + a 5 b 5 ) × e y-x2 , = 0.254829592; a 2 = -0.284496736; a 3 = 1.42141741; a 4 = -1.453152027; a 5 = 1.061405429; p = 0.3275911. . Considering now the expression of the memory kernel in equation (2.5), we see that the corresponding x and y of equation (7.1) have a direct relation: y = x 2 , that simplifies even more the expression of the memory kernel because we get rid of every exponential terms as shown in equation (7.3) preventing overflow errors: f = (a 1 b + a 2 b 2 + a 3 b 3 + a 4 b 4 + a 5 b 5 ). (7.3)

  ′ a(a 1 b + a 2 b 2 + a 3 b 3 + a 4 b 4 + a 5 b 5 )dt ′ , (7.4) with: 7.1. APPROXIMATION OF THE KERNEL TO AVOID OVERFLOW ERROR 189

Figure 7 . 1 :

 71 Figure 7.1: Comparison of kernel computation with or without the approximation shown in Eq. 7.3. The curve without approximation is computed using the python library mpmath to ensure accuracy and avoid overflow errors.

7. 3 . 7 . 3

 373 SITUATIONS OF COUNTER-TRANSFER AND P ÉCLET NUMBER 191 Situations of counter-transfer and Péclet number (a) Contour of concentration at P e = 10 (b) Contour of concentration at P e = 100

Figure 7 . 2 :

 72 Figure 7.2: Streamlines and concentration contours in the case "Water in Oil" (µ * = 0.2). Effect of the external Reynolds and Péclet numbers for k = 0.1 (k √ D * = 0.22).

Figure 7 . 3 :CHAPTER 7 . APPENDIX 7 . 4 Figure 7 . 4 :

 7377474 Figure 7.3: Concentration profiles along (green) and perpendicular (red) to the flow direction extracted from the contours in Fig. 4.7a and 7.2 "Water-in-Oil", Re = 20, µ * = 0.2, k √ D * = 0.22 with k = 0.1 (Lines: P e = 10, Symbols: P e = 100.

(a) Re = 20 ,

 20 P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 7 . 5 :

 75 Figure 7.5: Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.1. "Water-in-oil", k √ D * = 0.22.

Figure 7 . 6 :

 76 Figure 7.6: Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.1. "Oil-in-Water", k √ D * = 0.22 with k = 1.

  (a) Re = 20, P e = 10 (b) Re = 20, P e = 1000 (c) Re = 200, P e = 10 (d) Re = 200, P e = 1000

Figure 7 . 7 :

 77 Figure 7.7: Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Water-in-Oil", k √ D * = 22.36.

Figure 7 . 8 :

 78 Figure 7.8: Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Oil-in-Water", k √ D * = 4.47.

Figure 7 . 9 :

 79 Figure 7.9: Time evolution of the internal (Sh i ), external (Sh e ) and apparent (Sh ) Sherwood numbers for the case discussed in Sec. 4.3.2. "Water-in-Oil", k √ D * = 4.47 with k = 1.

6 )C

 6 By solving those differential equations and considering the following boundary condition: C(∞) = 0 and C(r = R) = C s , we obtain the solution in respectively Eq. (7.7) and (7.7) for non-reactive and reactive transfer. Ha -sinh r R Ha cosh (Ha) -sinh (Ha)   . (7.8) Or in normalized form (C * = C Cs and r * = r R ): * = exp Ha exp -r * Ha r * . (7.10)

E = 1 +Figure 7 . 10 :

 1710 Figure 7.10: Time evolution of the mean internal concentration in the droplet (line) and interfacial (point), for Re = 20, µ * = 5, k = 10 and P e e = 1000 at Da : Da = 0, Da = 1, Da = 10, Da = 100

Figure 7 . 12 : 100 7. 7 Figure 7 . 11 :

 7121007711 Figure 7.12: Time evolution of the mean internal concentration for Re = 200, µ * = 0.2, k = 0.1 and P e e = 1000 at Da : Da = 0, Da = 1, Da = 10, Da = 100

  1.1 Main correlations of the drag coefficient for a moving solid particle, droplet or bubble in a fluid. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . Error metrics from history and total force model compared to simulation results in a "turbulent" flow. . . . . . . . . . . . . . . . . . . . . . . . . 3.1 Measured fluid properties and corresponding dimensionless numbers of interest for the study. All values are calculated for droplet rising in water. 3.2 Definition of the parameters of the mesh. . . . . . . . . . . . . . . . . . 3.3 Comparison of the steady drag coefficients obtained with axisymmetric DNS simulations. Bold font: this study, Quotation marks: Rachih et al.

2.1 Range of dimensionless parameters (Re, µ * , ρ * , f * ) explored in numerical simulations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.2

  3.5 Drag coefficient evolution with Reynolds number linear interpolation data at ρ .6 t-stat and p-value for different hypothesis. . . . . . . . . . . . . . . . . 4.1 Evaluation of mesh dependence, comparison between coarse, normal and refined mesh. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4.2 Error metrics for isolated problem correlation results compared to coupled problem simulations results. . . . . . . . . . . . . . . . . . . . . . . . . 4.3 Error metrics for the proposed correlation results compared to coupled problem DNS results. . . . . . . . . . . . . . . . . . . . . . . . . . . . . and k on the internal Sh i , external Sh e and apparent Sh Sherwood numbers. The corresponding asymptotic values Sh I and Sh E given by Eq. 5.31 and 4.27 are given for comparison. . . . . . . . . 190 7.2 Effect of Re, P e, D * k and Da on the internal external sherwood Sh

	20	LIST OF TABLES
	À Mathilde,

* = 1 . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35.1 Maximum relative Difference, Mean Absolute Difference (MAD), Mean Absolute Percentage Difference (MAPD) and Determination coefficient R2 for the internal, external and total Sherwood numbers. . . . . . .

18 LIST OF TABLES 7.1 Effet of Re, P e, D * e and enhancement factor. . . . . . . . . . . . . . . . . . . . . . . . . . . 199 Merci à ma famille, qui dans l'écriture de ma thèse comme dans mes compétitions d'aviron, ont toujours été mes plus fidèles et meilleurs supporters. Même si j'ai échoué à vous faire comprendre les détails de mon sujet, votre soutien a été essentiel à ma réussite. J'aimerais également remercier ma belle famille pour leur aide et leur soutien. Enfin, c'est avec une certaine émotion que je tiens à remercier mon épouse Mathilde Mongeau. Durant ces trois années, j'ai eu l'occasion de partager avec elle de nombreux moments inoubliables, en particulier notre mariage. Mathilde m'a aidé à surmonter tous les moments difficiles de cette thèse et m'a permis de franchir tous les obstacles, qui, grâce à elle, semblaient être d'une facilité déconcertante. Résumé Les procédés d'extraction liquide-liquide sont largement exploités dans les industries chimiques et connaissent aujourd'hui un regain d'intérêt dans le secteur hydrométallurgique pour la mise en place des filières de recyclage. Ils sont également au coeur des étapes de traitement-recyclage des combustibles nucléaires irradiés. Face à ces nouveaux enjeux, la simulation numérique à un rôle majeur à jouer dans la mise au point ou l'adaptation des procédés. Dans ce cadre, nous devons disposer de modèles prédictifs, représentant l'impact de tous les phénomènes physiques concernés, aux différentes échelles. Cette thèse est consacrée à l'étude des couplages entre l'hydrodynamique et la cinétique de transfert de matière au niveau d'une goutte sphérique dans un écoulement uniforme d'une phase immiscible.

  extensive, corroborée par des expériences originales, nous avons identifié un critère de bifurcation et modélisé l'impact de la transition 2D/3D sur le coefficient de trainée stationnaire. Nous avons également mis en évidence l'effet favorable de cette bifurcation vers des écoulements 3D sur le transfert de matière, en promouvant le transport convectif du soluté dans et autour de la goutte.Ainsi, grâce à la simulation numérique directe, nous avons décrit un large panel d'écoulements uniformes ou oscillants, dans une plage de Reynolds (susceptible d'être rencontrée dans les procédés d'extraction) et de rapports de viscosité (bulles, gouttes, particules solides) et traduit leur impact sur le transfert de matière, le coefficient de trainée et la force d'histoire. Ceci a conduit à des corrélations directement applicables pour la simulation des procédés.

Ces forces sont déterminantes sur le temps de séjour des gouttes dans les contacteurs liquide-liquide, et pendant lequel les deux phases sont susceptibles d'échanger. Une revue bibliographique nous a permis d'identifier deux termes encore mal décrits dans la littérature : la force d'histoire et la trainée stationnaire en écoulement 3D. L'étude de la force d'histoire est très complexe dans le cas des gouttes à cause du couplage entre les écoulements interne et externe. Nous avons proposé et validé un modèle qui tire bénéfice de solutions analytiques pour une condition à la limite de type glissement partiel et donc applicable aux bulles et aux particules solides pour prédire l'évolution de cette force d'histoire pour le cas des gouttes. Enfin, nous avons étudié les bifurcations 3D de l'écoulement dans et autour d'une goutte sphérique. Grâce, à

21

LIST OF TABLES une étude paramétrique

Table 1

 1

.1. 32CHAPTER 1. INTRODUCTION Figure 1.4: Hill's vortex visualization.

Author

Analytical solution/ Correlation Validity Stokes

[START_REF] Gabriel | On the Effect of the Internal Friction of Fluids on the Motion of Pendulums[END_REF]

Table 1 .

 1

1: Main correlations of the drag coefficient for a moving solid particle, droplet or bubble in a fluid.

Boussinesq history force acting on a droplet in an oscillatory flow Contents 2.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . 36 2.2 Problem statement . . . . . . . . . . . . . . . . . . . . . . . 37 2.3 Numerical method . . . . . . . . . . . . . . . . . . . . . . . . 40 2.4 Flow field description . . . . . . . . . . . . . . . . . . . . . . 42 2.5 Discussion on force contributions . . . . . . . . . . . . . . . 48 2.6 History force in the frequency domain . . . . . . . . . . . . 50 2.7 Kernel for history force in the time domain . . . . . . . . 52 2.8 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57 2.9 Outlook, accuracy of history force kernel in a turbulent flow 58

	36CHAPTER 2. BASSET-BOUSSINESQ HISTORY FORCE IN AN OSCILLATORY FLOW

Table 2 .

 2 1: Range of dimensionless parameters (Re, µ * , ρ * , f * ) explored in numerical simulations

	15)

Table 2 .

 2

	= 0.2	0.0106	0.946	0.00516	0.995
	µ * = 0.5	0.00480	0.995	0.00999	0.983
	µ * = 2	0.00846	0.981	0.00859	0.986
	µ * = 5	0.0112	0.942	0.0107	0.980

The table 2.2 presents error metrics, revealing that the model accurately predicts the simulation results.

History force MAE History force R2 Total force MAE Total force R2 µ * 2: Error metrics from history and total force model compared to simulation results in a "turbulent" flow.

Table 3 .

 3

	18)

1: Measured fluid properties and corresponding dimensionless numbers of interest for the study. All values are calculated for droplet rising in water.

Table 3 . 3

 33

	µ * Re	1	10	20	50	100	150	200
	0.05	17.4	2.43	1.43 0.697 0.395 0.279 0.218
		'18.0' '2.45' '1.45' '0.700' '0.390' '0.280' '0.230'
	0.1	17.8	2.51	1.48 0.727 0.413 0.293 0.228
	0.2	18.6	2.65	1.57 0.781 0.449 0.320 0.251
		'19.1' '2.72' '1.6'	'0.79'	'0.45'	'0.32'	'0.25'
	0.33	19.5	2.80	1.67 0.846 0.492 0.353 0.278
		[19.9] [2.87] [1.71] [0.89]
		'20.0' '2.87' '1.7'	'0.85'	'0.49'	'0.35'	'0.28'
	0.5	20.3	2.96	1.79 0.919 0.542 0.393 0.311
		(20.7) (3.03) (1.82) (0.939) (0.552)		(0.317)
		'20.8' '3.03' '1.81' '0.92'	'0.54'	'0.39'	'0.31'
	1	23.4	3.29	2.02	1.08	0.659 0.487 0.391
		(22.4) (3.34) (2.04) (1.10) (0.666)		(0.397)
		[22.4] [3.33] [2.05] [1.12]
		'22.4' '3.34' '2.04' '1.08'	'0.66'	'0.49'	'0.39'
	2	23.9	3.62	2.25	1.25	0.802 0.614 0.503
		(24.0) (3.66) (2.26) (1.25) (0.803)		(0.504)
		'24.1' '3.66' '2.26' '1.25'	'0.8'	'0.61'	'0.5'
	3	24.7	3.79	2.37	1.33	0.879 0.687 0.573
		[24.8] [3.80] [2.38] [1.36]
	5	25.6	3.96	2.48	1.42	0.954 0.762 0.648
		(25.7) (3.97) (2.48) (1.41) (0.955)		(0.646)
		'25.8' '3.99' '2.49' '1.42'	'0.95'	'0.76'	'0.64'
	10	27.1	4.12	2.59	1.49	1.02	0.824 0.709
		(26.4) (4.12) (2.58) (1.48)	(1.01)		(0.716)
		'26.5' '4.14' '2.6'	'1.49'	'1.02'	'0.82'	'0.71'

3.26) the correlations of Mei et al. [3], Eq. (3.27), and Schiller and Naumann : Comparison of the steady drag coefficients obtained with axisymmetric DNS simulations. Bold font: this study, Quotation marks: Rachih et al. [7], Parentheses: Feng and Michaelides

  Working with Eq. (3.29), we can obtain the Reynolds number as a function of the separation angle and the viscosity ratio, see Eq.(3.30). If we put a small separation angle (small not 0 because otherwise, the Reynolds number should always be equal to 0 according to the equation), and the different coefficients as variable, we obtain a base function that we could fit to our critical bifurcation points data. But to respect the analytical solution Re c (µ * → ∞) = 212, we know that the coefficient that multiply with the exponential should be equal to this solution, thus simplifying the optimization problem with only two variables, see Eq.(3.31).

	96CHAPTER 3. 3D BIFURCATIONS STUDY
	Re c (µ * ) = 20e	θ d	1+µ * 42.5µ *	1 0.483	(3.30)
	42.5 ln	Re 20		0.483	µ * 1 + µ *	(3.29)

Table 3 .

 3

		t-stat p-val
	n-pentane H 0 0.260 0.780
	heptane H 0	-0.300 0.765
	n-pentane H 1 -0.1659 0.916
	heptane H 1	-0.1756 0.923
	n-pentane H 2 -0.600 0.565
	heptane H 2	-0.336 0.739

• Hypothesis one (H 1 ): Experimental drag coefficient results at Re < Re c could be represented by corresponding axisymmetric simulations drag coefficient • Hypothesis two (H 2 ): Experimental drag coefficient results at Re > Re c could be represented by corresponding 3D simulations drag coefficient We only have enough data to have significant t-stat and p-val results for n-pentane and heptane experiments, values are reported in Tab. 3.6. 6: t-stat and p-value for different hypothesis.

Table 4 .

 4

		Coarse	Normal	Refined	extrapolated
		Sh	Error	Sh	Error	Sh	Error	Sh
	k = 10, P e e = 100	0.217 3.33% 0.212 0.95% 0.210	0%	0.0210
	k = 10, P e e = 10	0.0812 1.5% 0.0808	1%	0.0802 0.25%	0.0800
	k = 0.1, P e e = 1000	15.1	3.2%	15.4	0.91%	15.5	0.64%	15.6

1: Evaluation of mesh dependence, comparison between coarse, normal and refined mesh. 122CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER

  →∞,P e,Re) = 0.852P e 1/3 1 + 0.233Re 0.257 + 1.3 -0.182Re 0.355 , Sh E(µ * =2,P e,Re) = 0.64P e 0.43 1 + 0.233Re 0.257 + 1.41 -0.15Re 0.247 .

		1.032 +	0.61Re Re + 21	+ 1.6 -	0.61Re Re + 21	,
	Sh E(µ (4.28)
	First case-study: k	√ D * = 22.36

*

  4.4. NEW CORRELATIONS FOR SH I AND SH E VALID IN CONJUGATED PROBLEMS141

		Max error MAE MAPE	R 2
	Sh I Colombet et al. [83], Eq. (4.22)	10.5	2.98	0.22% 0.605
	Sh E Michaelides [35], Eq. (4.27)	9.01	2.09	0.91% 0.859
	Deduced Sh, Eq. (4.29)	20.6	7.39	2.18% -0.900
	Table 4.2: Error metrics for isolated problem correlation results compared to coupled
	problem simulations results.

  As reported in Sec. 4.3.1, the Sherwood number of the internal problem, Sh I , has two asymptotic values depending on the range of Péclet (P e → 0 and infinitely large Péclet).

	144CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER
	New asymptotic value of Sh i in coupled problems
	1/2 + 3Re 16 + 3.315Re 1/2 + Re	U ∞ 1 + exp 0.947(ln(µ * )-0.931) .	(4.30)
	Eq. (4.30) satisfactorily correlates our simulation results (R 2 = 0.958, MAE = 0.053
	and MAPE = 0.119 %).

  At high P e, the trend in Fig.4.19b indicates the existence of two asymptotic values for Sh i , one at high viscosity ratio, corresponding to Kronig (Sh i,P e→∞ = Sh Kronig = 17.9 for µ * → ∞), reached all the more quickly as k is large, and another one when µ * → 0, with a value around Sh i,P e→∞ = 25, which yields the following correlation: Sh i,P e→0 and Sh i,P e→∞ are expressed by Eq. (4.31) and (4.32) respectively. Eq. (4.33) is valid both for the internal problem (where Sh i → Sh I ) and for coupled problems (where Sh i ̸ = Sh I ).

	146CHAPTER 4. CORRELATIONS FOR CONJUGATED MASS TRANSFER
	where
			Sh i,P e→0 = 6.58 +	4.37 1 + (P e max ) 0.426 .	(4.31)
	Sh i,P e→∞ =	25 -17.9 1 + exp (0.684 ln(µ * ) -ln(0.496k))	+ 17.9 (1 + Re) 0.044 .	(4.32)
	From there, we deduced an expression of the internal Sherwood Sh i :
	Sh i =	Sh i,P e→∞ -Sh i,P e→0 0.762 + exp (-1.56 ln (P e max + 6.69))	+ Sh i,P e→0 ,	(4.33)

  Sh e(µ * =2,P e,Re) 0 ≤ µ * ≤ 2, →∞,P e,Re) = 1.000P e 1/3 1.32 + 0.177Re 0.308 -0.303Re 0.324 , Sh e(µ * =2,P e,Re) = 0.333P e 0.528 0.891 + 0.518Re 0.175 -0.110Re 0.346 . NEW CORRELATIONS FOR SH I AND SH E VALID IN CONJUGATED PROBLEMS149 As for Eq. (4.33), Eq. (4.34) is valid both for the external problem (where Sh e → Sh E ) and for coupled problems (where Sh e ̸ = Sh E ).

	4.4.
	µ * 6 + µ Sh e = 2 4µ * Sh e(µ * =0,P e,Re) + 1 µ * + 2 Sh e(µ * =2,P e,Re) + µ * + 2 µ * -2	Sh e(µ * →∞,P e,Re)	µ * > 2.	(4.34)
	With:
	Sh e(µ * =0,P e,Re) = 0.0158P e 1/2 21.8 + 10.6Re 0.170 + k 0.748 ,
	Sh e(µ (4.35)

* *

  The corresponding error metrics are reported in Tab. 4.3, again assessing the drastic improvement of the quality of the correlation, by comparison with Tab. 4.2. The proposed correlations can readily be used with the resistance additivity rule Eq. (4.3).

		Max error MAE MAPE	R 2
	Sh i , this work, Eq. (4.33)	4.98	0.814 0.065 % 0.971
	Sh e , this work, Eq. (4.34)	4.04	0.964 0.299% 0.974
	Deduced Sh, Eq. (4.3)	1.77	0.376 0.143% 0.994
	Table 4.3: Error metrics for the proposed correlation results compared to coupled
	problem DNS results.

  . It is corrected here by5.4. PREDICTING THE GLOBAL SHERWOOD NUMBER FOR CONJUGATE TRANSFER ENHANfitting on the new database provided by DNS in this study. Among other improvements, this relation also provides better account for the effect of the thermodynamic equilibrium constant k. Sh e(µ * =2,P e,Re) 0 ≤ µ * ≤ 2,

	e→∞ -Sh i,P e→0 0.762 + exp (-1.56 ln (P e max + 6.69)) where Sh i,P e→0 and Sh i,P e→∞ are expressed by Eq. (5.33) and (5.32) respectively: + Sh i,P e→0 , (5.31) Sh i,P e→∞ = 25 -17, 9 1 + exp (0.684 ln(µ * ) -ln(0.496k)) + 17, 9 (1 + Re) 0.044 . (5.32) and Sh i,P e→0 = 6.58 + 4.37 1 + (P e max ) 0.426 (5.33) with P e max a Péclet number defined using the maximum velocity at the droplet surface as: P e max = d D e 1 2 16 + 3.315Re 1/2 + 3Re 16 + 3.315Re 1/2 + Re U ∞ 1 + exp 0.947(ln(µ * )-0.931) = U max U ∞ P e. (5.34) 2 -µ * 2 Sh e(µ * =0,P e,Re) + 4µ * 6 + µ Sh e = 1 µ * + 2 Sh e(µ * =2,P e,Re) + µ * -2 µ * + 2 Sh e(µ * →∞,P e,Re) µ * > 2. (5.35) with Under this form, Eq. (5.31) is valid both for the internal problem (where Sh Sh e = Sh e(µ

i → Sh I ) and for coupled problems (where Sh i ̸ = Sh I ).

The correlation Eq.

(5.35)

for Sh e at Da = 0 is based on the work of Michaelides

[START_REF] Efstathios | Hydrodynamic force and heat/mass transfer from particles, bubbles, and drops -the freeman scholar lecture[END_REF]

who distinguished several regimes depending on the viscosity ratio µ * * * =0,P e,Re) = 0.0158P e 1/2 21.8 + 10.6Re 0.170 + k 0.784 , Sh e(µ * →∞,P e,Re) = 1.000P e 1/3 1.32 + 0.177Re 0.308 -0.303Re 0.324 , Sh e(µ * =2,P e,Re) = 0.333P e 0.528 0.891 + 0.518Re 0.175 -0.110Re 0.346 .

(5.36)

Table 5 .

 5 1: Maximum relative Difference, Mean Absolute Difference (MAD), Mean Absolute Percentage Difference (MAPD) and Determination coefficient R2 for the internal, external and total Sherwood numbers.

7.2 Steady-state Sherwood numbers from DNS pre- sented in Sec. 4.3.1 and Sec. 4.3.2

  Sh e Sh Sh I Sh E Sh(Sh I , Sh E )

	k Re Sh i 0.04 √ D * D * k P e 0.2 0.1 10 20 7.9	1.9	7.3 14.0	3.7	13.0
	0.04	0.2 0.1	10	200 9.5	2.0	8.7 15.7	4.1	14.6
	0.04	0.2 0.1 1000 20 21.4 18.9 20.9 18.1 17.7	17.8
	0.04	0.2 0.1 1000 200 20.5 23.6 20.2 18.3 22.7	18.0
	0.22	5.0 0.1	10	20	7.2	0.3	0.5	6.6	4.1	3.7
	0.22	5.0 0.1	10	200 6.6	0.4	0.7	6.7	4.4	3.8
	0.22	5.0 0.1 100	20	8.0	3.8	3.9	9.4 10.02	6.4
	0.22	5.0 0.1 1000 20 19.4 23.0 13.6 17.7 28.4	13.5
	0.22	5.0 0.1 1000 200 19.2 31.3 14.7 18.1 33.7	14.3
	0.22 0.05 1	10	20	8.4	3.3	7.4 17.7	3.6	14.2
	0.22 0.05 1	10	200 10.6 3.3	9.1 18.1	4.0	14.7
	0.22 0.05 1 1000 20 21.6 17.3 20.4 18.1 15.2	17.1
	0.22 0.05 1 1000 200 21.1 21.4 20.1 18.3 19.3	17.5
	4.47	0.2 10	10	20	9.9	3.5	1.5 13.2	3.6	1.6
	4.47	0.2 10	10	200 10.7 3.8	1.6 15.7	4.1	1.8
	4.47	0.2 10 1000 20 27.7 19.8 7.3 18.1 17.7	5.9
	4.47	0.2 10 1000 200 27.2 23.0 8.1 18.3 22.7	7.0
	4.47	20	1	10	20 10.4 1.6	0.1	6.6	4.0	0.2
	4.47	20	1	10	200 9.7	2.1	0.1	6.6	4.3	0.2
	4.47	20	1 1000 20 12.1 27.5 1.2 14.0 28.6	1.3
	4.47	20	1 1000 200 13.8 32.7 1.5 15.7 33.7	1.5
	22.36 5.0 10	10	20 10.2 3.9	0.1	6.6	4.1	0.1
	22.36 5.0 10	10	200 9.8	4.0	0.1	6.6	4.2	0.1
	22.36 5.0 10 1000 20 20.0 28.7 0.6 16.3 28.4	0.7
	22.36 5.0 10 1000 200 21.0 33.9 0.7 18.1 33.7	0.7

Table 7 .

 7 1: Effet of Re, P e, D * and k on the internal Sh i , external Sh e and apparent Sh Sherwood numbers. The corresponding asymptotic values Sh I and Sh E given by Eq. 5.31 and 4.27 are given for comparison.

7.6 Steady-state Sherwood numbers from DNS pre- sented in Chapter 5.

	Re µ *	k	P e Da Sh Juncu)
	20	5	10	10	1	4.67	3.59	1.30	1	1.31	0.557	1.10
	20	5	10	10	10	10.2	3.59	2.84	3.16	3.17	1.76	1.87
	20	5	10	10	100 29.3	3.59	8.16	10.0	10.0	5.57	5.57
	20	5	10 1000	1	19.4	19.3	1.01	1	1.31	0.103	1.01
	20	5	10 1000 10	19.8	19.3	1.03	3.16	3.17	0328	1.04
	20	5	10 1000 100 31.8	19.3	1.65	10.0	10.0	1.04	1.34
	200 0.2 0.1 1000	1	31.4	31.3	1.00	1	1.31	0.0638	1.00
	200 0.2 0.1 1000 10	32.2	31.3	1.01	3.16	3.17	0.202	1.01
	200 0.2 0.1 1000 100 39.5	31.3	1.04	10.0	10.0	0.639	1.13
	20 0.2 0.1	10	1	0.393	0.256	1.53	1	1.31	7.80	7.80
	20 0.2 0.1	10	10	1.95	0.256	7.62	3.16	3.17	24.7	24.7
	20 0.2 0.1	10	100 25.3	0.256	98.7	10.0	10.0	78.0	78.0

e Sh e (Da = 0) EF Hatte (film) EF (film) hatta (Juncu) EF (

Table 7 .

 7

2: Effect of Re, P e, D * k and Da on the internal external sherwood Sh e and enhancement factor.

Remerciements