Hassan El Itawi

Benjamin Lalanne

Gladys Massiera

Nathalie Le Sauze

Olivier Masbernat

Numerical simulation of the crossing of a liquid-liquid interface by a droplet

Numerical simulations of a drop crossing a plane liquid-liquid interface in a centrifugal field are performed using the level-set method. The objective is to characterize the hydrodynamical parameters controlling the coating volume of the droplet, which results from the rupture of the liquid column of the lighter phase entrained by the droplet during the crossing of the interface in the tailing regime. The numerical method is validated first in two-phase flow simulations of a drop rising in a stagnant liquid and then in three-phase flow configurations to reproduce the theoretical critical condition for a drop to cross an interface in static conditions (without initial velocity). Then, in inertial conditions, extensive simulations of crossing droplets are performed in a wide range of flow parameters and phase properties for two types of drop: solidlike droplets (mimicking rigid particles) and deformable drops. The crossing criterion is found to remain very close to that given by the theory in static conditions, for either a spherical or a deformed droplet. For each case studied, the crossing time, the maximum length of the column of liquid pulled by the droplet, and the volume encapsulating the drop after crossing are computed and scaled as a function of an inertia parameter, which is the ratio F * between the inertial stresses pushing on the interface, defined from the minimum drop velocity reached during crossing as the characteristic velocity, and the opposite stress pulling back the entrained column towards the interface. The maximal column length increases with F * (when rescaled by the minimal inertial required for crossing) under two distinct growth rates according to the flow regime in the column. For solidlike drops, the final coating volume is a unique function of F * and grows with F * at large inertia. In the case of deformable droplets, the coating volume evolution can also be well predicted by F * but corrected by the drop-to-film viscosity ratio, which strongly affects the drainage rate of the film along the drop surface during the encapsulation process.

I. INTRODUCTION

Encapsulation is a process by which an active ingredient is isolated from its external environment within a shell, in order to protect it and trigger and control its release in time and space. Several encapsulation processes already exist in the industry (pharmaceutical, food industry, cosmetics, etc.) [START_REF] Dubey | Microencapsulation technology and application[END_REF][START_REF]Handbook of Encapsulation and Controlled Release[END_REF]. The context of this work is the development of a microfluidic encapsulation process allowing the coating of submillimeter-sized droplets (the target diameter d being of order 100 μm) by a liquid phase immiscible with the droplet phase by crossing a liquid-liquid (typically oil-water) interface in a centrifugal field [START_REF] Massiera | [END_REF]. The same principle was applied at very low inertia for the production of giant vesicles [START_REF] Abkarian | Continuous droplet interface crossing encapsulation (cDICE) for high throughput monodisperse vesicle design[END_REF]. However, in the latter investigation, the crossing of the interface and the encapsulation process were driven by a different mechanism: the zipping of two amphiphilic monolayers adsorbed on both the droplet surface and the plane interface. In the present study, inertia forces, such as the FIG. 1. Scheme of the studied three-phase system in a centrifuged cell (ω is the angular speed): phase 1, lighter (organic) continuous phase; phase 2, droplet (aqueous); and phase 3, heavier (aqueous) continuous phase.

apparent weight or dynamic pressure force, are responsible for the crossing of droplets and result in the entrainment of a lighter liquid column into the heavier fluid (tailing regime) [START_REF] Maru | Behaviour of a rigid sphere at a liquid-liquid interface[END_REF][START_REF] Magnaudet | Particles, drops, and bubbles moving across sharp interfaces and stratified layers[END_REF], which, after rupture, leaves the droplet encapsulated.

The development of such a process first requires one to investigate the physics of the crossing of a liquid-liquid interface by a droplet in the inertial regime, including the determination of the crossing conditions and of the mechanism of encapsulation. In the process considered, a centrifugal field drives the droplet (phase 2) towards an interface separating two immiscible liquids, phases 1 and 3 (liquid 1 being the lighter phase), as sketched in Fig. 1. In order to ensure an inertial regime for submillimeter-sized drops and ordinary liquids, centrifuge accelerations as large as 5000g are considered.

The system is composed of ten physical parameters, which are the densities of the three phases ρ i (i = 1, 2, 3), their dynamic viscosities μ i (i = 1, 2, 3), the surface tensions γ 12 and γ 13 between phases 1 and 2 and between phases 1 and 3, respectively (assuming that phases 2 and 3 will never be in direct contact allows us to disregard γ 23 ), the droplet diameter d (or radius R), and the centrifugal acceleration a c = r i ω 2 evaluated at the interface position, r i being the radial position of the interface in the cylindrical coordinate system centered on the axis of rotation as sketched in Fig. 1 and ω the angular speed of rotation. As these ten quantities involve three fundamental units (mass, time, and length), the theorem of Vaschy-Buckingham states that the problem can be described with seven independent dimensionless parameters. In this study, we choose the density ratios ξ 12 = ρ 2 ρ 1 -1 and ξ 13 = ρ 3 ρ 1 -1, the viscosity ratios λ 12 = μ 2 μ 1 and λ 13 = μ 3 μ 1 , the Bond number Bo 12 = (ρ 2 -ρ 1 )a c R 2 γ 12 of the droplet immersed in phase 1, which describes its deformability by comparing the centrifugal force applied on the droplet to its interfacial tension, the Bond number Bo 13 = (ρ 3 -ρ 1 )a c R 2 γ 13 of the interface, and the Archimedes number Ar = ρ 1 (ξ 12 a c ) 1/2 R 3/2 μ 1 , which describes the flow regime in phase 1, which is merely a Reynolds number based on the gravitational velocity (ξ 12 a c R) 1/2 . Other useful nondimensional numbers can be used, such as the drop Reynolds number in phase 1, Re(t ) = ρ 1 V (t )d μ 1 , where V (t ) is the velocity of the droplet center of mass at a given time, and its Weber number in phase 1, which compares the inertial stress responsible for its deformation over the interfacial tension stress, We(t ) = ρ 1 V 2 (t )d γ [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF] .

The problem of a sphere settling through a liquid-liquid interface (with ρ 1 and ρ 3 of the same order of magnitude, contrary to the problem of splashing, where a sphere impacts a free surface) was addressed in the literature through theoretical models, numerical simulations, and experiments, gravity generally being the driving force in most of the studies. The different investigations have generally addressed one of the following questions: (i) the film drainage problem, which corresponds to a sphere standing at the interface in a quasistatic equilibrium, or (ii) the tailing regime leading to the entrainment of a column of fluid, where the sphere crosses the interface with a finite velocity, its crossing being possible even though the film drainage process is not achieved.

In the case of solid particles, the two configurations of (i) the film drainage mode (see the study of Hartland [START_REF] Hartland | The approach of a rigid sphere to a deformable liquid/liquid interface[END_REF]) or (ii) the tailing mode (see the study of Dietrich et al. [START_REF] Dietrich | Dynamical deformation of a flat liquid-liquid interface[END_REF]) were reproduced in the simulations of Geller et al. [START_REF] Geller | The creeping motion of a spherical particle normal to a deformable interface[END_REF] and presented in a regime diagram experimentally obtained by Jarvis et al. [START_REF] Jarvis | Experiments on the low-Reynolds-number settling of a sphere through a fluid interface[END_REF] and depending on both the interface Bond number Bo 13 and the viscosity ratio λ 13 . Generally, small and light spheres remain trapped at the interface. To analyze the quasistatic problem, Maru et al. [START_REF] Maru | Behaviour of a rigid sphere at a liquid-liquid interface[END_REF] proposed a force balance model, between the particle weight, driving its crossing, and surface tension and buoyancy forces opposing it, leading to the prediction of a minimal sphere radius as the crossing condition, validated by experimental observations. A similar critical condition was established by Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF] in the form of a minimum relative density contrast ξ 12 ξ 13 required to cross an interface under static conditions, at a given Bond number Bo 13 . From experimental observations with several particles and pairs of fluids, these authors found that this relationship successfully predicted, in most of the cases, if a particle would cross the interface or be trapped on it. This criterion was validated in quasistatic regimes, where a settling particle is stopped for a long time at the interface, and in dynamic conditions, where the sphere velocity remains finite during crossing. In only a few cases was the crossing possible, due to dynamic effects, whereas the quasistatic balance predicted that the particle should not cross, associated with a strong deceleration of the particle at the interface; the crossing mechanism was due in these cases to a supplementary history force [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF] issued from the collapse of the particle wake at the interface, which made it possible to overcome the resisting forces to crossing. In the tailing regime, Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF] have shown that the entrained liquid column geometry strongly depends on λ 13 and the Archimedes number Ar based on properties of phase 1, due to the fact that the column keeps the footprint of the wake developed in phase 1. The entrained tail stretches inside phase 3 and eventually breaks. Different column rupture modes have been classified by Aristoff and Bush [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] and characterized by Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF]: deep seal pinch-off, corresponding to a column detachment close to the particle in the case where buoyancy effects are dominant, and shallow pinch-off, where the column detachment first occurs close to the position of the interface when driven by the interfacial tension. Column axisymmetry is ensured provided that Ar < 55 [START_REF] Fabre | The steady oblique pass of buoyancy-driven disks and spheres[END_REF]. In very high inertial conditions (Ar > 55), the development of instabilities has been observed [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF][START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF] and leads to configurations of tailing with surrounding corollas or even full three-dimensional tails with possible fragmentation due to the strong shearing from the boundary layer developing around the tail when Ar is O(100) and λ 13 < 1. In all cases, once the entrained column breaks up, a part of the fluid of phase 1 covers the settling particle. At small ξ 13 ( 0.1) and with continuous phases of contrasted viscosities λ 13 0.02, the coating volume experimentally measured by Pitois et al. [START_REF] Pitois | Franchissement d'interface et enrobage d'une sphère[END_REF] (particles falling in Stokes regime in phase 1) and Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF] scales linearly with a modified Bond number ( ξ 12 ξ 13 -1)Bo 13 obtained from a force balance on the coated particle in the limit of small Ar; note that, in these cases, the encapsulating volumes are of same order of magnitude as that of the particle. At larger liquid density contrasts (ξ 13 > 0.1), the coating volume was observed to be highly sensitive to both Bo 13 and ξ 12 /ξ 13 in a nontrivial way.

The crossing of a liquid-liquid interface by bubbles is more complex due to bubble deformation. Shopov and Minev [START_REF] Shopov | The unsteady motion of a bubble or drop towards a liquid-liquid interface[END_REF] performed numerical simulations of deformed bubbles crossing an interface, at small and moderate Reynolds numbers, focusing on the bubble and interface deformations and on the film drainage dynamics rather than on the tailing mode. They have shown that, at very low Weber and Reynolds numbers, a bubble can be deformed into a prolate shape (i.e., elongated in the direction of motion) during the passage through the interface in cases at low λ 13 , a result confirmed by the simulations of Manga and Stone [START_REF] Manga | Low Reynolds number motion of bubbles, drops, and rigid spheres through fluid-fluid interfaces[END_REF], whereas at higher Weber and Reynolds numbers, the action of inertial forces and the interaction with the interface lead to oblate bubbles with the formation of a concavity at the rear and spherical cap shapes during the crossing. Interface crossing in inertial regimes was also experimentally investigated by Dietrich et al. [START_REF] Dietrich | Passage of a bubble through a liquid-liquid interface[END_REF], who provided the relationship of the crossing time of the interface by a bubble as a function of the ratio of the terminal velocities inside the two continuous phases. Once the tail is formed in the column, their particle image velocimetry measurements have shown the coexistence of a fluid motion entrained in the bubble wake and an opposite flow driven by gravity. The crossing in the tailing mode in the case of a bubble was also considered by Bonhomme et al. [START_REF] Bonhomme | Inertial dynamics of air bubbles crossing a horizontal fluid-fluid interface[END_REF], both experimentally and numerically under inertial conditions, in a wide variety of bubble shape configurations (from spherical to toroidal). A mapping of the bubble shapes and entrained column geometries was provided as a function of both Bo 13 and Ar. Small bubbles are slowed down and even stopped at the interface without crossing (film drainage configuration, the crossing being possible only provided the drainage process is achieved), whereas larger spherical cap bubbles at larger Bo 13 generally cross the interface easily. It has been found that the volume entrained is also larger with such spherical cap bubble shapes since they offer a larger cross section to the interface than spheroidal bubbles, for a given gas volume. In the experiments of Emery et al. [START_REF] Emery | Flow regimes and transition criteria during passage of bubbles through a liquid-liquid interface[END_REF], the crossing by a single bubble shows additionally that, in the tailing mode, the column of liquid entrained is longer in the case where the bubble velocity does not change much during the crossing of the interface; the tail is observed to remain connected a long time before its rupture, in some cases the liquid shell covering the bubble breaking before the column. The latter study has investigated the crossing of a stream of bubbles, giving a map of the different flow regimes with the possible formation of clusters. The experimental results of Manga and Stone [START_REF] Manga | Low Reynolds number motion of bubbles, drops, and rigid spheres through fluid-fluid interfaces[END_REF] also illustrate the case of vertically aligned bubbles, by showing that the following bubble in a train experiences less resistance during the crossing, leading to a more elongated shape as compared to the previous bubble.

Finally, the problem of droplets crossing a liquid-liquid interface is even more complex because it involves more parameters as compared to the case of particles or bubbles, due to the droplet viscosity and deformability. Only a few studies are available, most of them addressing the film drainage problem in the context of drop-interface coalescence under quasistatic conditions. In this configuration, Hartland [START_REF] Hartland | The effect of circulation patterns on the drainage of the film between a liquid drop and a deformable liquid-liquid interface[END_REF][START_REF] Hartland | The profile of the draining film between a fluid drop and a deformable fluid-liquid interface[END_REF] studied experimentally the profile of the draining film beneath a droplet approaching a liquid-liquid interface and derived expressions of the drainage dynamics in the case of either mobile or immobile interfaces, a more complete theoretical analysis being further provided by Jones and Wilson [START_REF] Jones | The film drainage problem in droplet coalescence[END_REF], who clarified the possible narrowing of the drainage film at its periphery (dimple formation). Simulations by Chi and Leal [START_REF] Chi | A theoretical study of the motion of viscous drop toward a fluid interface at low Reynolds number[END_REF] at low Reynolds number were also able to consider the influence of the drop internal circulation on the film drainage dynamics, through the viscosity ratio λ 12 (=λ 13 in their particular conditions). They confirmed that different types of film geometry can be observed upon the arrival of the drop at the interface: a film with a minimum thickness at the top indicating a rapid drainage at low λ 12 , a film with a uniform thickness when λ 12 is of the order of unity, and a film with minimum thickness at its periphery involving the formation of a dimple at high λ 12 . When the droplet travels a distance of several radii beyond the interface, i.e., when it actually crosses the interface, Manga and Stone [START_REF] Manga | Low Reynolds number motion of bubbles, drops, and rigid spheres through fluid-fluid interfaces[END_REF] did a parametric numerical study in low-Reynolds-number regimes in the first phase. The influence of four nondimensional numbers λ 12 , λ 13 ,Bo 13 , and Bo 12 on the drop and interface deformations and film drainage rate was addressed. As for bubbles, at low λ 13 , drops have been observed to undertake a prolate shape during the interface crossing; at large Bo 12 , they deform by developing either an elongated tail or a small cavity at the back. Concerning the film drainage, it has been found that the higher the droplet viscosity (the higher λ 12 ), the thicker the coating film of phase 1 around the drop, the slowest drainage rate occurring when λ 12 1 λ 13 . However, due to the limitation of the numerical approach, calculations could not capture the long tail dynamics and its rupture. Studies dedicated to investigating the tailing mode with droplets until breakup of the column are lacking.

In this paper, the crossing of a liquid-liquid interface by a droplet submitted to a centrifugal field is numerically investigated, by focusing on the tailing regime. Two types of drops are considered, both solidlike drops mimicking rigid particles (with a high internal viscosity and surface tension) and deformable drops at different viscosity ratios λ 12 , in a wide range of flow parameters and phase properties, which however correspond to common liquid phases and oil-water interfaces. In all cases, droplet motion lies in inertial regime (7 Re 160, based on the settling velocity in phase 1).

The paper is structured as follows. In Sec. II the numerical level-set method used in the three-phase system is briefly presented and some validation test cases are discussed to compute the terminal velocity of nondeformed droplets, the deformation of drops in inertial regimes, and, in three-phase configurations, the crossing criterion for solidlike drops in static conditions. Then the results obtained in dynamic crossing conditions (i.e., with an impacting velocity) are presented and discussed in Sec. III, starting with the effect of the drop velocity at the interface on the condition for crossing, followed by an analysis of the drop shape and velocity during crossing. Finally, the length of the entrained column and the volume encapsulating the droplet are scaled as a function of parameters characterizing drop inertia.

II. NUMERICAL METHODS AND VALIDATIONS A. Simulation method

Axisymmetric direct numerical simulations of droplets traveling through liquid-liquid interfaces are carried out using the numerical code DIVA. In view of solving a three-phase flow involving the dynamics of two interfaces, the numerical approach is based on a one-fluid formulation of the Navier-Stokes equations, and the level-set method is used to capture the two interfaces on a Cartesian grid that does not follow the complex interface shapes. This numerical method is the natural extension for three phases of that described by Tanguy et al. [START_REF] Tanguy | A level set method for vaporizing two-phase flows[END_REF], Lalanne et al. [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF], or Rueda-Villegas et al. [START_REF] Rueda-Villegas | A ghost fluid/level set method for boiling flows and liquid evaporation: Application to the Leidenfrost effect[END_REF] in the case of two-phase flows, with several validation tests of the dynamics of deformed droplets and bubbles. The method is briefly presented here.

The method consists in solving the Navier-Stokes equations for an incompressible flow to compute the velocity field u and the pressure field P,

ρ ∂ u ∂t + ( u • ∇ ) u = -∇P + ∇ • (2μD) + F , (1)
∇ • u = 0, (2)
where ρ and μ are the density and dynamic viscosity of the fluid which is present at the considered mesh cell center, D is the rate of the deformation tensor defined as

D = ∇ u+ ∇ u T 2
, and F is the volume force. These equations are used together with a level-set method [START_REF] Osher | Fronts propagating with curvature-dependent speed: Algorithms based on Hamilton-Jacobi formulations[END_REF][START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] for capturing the two interfaces 1 and 2 , which are the droplet and the liquid-liquid interfaces, respectively. Numerically, these interfaces are defined as the zero level of two respective distance functions φ 1 and φ 2 that represent the algebraic distance between a mesh point and each interface, allowing us to locate in which phase any mesh cell center lies. The evolution of φ 1 and φ 2 is computed from two advection equations, for i = 1, 2,

∂φ i ∂t + u • ∇φ i = 0. ( 3
)
For each of these two level-set functions, an algorithm of redistancing [START_REF] Sussman | A level set approach for computing solutions to incompressible two-phase flow[END_REF] is used to ensure that φ i (i = 1, 2) remain distance functions at each time step. Different jumps of density, viscosity, and pressure, denoted by square brackets in the equations, exist at both interfaces; in particular,

[P] 1 = γ 12 κ 1 + 2[μ] 1 ∂u n ∂n , ( 4
) [P] 2 = γ 13 κ 2 + 2[μ] 2 ∂u n ∂n (5)
give the balance of the normal stresses at each interface involving pressure and viscous normal stress discontinuities, where κ i is the local curvature of the interface and u n the normal velocity of the fluids at the interface. The ghost fluid method [START_REF] Fedwik | A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method)[END_REF] is employed to take into account these discontinuities, as jumps of pressure during the resolution of the Navier-Stokes equations [START_REF] Lalanne | On the computation of viscous terms for incompressible two-phase flows with level set/ghost fluid method[END_REF].

In the following simulations, the force F driving the drop motion will be either normal gravity ( F = g) or a higher acceleration similar to that existing in a centrifugal field, i.e., linearly increasing with the drop position. In the latter case, it is expressed as F = ω × ω × r, where ω is the rotation vector and r the position vector regarding the axis of rotation (see Fig. 1). Note that such an inertial force is not sufficient to fully simulate the flow in a centrifugal device from a rotation frame, as a Coriolis force should also be considered for that purpose, which would require us to use three-dimensional simulations; however, we can consider that an axisymmetric simulation, which includes only the centrifugal force, will be a good approximation of the flow in the case where the Coriolis force has a negligible impact on drop deformation, influencing only slightly its trajectory. All simulations are thus carried out under the axisymmetric assumption, in a frame moving with the droplet, which allows us to keep it at a fixed position in the simulation domain so as to maintain a high grid resolution around the drop.

B. Post-treatment of film volume

In the simulation of the dynamic interface crossing resulting in drop encapsulation, the film volume covering the droplet is evaluated as the post-treatment at the first instant where the capsule is formed. The method consists in evaluating the volume V e enclosed by the interface between phases 1 and 3 around the droplet (see Fig. 13 for screenshots from the simulations), based on the level-set function φ 2 , from which the droplet volume V drop is subtracted in order to obtain the film volume V f = V e -V drop . As the post-treatment, V e is computed through the volume integration of a regularized Heaviside function H, equal to 1 inside the film contour defined by a minimum numerical thickness 0 = 1.5 x (with x the smaller mesh cell size) as

H (φ 2 ) = ⎧ ⎨ ⎩ 0 i f φ 2 < -0 1 2 [1 + φ 2 0 + 1 π sin( πφ 2 0 )] if |φ 2 | 0 1 i f φ 2 > 0 .
As explained in Sec. III B 4, for most of the cases of dynamic crossing, the contribution to V f from the bottom part of the droplet is generally dominant, which is a well-resolved zone, contrary to the contribution at the top, which is of negligible thickness (see Figs. 13 and14). The uncertainty of the computed V f values needs to be assessed and will now be discussed.

First, concerning the level-set method accuracy on the mass conservation, whereas the numerical method ensures the conservation of the droplet phase in time, even in the case of high deformation, that of the covering film volume is slightly affected, in particular because of the underresolution of the film at the top of the droplet. In order to get an estimate of the resulting mass loss, some computations have been continued for a longer time than the time required to encapsulate the drop, allowing us to compute V f at several instants when it should remain constant. The computations show that the mass loss of the film reaches 2% after a supplementary travel distance of six radii for the capsule. Such a mass loss is not representative of the error of the film topology prior to detachment, but gives an idea of the uncertainty. To minimize it, in the following simulations, V f is computed as soon as the droplet is encapsulated.

Second, in order to test the sensitivity of V f on the grid resolution (which has an influence on the mass conservation in the numerical method, the resolution of the flow in the film, and the posttreatment used to compute V f ), a mesh convergence study has been carried out, in case D 2 from Table II, by using grids with 48, 72, and 96 points per radius. The difference in volume computation was less than 0.7% between the different grids; therefore, a grid resolution of 48 points per radius has been chosen for the following simulations. Note, however, that, despite this very good accuracy of V f , all these grid resolutions are still insufficient to accurately capture the very thin thickness of the film at the front of the droplet. This is a limitation of such Eulerian numerical simulations to study this multiscale problem, the film thickness becoming several orders of magnitude lower than the characteristic length scales of the flow around the drop, as also discussed in the numerical study of Bonhomme et al. [START_REF] Bonhomme | Inertial dynamics of air bubbles crossing a horizontal fluid-fluid interface[END_REF].

Finally, the uncertainty of V f comes mainly from the contribution of the film volume from the top part of the droplet. In order to quantify an error of the computed V f , a criterion is defined

= 2π R 2 x V f
based on the size of one mesh cell x = R/48 around the northern hemisphere of the drop where the mesh resolution is not sufficient. Such a criterion probably overestimates the error from the numerical simulation and the post-treatment process, but allows us to include error bars in Fig. 13 for the computed V f . The accuracy will be higher for the highest V f .

C. Validation of the numerical model

In this section we validate the numerical method by (i) performing simulations of the gravitydriven motion of a liquid droplet in a stagnant liquid phase and comparing the results to predictions of the drop velocity and its deformation and (ii) performing simulations of a three-phase flow configuration of a solidlike droplet crossing or floating at a liquid-liquid interface in quasistatic conditions and comparing the result with existing data for solid particles.

Two-phase flow validations

As a validation of our methodology, we compute the rise of a droplet in a stagnant liquid (ρ 2 /ρ 1 = 1.29 and μ 2 /μ 1 = 1.67), under constant gravitational acceleration g at two different Reynolds numbers Re T = ρ 1 u T d μ 1 , where u T is the terminal velocity: Re T = 20 (case 1) and Re T = 180 (case 2). For each case, the value of g is adjusted so as to reach the desired u T and the Weber

number We 12 = ρ 1 u 2 T d γ 12
based on u T is maintained low enough to keep the droplet spherical during its motion. The terminal velocity can be predicted through a force balance between the drag force and the buoyancy force in the steady-state regime

0.5πR 2 ρ 1 C D u 2 T = (ρ 2 -ρ 1 ) 4π R 3 3 g, ( 6
)
where C D is the drag coefficient for spherical droplets, expressed as, according to Rivkind and Ryskin [START_REF] Rivkind | Flow structure in motion of a spherical drop in a fluid medium at intermediate Reynolds numbers[END_REF], ( 24 Re + 4Re -1/3 ) + 14.9Re -0.78 ], valid with a maximal accuracy of 7% for Reynolds numbers + lower than 200. The mesh is uniform in both cases and grid resolutions of 8, 16, and 32 nodes per radius are evaluated. Figure 2 shows the Reynolds-number evolution as a function of the normalized distance traveled by the droplet for each case and with different mesh resolutions. In both cases, the numerical code shows a very good convergence above 16 points per radius, as demonstrated by the overlap of the curves obtained for more refined grids. The Reynolds number converges towards a limit value Re T at distances z large compared to the droplet diameter. For case 1, the prediction of u T by Rivkind and Ryskin's correlation is perfectly consistent with the results of numerical simulations with a discrepancy of 0.4%. For case 2 it is around 9%. Since the numerical results are converged in terms of mesh resolution, it can be concluded that the boundary layers around the droplet are well resolved in this range of drop Reynolds number with our numerical tool, which is the range of interest in this study.

C D = 1 1 λ 12 [λ 12
In order to validate the computation of drop deformation by the numerical tool, simulations of rising droplets in another liquid (1.11 ρ 2 /ρ 1 1.62 and 0.12 μ 2 /μ 1 50) are carried out under a constant gravitational acceleration. Results are obtained with mesh resolutions of 48 mesh points per drop radius. The drop Reynolds number is large (7 Re T 163), leading to a drop deformation controlled by the inertial forces inside the continuous phase, then characterized by the Weber number We 12 . The drop deformation is described through the aspect ratio χ at steady state, once the deformation no longer evolves, which is defined as χ = a/b, where a is the larger axis of the ellipsoidal drop shape and b the smaller axis. Figure 3 shows the measured drop deformation as a function of We 12 : χ varies in a large range in the simulations, until 3.3, and is found to increase with We 12 . For weakly deformed droplets (χ 1.3), the results are in good agreement with the correlation proposed by Wellek et al. [START_REF] Wellek | Shape of liquid drops moving in liquid media[END_REF], where there is not any influence of the viscosity ratio λ 12 on the steady deformation. In the case of more deformed droplets, our numerical results are also fully consistent with the simulation predictions of Bäumler et al. [START_REF] Bäumler | Drop rise velocities and fluid dynamic behaviour in standard test systems for liquid/liquid extraction-experimental and numerical investigations[END_REF] at a viscosity ratio λ 12 = 2.4 for oblate droplets with an aspect ratio as high as χ = 2.8. At such large deformations, it is observed in Fig. 3 that the drop viscosity strongly influences the deformation: For a given We 12 , the higher the λ 12 , the lower the deformation. This was explained by Hinze [START_REF] Hinze | Fundamentals of the hydrodynamic mechanism of splitting in dispersion processes[END_REF]: In contrast to bubbles, the viscous dissipation associated with the internal circulation in droplets consumes a part of the outer flow energy responsible for deformation. On the basis of these quantitative and qualitative comparisons, axisymmetric simulations of drop shapes are validated, even in the case of highly deformed droplets. FIG. 3. Simulation of droplets with different values of λ 12 : aspect ratio as a function of the Weber number We 12 , including numerical results obtained with the same numerical tool from Lalanne et al. [START_REF] Lalanne | Effect of rising motion on the damped shape oscillations of drops and bubbles[END_REF], experimental results from Wellek et al. [START_REF] Wellek | Shape of liquid drops moving in liquid media[END_REF], and numerical results from Bäumler et al. [START_REF] Bäumler | Drop rise velocities and fluid dynamic behaviour in standard test systems for liquid/liquid extraction-experimental and numerical investigations[END_REF] with another simulation tool.

Three-phase flow: Static crossing

For validation purposes, we consider here a three-phase flow configuration of a solidlike droplet crossing a liquid-liquid interface in static conditions. We call a condition static when the solid particle is directly placed at the position of the interface without initial velocity. A force balance model able to predict the critical crossing conditions under such a particular assumption for a rigid sphere was developed by Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF] with several parts common with previous publications (e.g., [START_REF] Maru | Behaviour of a rigid sphere at a liquid-liquid interface[END_REF][START_REF] Bonhomme | Inertial dynamics of air bubbles crossing a horizontal fluid-fluid interface[END_REF][START_REF] Malmazet | Coalescence of contaminated water drops at an oil/water interface: Influence of micro-particles[END_REF][START_REF] Vella | The load supported by small floating objects[END_REF]). The model predicts if the solid particle will float or cross the interface, by considering the balance between the force favoring the crossing (apparent weight) and the resisting forces (surface tension force of the liquid-liquid interface and hydrostatic pressure). In two different asymptotic regimes depending on Bo 13 , the critical conditions are given in a nondimensional form according to

ξ 12 ξ 13 f (Bo 13 ) ( 7 )
for Bo [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] 1, with

f (Bo 13 ) = 3 2 Bo 13 + 1 2 + 3 4 log 10 4 √ Bo 13 -0.577 , ( 8
)
and for Bo [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] 1, with

f (Bo 13 ) = 3 2 Bo 13 sin 2 k + 1 4 (2 + 3 cos k -cos 3 k) + 3 4 2 Bo 13 1/2 sin 2 k, ( 9
)
where k = 2(2 Bo 13 ) -1/4 . In order to simulate the crossing of a rigid particle with the numerical tool being used in our simulations, configurations of static crossing are simulated with solidlike droplets [λ 12 = 50 and small droplet Bond number Bo 12 ∼ O(10 -2 )] for different values of ξ 12 ξ 13 and Bo 13 . In these simulations, the mesh resolution is 48 grid points per drop radius. The numerical results show excellent consistency with theory, as illustrated in Fig. 4, where the theoretical criterion for crossing and the corresponding numerical simulation cases are reported. These results validate both the three-phase flow extension of the numerical code and the use of the solidlike droplet approximation as a model of the solid particle.

III. INTERFACE CROSSING IN DYNAMIC CONDITIONS

In this section we consider axisymmetric simulations of interface crossing in inertial conditions by either solidlike droplets S i , simulated as highly viscous (λ 12 = 50) and nondeformed droplets (both Bo [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF] 1 and We 12 1), or deformable droplets D i with variable λ 12 . As shown in Fig. 1, the driving force is induced by a high acceleration, which increases linearly with the drop position like in a centrifugal field. Simulations are carried out in the frame moving with the droplet, of size 6R × 13R or 19R, with a resolution of 48 mesh points per radius based on the validation tests. In what follows, the conditions and dynamics of interface crossing and film entrainment are investigated.

In terms of nondimensional parameters, several values of ξ 12 are considered (1.67 ξ 12 21,

ξ 13 ξ 13
which correspond to 0.0105 ξ 12 0.2941 and 0.0005 ξ 13 0.176). For each case, the Bond number Bo 13 is taken to be high enough for interface crossing, the study being performed from moderate to high inertia 5.3 Ar 70 (the maximal value of Ar is 55 for solidlike droplets). For deformable drops, the drop Bond number Bo 12 lies in between 1.3 and 6.25 and the viscosity ratio λ 12 is varied in a wide range (0.1 λ 12 50). In most cases, λ 13 = 1 and in a few cases λ 13 < 1 (down to 0.05). For solidlike droplets, it has been verified that the tangential velocity along the interface is always small compared to the drop rise velocity u T ; the maximal value never exceeds 10% of u T , leading to a negligible fluid velocity both at the interface and inside the solidlike droplet. FIG. 4. Mapping of the crossing configuration for solidlike droplets. The red solid line corresponds to the expression (8) and the green solid line corresponds to the expression [START_REF] Geller | The creeping motion of a spherical particle normal to a deformable interface[END_REF], with an assumed red dashed line that connects the two theoretical predictions. Closed squares correspond to a success in crossing, while open squares are associated with a fail in crossing, in the static case.

The physical parameters and nondimensional numbers of all simulations are reported in Tables I and II, respectively.

A. Crossing conditions with additional inertia

Here the effect of drop impact velocity on the crossing condition is examined. In a quasistatic regime, the critical crossing condition is given by Eqs. ( 8) and ( 9) in two limits Bo [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] 1 and Bo [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] 1. In dynamic conditions, the droplet is released far from the interface (at a distance of three diameters in the simulations) and has time to develop its wake in phase 1, the drop motion being inertial in the conditions studied. When the droplet approaches the interface, it decelerates due to the presence of the interface, reaching a minimum value U min while crossing. For each case, the minimal Reynolds number Re min computed from this velocity is reported in Table II. In most cases, Re min is positive, but in some cases, Re min falls to 0 when the drop (deformable or not) is stopped by the interface. Further, U min can even reach negative values for bouncing drops. For a noncrossing case in static conditions, inertia could be expected to favor the crossing. Figure 5 compares the theoretical critical condition for crossing given by the static theory with the simulation results in dynamic cases at different Ar values for both solidlike and deformable droplets.

For the simulations at different values of ξ 12 ξ 13 in the range investigated, the critical condition for crossing is found to be consistent with the quasistatic theoretical predictions. This conclusion stands for both the solidlike and the deformable droplets, showing the strong interest of the static theory, given by Eqs. ( 7)-( 9) and developed for spherical particles to make a reliable prediction of whether or not a droplet will cross an interface. Moreover, from the numerical results of Fig. 5, the crossing condition is shown to be very slightly affected by the impact velocity. Actually, the value of Ar is probably not crucial to understand the role of additional inertia. Indeed, in the inertial regime and for liquid-liquid systems of close density, the dynamic pressure force exerted by the drop on the interface is approximately equal to its apparent weight; whatever the value of Ar, the maximum force brought by inertia therefore simply doubles the particle weight. Hence, an equivalent Bond number Bo 13eq can be defined to take into account the maximum dynamic force due to inertia, which finally turns out to consider a drop of the same properties but with a volume multiplied by 2, i.e., a diameter multiplied by 2 1/3 , leading to Bo 13eq = 2 2/3 Bo 13 = 1.59Bo 13 (note that the correction factor is small). Making use of this equivalent apparent weight, the solidlike case at Ar = 22 in Fig. 5, which is subcritical (i.e., the noncrossing case), has a Bo 13 approximately equal to 0.6. Multiplying Bo 13 by 1.59 leads to Bo 13eq = 0.95, which still stands below the critical condition on the static crossing condition curve: The particle is not crossing. The same evaluation can be achieved for the deformable drop case at Ar = 28 and Bo 13 < 0.04. In this case, Bo 13eq = 0.06, which remains smaller than the critical value, close to 0.08; here again, the droplet is not crossing.

The case Ar = 34, corresponding to a solidlike droplet, is interesting because its equivalent Bond number (Bo 13eq = 0.087) is close to (slightly above) the critical value (around 0.08) predicted by the static crossing theory. In this case, the drop is submitted to several oscillations and it does not cross the interface at the end of the simulation. Such a behavior indicates that this point is close to the critical condition of crossing as predicted by the equivalent Bond number. Running a simulation over a very long time could have shown that crossing was finally possible. Finally, these different examples of crossing conditions in dynamic cases emphasize the strong interest of the static theory, which is shown to be able to predict quite well if a droplet succeeds or fails to cross the interface, even in the case of a deformable droplet which arrives at the interface with a nonzero velocity. For such liquid droplets translating in an inertial regime, a rough approximation to take into account the dynamic force that could favor the crossing allows us to consider a higher equivalent Bond number Bo 13eq , the crossing being possible when the latter overcomes the minimal value required which is predicted by the static theory. Note that the effect of inertia discussed here only concerns the final state, i.e., crossing or noncrossing, and disregards the dynamics of the interaction of the drop with the interface, which naturally is strongly dependent on Ar and involves transient inertial forces as discussed by Pierson and Magnaudet [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF].

In the following section, simulation results of crossing cases S i (solidlike drops) and D i (deformable drops) are presented and discussed.

B. Dynamics of interface crossing

Phenomenology

For solidlike droplets, Fig. 6 shows the evolution of the velocity of the particle centroid as a function of its travel distance across the interface and the corresponding image sequence, for cases S 1 (Ar = 53) and S 2 (Ar = 20), with the same density and viscosity ratios (Table II). In both cases, the droplet accelerates in phase 1 and the droplet wake develops with recirculations, as expected [START_REF] Johnson | Flow past a sphere up to a Reynolds number of 300[END_REF] at these Reynolds numbers (the maximum values of the Reynolds number in phase 1, Re max , are 160 8)and ( 9). The squares and circles correspond to simulations of solidlike droplets and deformable droplets, respectively, in the dynamic case, with the points labeled by the Archimedes number Ar. Symbols are closed in the case of crossing in the simulation and open otherwise. and 48, respectively). In these cases, the distance of the drop release to the interface is not sufficient to reach the point where viscous effects are fully established (note that, even with a larger initial distance between the droplet and the interface, a steady state is not expected because the centrifugal acceleration increases linearly with the distance from the rotation axis; only a quasisteady state could be reached, with a drop acceleration that becomes low but does not vanish). At a distance from the interface of the order of the drop radius (close to point 1), the droplet velocity reaches a maximum U max and then decelerates due to the presence of the interface (point 2 corresponds to the drop arrival at the position of the plane interface). Between points 2 and 3, drainage of phase 1 develops in the thin film between the top of the droplet and the interface, giving rise to a lubrication flow in this film, as the droplet continues its rising motion, pulling a column of phase 1 in its wake. At point 3, the droplet velocity reaches a minimum value U min (from which Re min = ρ 1 U min d/μ 1 is computed). In both cases, while the droplet velocity reaches an extremum (U max and U min close to point 1 and at point 3, respectively), all of the forces exerted on the droplet center of mass are at equilibrium. Then, from point 3 to point 5, the droplet accelerates inside phase 3 with a different rate in the two cases considered, entraining a liquid column of phase 1 which extends up to a maximum length L max before it breaks. For the two cases considered here, the column first detaches at the bottom close to the interface position (point 6) and then at the rear of the droplet (point 7), leaving the droplet coated by a volume of phase 1 rising in phase 3. Due to gravitational effects, the volume of the lighter phase 1, coating the droplet, moves towards its rear, thinning up strongly at the top. It can be noticed that L max is larger in case S 1 than in case S 2 , i.e., when inertia is higher.

Let us now consider the case of deformable droplets at Ar = 53 and the same Bond number of the interface Bo 13 = 0.134, similarly to S 1 (solidlike droplet), for two different viscosity ratios λ 12 = 50 (case D 1 ) and λ 12 = 5 (case D 2 ). Case D 1 corresponds to a very viscous but deformable droplet and case D 2 corresponds to a deformable droplet of lower internal viscosity.

The droplet dynamics during its rising motion towards the interface is first analyzed. Simulation results of cases D 1 and D 2 both in the presence of the liquid-liquid interface (three-phase system) and without the interface (two-phase system) are compared in order to evaluate the coupling between the drop velocity and its shape.

Figure 7 displays the evolution of the instantaneous Reynolds number and aspect ratio of the two droplets. Quasisteady conditions are not reached in the three-phase flow simulation before the liquid-liquid interface. In particular, if the velocity is close to its terminal value, the aspect ratio is still increasing, based on the results from the two-phase flow simulations. The slow drop deformation dynamics is due to viscosity effects from both phases 1 and 2, in the same way as the damping rate of eigenmodes of drop shape oscillations [START_REF] Miller | The oscillations of a fluid droplet immersed in another fluid[END_REF]: For the same μ 1 , the higher the drop viscosity, the slower the drop shape response to the deforming stress. By analyzing the steady-state conditions, the D 1 droplet is found to have a larger velocity (larger Re and We 12 )but lower deformation than the D 2 droplet. Indeed, as already emphasized in Fig. 3, deformation of a viscous droplet in inertial conditions is a function not only of We 12 but also of the viscosity ratio λ 12 , which is ten times larger for D 1 , resulting in a less flattened drop shape in quasisteady conditions. Concerning the drop rising velocity, even if Ar is the same in both cases (same acceleration), D 1 rises faster than D 2 , a result which probably arises from two competing effects acting on the drag force: As μ 2 is higher for D 1 , dissipation of the energy provided by gravity is increased for D 1 , whereas the aspect ratio, of 2.4 for D 1 lower than that of D 2 (3.4), results in a larger drag on D 2 than on D 1 , the influence of drop shape being observed to dominate since the velocity of D 2 is finally smaller. Note also that some velocity and shape oscillations are observed in the quasisteady regime, of larger amplitude in case D 2 , probably due to the fact that Ar is close to the critical value of 55 where a path instability appears (known for a solid particle motion [START_REF] Fabre | The steady oblique pass of buoyancy-driven disks and spheres[END_REF]), even though such a threshold could be different in a computation which assumes axisymmetry. All these comments make us understand that when these droplets arrive close to the liquid-liquid interface in the three-phase flow simulation, their shapes can strongly differ and can still be deforming despite close velocities. Now, in the presence of the liquid-liquid interface, Fig. 8 shows the variation of the Reynolds number along with that of the drop aspect ratio as a function of the drop position, for cases D 1 and D 2 as well as for the solidlike droplet S 1 .

Before the interface, due to the strong effect of deformation on the rise velocity, the highly deformed droplets D 1 and D 2 both have a smaller maximum velocity (at point 1) than the solidlike one S 1 . Between points 1 and 3, the presence of the interface makes the velocity decrease at a distance of about 1R before the interface, as clearly visible from the curves of Fig. 7. The droplet shape continues flattening until point 2, where the action of external stresses deforms the droplet in the direction perpendicular to the acceleration, hence towards a spherical shape. At point 3, the interface is crossed and a column of phase 1 is entrained by the drop motion, like in the solidlike drop case. Between points 3 and 4, droplet D 2 , which reacts faster in deformation to the external flow due to a lower μ 2 , sees its deformation nearly vanishing (the minimum value of the aspect ratio is close to 1) before flattening again due to the acceleration in phase 3; the same dynamics is observed for D 2 , except that the droplet never retrieves a spherical shape during its deformation history. The drop acceleration at point 4 allows the column of phase 1 to be strongly extended.

Figure 9 shows the velocity fields of cases S 1 , D 1 , and D 2 at point 5, which is the last instant before the column breaks. In all cases, the entrained fluid column behind the droplet is lighter than the surrounding phase 3, so a part of the fluid in the column moves back towards the interface whereas another part is entrained in the droplet wake. This wake is already developed from the drop rising in phase 1, and in these examples, it is not really modified by the presence of the interface because the viscosity of phases 1 and 3 is the same (λ 13 = 1). The column is stretched and thinned by an upward and a downward flow, which leads to the presence of a zero velocity point. The rupture takes place at point 6 for cases D 1 and D 2 . In case D 1 , the detachment mechanism is similar to that described for cases S 1 and S 2 of Fig. 6, taking place first at the bottom of the column and then in the drop wake. However, in case D 2 , the column detachment first occurs behind the droplet. In all cases, the encapsulated droplet, detached from the interface, rises in phase 3 and a part of the column of phase 1 moves back to the interface, possibly breaking up into several droplets. The detachment in cases S 1 and D 1 is similar to the shallow pinch-off detachment mode, while that of case D 2 is similar to deep seal detachment mode, both reported by Aristoff and Bush [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF].

The velocity field of case D 2 shows a high internal circulation inside the droplet and a high tangential velocity at the interface comparable to the droplet velocity, which forces the interface to follow the deformation of the droplet and to stay close to it up to a high separation angle. Such a high tangential velocity and drop internal circulation (lower value of λ 12 ) favors a rapid film drainage along the interface and forces the column to detach earlier behind the droplet rear. In cases S 1 and D 1 , the droplet has a stronger internal viscosity leading to an approximately null tangential velocity and internal circulation, due to the continuity of tangential stresses at the interface. This does not favor the film detachment close to the droplet rear as in case D 2 and tends to increase the film volume finally covering the droplet. In case D 2 , the resulting film volume (V f = 0.58V drop )is smaller than in case S 1 (V f = 0.67V drop ) and much less than in case D 1 , where the film volume is approximately equal to the droplet volume.

Regarding the film drainage dynamics during the crossing, in the case where the drainage is slow (high values of λ 12 ), the grid resolution is fine enough to capture the drainage flow (see points 2 and 3inFigs.6 and 8). Obviously, this is not the case when fast drainage occurs. However, the film at the top of the droplet is expected to have a negligible influence regarding the interface crossing problem, since it has a very small thickness and the pressure across this film can be considered as constant (the vertical pressure gradient in the film is only hydrostatic based on the lubrication theory), and thus does not impact the force balance acting on the droplet, as already confirmed by the excellent agreement, in static conditions, between the theoretical predictions of Eqs. ( 8) and ( 9) (also neglecting the influence of the thin film) and experimental data with solid particles [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF]. This allows us to be confident about the results obtained by the simulations on the interface crossing problem in dynamic conditions.

To conclude this section, in crossing conditions at the same Ar and Bo 13 , we observe that a very viscous droplet (λ [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 2. Sphere and tail dynamics[END_REF] 1), even deformed, has a behavior similar to a solidlike droplet regarding column entrainment during the crossing and further detachment that takes place at the bottom of the column, whereas a moderate λ 12 results in a more deformed droplet and enhances the drainage rate in the film due to the nonvanishing tangential velocity at the interface, leading to a column detachment first occurring at the droplet rear and to a smaller volume encapsulating the droplet.

Crossing time

When the drop reaches the interface with a nonzero velocity, it is slowed down (sometimes significantly), as illustrated in Fig. 6, due to the resistance of the plane interface deformation. In view of the development of a continuous encapsulation process, the slowing time is an important quantity to take into consideration because it fixes an upper limit for the frequency of droplet arrival at the interface. Indeed, the time between two successive drops arriving at the interface must be significantly longer than a characteristic time related to the crossing dynamics.

In the case where the condition of crossing is not achieved, i.e., when the resisting forces of the interface overcome the drop inertial forces, the drop will bounce and its velocity will become negative before canceling; then the drop will finally stand below the interface. In such a condition, we have observed that the distance traveled by the drop above the level of the plane interface is always smaller than its diameter. We define a crossing time t cr as the time required for the droplet to travel a distance equal to its diameter d, after its center of mass reached the interface position, as illustrated in Fig. 10.

As shown in Fig. 11, for either solidlike or deformable droplets, it is found that t cr is well scaled by the arithmetic average of U max , the drop maximal velocity prior to the interface (close to the terminal velocity in phase 1), and U min , the minimal velocity reached after the slowing down when crossing the interface. This result is used to consider the linear behavior of the velocity in between the two extrema. The crossing time can therefore be predicted from the scaling of the two velocities U max and U min . This question will be addressed in Sec. III B 5.

Column maximal length

As clearly observed in Figs. 6 and8, inertia has an impact on the entrained column of phase 1 in phase 3 during the droplet crossing: The higher the velocity at the interface, the longer the entrained column. We define the maximal column length L max as the distance between the droplet center of mass and the position of column detachment at the bottom of the column (see Fig. 6), which always occurs even in cases for which the column rupture takes place before, at the droplet rear. The scaling of this quantity is important because it scales the volume of phase 1 which does not stay attached to the drop when it breaks but is rather sent back to the interface under the form of drops. In real continuous process conditions, this phenomenon can be limiting since it leads to the formation of an emulsion that keeps on growing on the interface as the drops are continuously crossing the interface.With the objective to characterize the driving force responsible for the film entrainment during drop crossing, two force ratios are defined. The first one, F * , compares the importance of forces that I for all cases). push on the interface, which are both the drop apparent weight (ρ 2 -ρ 1 )a c π d 3 /6 and the dynamic pressure due to drop inertia ρ 2 U 2 min πd 2 , to the stress which tends to pull back the fluid entrained of phase 1 towards the 1-3 interface and is defined by a gravity force at the scale of the drop volume (ρ 3 -ρ 1 )a c π d 3 /6, leading to

F * = (ρ 2 -ρ 1 )a c d + 6ρ 2 U 2 min (ρ 3 -ρ 1 )a c d = ξ 12 ξ 13 1 + 6 ρ 2 ρ 2 -ρ 1 U 2 min a c d . ( 10
)
The parameter F * is based upon U min , taken as the reference velocity in the film entrainment process in phase 3. It does not include any restoring surface tension force exerted by the film on the droplet, this force being assumed to be negligible once the crossing condition is achieved, which is fulfilled for all cases under consideration [i.e., ξ 12 ξ 13 f (Bo 13 )]. Once the drop (or particle) is accelerated in phase 3, the entrained film phase will unavoidably break due to the counterflow developing in the column, as evidenced by the velocity fields (see Fig. 9). One part is entrained in the wake of the rising drop and one part is pulled back towards the plane interface due to gravity, leading to the column thinning and rupture and then to the detachment of the encapsulated droplet. The second ratio represents the inertia excess with respect to the critical value for crossing in the static case, denoted by f (Bo 13 ):

I ex = F * / f (Bo 13 ). (11)
The normalized maximal column length L * max is plotted against I ex in Fig. 12. For both drops and solidlike droplets, L * max increases with I ex , and despite the scattering of the data (mainly due to the limitation of the numerical resolution at the instant of column breakup, leading to an average uncertainty of 5% of the computed values of L * max on the basis of simulations with mesh grids of 48 and 96 nodes per drop radius), two contrasting trends are observed even at low values of I ex , which are not related to drop deformation but to the Reynolds number Re col of the column flow (Re col is based upon phase 1 properties, the drop velocity at the instant of detachment and the column average thickness). For cases with small values of Re col (typically less than or equal to 1), the growth rate of L * max is about 6 times larger than that observed with larger Re col values; this behavior is observed for both deformable and solidlike droplets. Note that, concerning the branch with a lower growth rate of L * max , a different behavior between deformable droplets and solidlike drops appears at large I ex > 5: Whereas the particles seem to remain on the same linear evolution as that for lower I ex , the increase becomes much pronounced for the deformed droplets.

Even at low values of I ex , it is not clear how these two different regimes of film extension develop, despite the linear evolution observed between Re col and Re min (U min begin involved in I ex ). At low inertia, the observation of a large growth rate of L * max with I ex is obviously due to a delay of the pinch-off formation behind the drop, possibly resulting from the viscous resistance of the column to deformation that would lead to longer and thinner columns before breakup. Such an effect is expected to be scaled by the film Ohnesorge number Oh 13 (given in Table II): Indeed, Re col globally decreases with Oh 13 . However, in all cases studied, Oh 13 remains small (the maximum value of Oh 13 is 0.2), suggesting a limited influence of phase 1 viscosity on the deformation and breakup process of the column (in addition to the fact that, in most cases, λ 13 = 1). Moreover, the use of Oh 13 to find a unique scaling law for L * max was not successful, leading to the exhibition of these two different behaviors, which are both growing functions of the inertia excess I ex and which depend on the flow regime in the column characterized by Re col .

Volume of the encapsulating film

Whatever the column extension is, its breakage occurs at the rear of the droplet (even when column breakage first occurs at the bottom of the column, it is followed by a breakup event in the droplet wake region), leading to a volume V f of film entrained around the droplet. This encapsulating volume of phase 1 rapidly migrates towards the rear of the droplet due to buoyancy effects. The volume V * f , normalized by the droplet volume, has been determined for each case based on the procedure explained in Sec. II B.

For solidlike droplets, Fig. 13 reports the evolution of V * f as a function of ln(F * ). We first observe that film volume values are smaller than the particle or drop volume in the range of parameters investigated. These values collapse remarkably well on a single curve which can be fitted by a second-order polynomial. It seems that the evolution goes through a minimum around ln(F * ) 2, with a sharp increase rate at high inertia and a slight decrease rate at low inertia. For deformable droplets, a similar trend is observed when plotting V * f as a function of the logarithm of F * corrected by the viscosity ratio λ 12 to the power 0.4. Such a correction is related to the film drainage rate kinetics around the droplet during rising: The higher the λ 12 , the slower the film drainage and the larger the remaining film volume at the instant of detachment. However, this effect of λ 12 seems to be of importance for deformed droplets only. Indeed, one case of deformable droplet does not follow this trend but that of solidlike droplets (case D 6 , orange circles on the graph of Fig. 13). This drop has a high viscosity ratio (λ 12 = 50) and is weakly deformed at the time of column detachment (aspect ratio equal to 1.13), possibly explaining why this drop follows the trend of solidlike droplets. Therefore, the correction of F * by λ 0. 4 12 probably reflects a more subtle coupling FIG. 14. Profile of axial velocity (normalized by the instantaneous drop velocity) along the drop equator (z * = 0), for case S 9 , in the droplet frame. The covering film (phase 1) is still described by four mesh points at this instant.

between film drainage and drop deformation. For both solidlike droplets and deformable drops, the film volume encapsulating the drop seems to be controlled by inertia and not by surface tension forces. Of course, this is valid for particles or droplets crossing the interface, i.e., fulfilling the criterion ξ 12 /ξ 13 f (Bo 13 ), and is confirmed by the fact that V * f is correlated neither to the interface Bond number Bo 13 or Bo 13 ( ξ 12 ξ 13 -1) nor to a droplet Weber number (based on γ 13 and drop velocity at the minimum point or at the film detachment instant). At high inertia, increasing inertia (i.e., increasing F * or F * λ 0.4 12 ) favors the growth of the encapsulating volume. It is interesting to note that the two curves of Fig. 13 seem to merge in the limit of high inertia. At low inertia, but for crossing solidlike or deformable droplets, their respective behaviors diverge, the encapsulating film being significantly smaller for deformed drops than for solidlike drops, and the presence of a minimum is much less pronounced for deformable droplets than for solidlike droplets. This observation calls for deeper insight into the encapsulation process in this regime.

To that end, we have reported in Fig. 13 some screenshots of the phase indicator function fields, immediately after which the film breakage occurred for a few cases corresponding to low-and highinertia regimes, for deformable and solidlike droplets. The encapsulated volume can be visualized on these fields by the yellow ring around the blue droplet. Series S 9 , S 8 , and S 11 correspond to increasing values of F * for solidlike droplets in the low-inertia regime [ln(F * ) < 2]. For case S 9 , the film volume is clearly composed of two significant contributions, one due to the film coating the top of the particle and one at the rear of the particle. With S 8 , the contribution of this film at the top has decreased compared to the bottom one and for case S 11 , the volume at the rear of the drop has increased due to increasing inertia and keeps on increasing in cases S 3 and S 1 as F * is increased. Hence, the presence of a minimum of V * f can be due to the relative weight of the film remaining at the top of the drop at the instant of detachment in low-inertia conditions. This can be explained by the fact that the rate of film drainage coating the particle is an increasing function of inertia at the front of the droplet. Note that this effect is even emphasized in case S 9 due to the lower viscosity ratio λ 13 = 0.1 compared to other solidlike drop case series where λ 13 = 1.

To illustrate this effect of λ 13 on the film drainage, we have reported in Fig. 14 the axial velocity profile in the radial direction along the equator for case S 9 , at the instant of detachment. The continuity of tangential stresses at the interface imposes a strong gradient of velocity in phase 3 close to the interface (within the external boundary layer), making the velocity at the interface between the film and phase 3 to be quite small. This is a condition allowing a slow film drainage in addition to the quasinull velocity at the drop surface. Thus, a larger film thickness can be observed in cases where λ 13 < 1 (which is consistent with the results of Manga and Stone [START_REF] Manga | Low Reynolds number motion of bubbles, drops, and rigid spheres through fluid-fluid interfaces[END_REF] on the drainage of the column at low Reynolds number).

In the deformable drop cases (D 12 , D 13 , D 9 , and D 1 series), such a contribution to V * f of the film coating the top of the drop seems to be strongly attenuated at low inertia, the film being drained to the rear of the drop before detachment occurred (see cases D 12 and D 13 ). Compared to the solidlike cases, smaller viscosity ratios λ 12 (equal to 1 and 5, respectively) accelerate the film drainage, leading to smaller encapsulation film volumes with a systematic small contribution of the drop top film. Note that for both cases D 12 and D 13 , λ 13 = 0.1asforS 9 . However, deformation also seems to play a role in the film drainage at low inertia. More precisely, the transition from an oblate to a prolate shape when the drop rises from the interface to the point of minimum velocity is thought to accelerate the drainage of the film from the top to the rear of the drop. This effect induced by deformation can be observed for D 13 in Fig. 15 at Ar = 10, where the drop velocity falls to zero when crossing the interface, taking even small negative values due to the resistance force of the interface. Then the droplet becomes elongated (prolate shape) when trapped at the interface and even more elongated due to the strong acceleration suddenly experienced in phase 3 (of lower viscosity than phase 1), until point 5 in the plot, where it begins to flatten again (oblate shapes) during its rise in phase 3. Even though λ [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF] 1, the screenshots clearly suggest that such a prolate shape will induce a strong downward flow inside the drop, which favors a faster drainage. Prolate shapes are always obtained in such cases of deformable droplets crossing the interface at low Reynolds number, which correspond to the lower values of F * in Fig. 13. Note that, in the simulations of Shopov and Minev [START_REF] Shopov | The unsteady motion of a bubble or drop towards a liquid-liquid interface[END_REF] and in experiments of Bonhomme et al. [START_REF] Bonhomme | Inertial dynamics of air bubbles crossing a horizontal fluid-fluid interface[END_REF] in the case of bubbles trapped at a liquid-liquid interface, crossing regimes at low Re and Ar exhibited similar elongated shapes. On the contrary, at high values of F * (cases D 9 and D 1 in Fig. 13), with the Reynolds and drop Weber numbers We 12 being large, the drop is oblate before reaching the interface and during the crossing, as already shown in Fig. 7. Consequently, there is no acceleration of film drainage due to a change of shape and the evolution of film volume with inertia (corrected by λ 0.4 12 ) is consistent with the trend observed with solidlike drops.

Scaling to predict U min

Because the minimum velocity U min has been identified as a reference velocity for the scaling of film entrainment dynamics during crossing, the scaling of this quantity to independent parameters of the problem is required. To this purpose, first we examine the decrease in kinetic energy during crossing, which involves U min and the maximal drop velocity U max before reaching the interface. Then U max is scaled as a function of the dimensionless parameters known a priori.

We define the normalized ratio E * k of kinetic energy variation between the minimum and maximum velocity values, accounting for the added mass coefficient dependence with deformation in the case of deformable droplets

E * k = E kmax -E kmin E kmax , ( 12
)
where

E kmax = 1 2 (m 2 + C M max m 1 )U 2 max and E kmin = 1 2 (m 2 + C M min m 1 )U 2
min are the total kinetic energy at the instant of the maximal velocity and minimal velocity, respectively, C M max and C M min being the added mass coefficients at U max and U min for a displaced mass m 1 of phase 1 and the mass m 2 of the drop. For spherical droplets, C M = 0.5, whereas for deformable droplets, C M is a function of the instantaneous droplet aspect ratio [START_REF] Lamb | Hydrodynamics[END_REF]

C M = α 0 α 0 -2 , with α 0 = 2(ζ 2 0 + 1)[1 -ζ 0 cot -1 (ζ 0 )] for prolate shapes and α 0 = 2(ζ 2 0 + 1)ζ 0 cot -1 (ζ 0 ) -ζ 2 0
for oblate shapes, with ζ 0 = (χ 2 -1) -1/2 . In Fig. 16, the kinetic energy ratio E * k is plotted as a function of the excess of inertia I ex for both rigid and deformable droplets. It can be observed that all points gather on a single master curve, which can be fitted by a double exponential decay, which tends to zero as the excess of inertia tends to infinity. This means that, as the inertia of the droplet increases with respect to the minimum needed to cross the interface, the droplet will not be slowed down at the interface (U min ≈ U max ) and the ratio E * k tends to vanish. On the contrary, at low I ex , the crossing dynamics will tend FIG. 17. Archimedes number as a function of 3 2 Re max (1 + 0.15 Re 0.687 max ) for solidlike droplets. Orange squares correspond to all solidlike droplets in the three-phase flow simulations and the yellow square corresponds to the S 1 case in a two-phase flow simulation (i.e., without the interface) where the terminal velocity can be reached and fits with Schiller and Naumann's correlation (black dashed line). towards the static configuration (i.e., E * k = 1), the initial kinetic energy of the drop motion below the interface being fully dissipated.

The fit of the curve of Fig. 16 then provides an implicit relationship allowing us to predict U min for a given U max . Note that it has not been considered that the added mass coefficient can increase depending on the distance between the drop and the interface [START_REF] Pierson | Inertial settling of a sphere through an interface. Part 1. From sphere flotation to wake fragmentation[END_REF][START_REF] Milne-Thomson | Theoretical Hydrodynamics[END_REF]. However, the exact value of the added mass coefficient has only a small effect on the value of E * k here, because of the small density difference between the liquid phases.

Finally, the maximum velocity U max reached by the drop during its rising in phase 1 needs to be scaled, by assuming that it is the terminal velocity. For solidlike droplets, the force balance between buoyancy and drag sufficiently far from the interface leads to an implicit relation between Re max and Ar:

Re max 3 32 C D (Re max ) = Ar. (13)
If we substitute in (13) the Schiller and Naumann correlation [START_REF] Schiller | A drag coefficient correlation[END_REF] for the drag coefficient, we obtain Ar = 3 2 Re max 1 + 0.15 Re 0.687 max . (

) 14
Figure 17 plots the Archimedes number as a function of the right-hand-side term of Eq. ( 14). It is a linear fit which deviates from the first bisector by a nearly constant factor of 14%. Even if Schiller and Naumann's correlation [START_REF] Schiller | A drag coefficient correlation[END_REF] has some finite accuracy (a few percent), this deviation is more likely to be due to the fact that terminal velocity is not reached by solidlike droplets before interacting with the interface (the drops travel a distance of 3d in the simulations). A proof of that is given by case S 1 , which has a maximum velocity below the interface of 0.8 m/s. When running the simulations on longer travel distances without the plane interface (two-phase flow simulation), the terminal velocity found is 0.97 m/s and the corresponding Re max value falls on the first bisector of Fig. 17. For deformable drops, the plot of Re max as a function of Ar displays a linear evolution in the range of parameters investigated (with a slope equal to 2.25). In this case, results cannot be compared with the well-known law of the drag coefficient inserted in Eq. [START_REF] Aristoff | Water entry of small hydrophobic spheres[END_REF]. Moreover, because of viscous effects (from both internal and external phases), the response in deformation to the stress is always delayed and when the drop center of mass velocity reaches a plateau, the drop shape is not always steady and keeps on deforming, as already shown in Fig. 8. However, the linear fit of Fig. 18 suggests that the drag coefficient in Eq. ( 13) is nearly constant, as a consequence of two competing effects when increasing Archimedes number: Increasing Re max decreases C D but also increases drop deformation, which increases C D . So the plateau reached by the drop velocity far from the interface corresponds to a regime where this balance between deformation and drag is reached and seems to be independent of the instantaneous value of the deformation, which strongly increases during the drop rising until the interface. It is interesting to note that this linear behavior does not seem to be influenced by viscosity ratios λ 12 or λ 13 and covers a wide range of drop deformation (aspect ratio varying between 1.2 and 2.5).

To summarize, U max can be scaled with Archimedes number Ar for both solidlike (Fig. 17) and deformable droplets (Fig. 18) and U min can be deduced from U max through a known function of I ex (Fig. 16). This scaling gives access to the prediction of the different quantities relevant for the encapsulating process in a wide range of flow parameters in the inertial regime.

IV. CONCLUDING REMARKS AND PROSPECTS

In this work, the study of a drop crossing a liquid-liquid interface under the action of a centrifugal field has been undertaken by means of resolved numerical simulations. The conditions for crossing investigated correspond to an inertial regime in phase 1 (10 < Re max < 200). Two types of drops have been studied: solidlike droplets with a high internal viscosity and surface tension, mimicking rigid particles, and deformable droplets. The level-set method used in the simulations was first validated in two-phase flow configurations (i.e., without the plane interface) by comparing the computed terminal velocity of undeformed drops at vanishing Weber number and the drop aspect ratio at higher Weber numbers with existing results in the literature. As three-phase flow validation, the crossing criterion of an interface by a solidlike droplet in static conditions (i.e., starting from the interface without initial velocity) was successfully compared to a theoretical model based upon the balance between gravity and surface tension forces. A large number of simulations were run in crossing conditions for both solidlike and deformable drops, covering a wide range of the seven nondimensional parameters describing this problem. The main results of this study can be summarized as follows.

(i) The static criterion derived for a sphere to cross an interface is still a relevant reference of the minimum inertia required for crossing, and this conclusion stands for both solidlike and deformable droplets. In dynamic conditions, i.e., under the influence of the drop impact velocity (close to its terminal velocity), the interface crossing criterion remains quite similar to the static criterion, because the effect of inertia (by neglecting the velocity decrease for the drop arriving at the interface) is to double the drop apparent weight in the best case. This can be inferred from the fact that with a two-liquid phase system of close density and in the inertial regime, the dynamic pressure force exerted by the drop on the interface is of the same order as its apparent weight. This results in an effective equivalent Bond number (Bo 13eq = 1.59 Bo 13 ) to be compared to the minimal value required for crossing at a given ξ 12 /ξ 13 for the system, given by the theoretical criterion in static conditions (then inertia can only slightly shift this criterion, which remains a valuable reference).

(ii) In dynamic crossing conditions, the drop velocity during crossing goes through a minimum velocity U min which is used as a reference velocity to scale inertia additional to the drop apparent weight.

(iii) During crossing, film entrainment of the lighter phase by the drop leads to the formation of a column that extends and thins as it is pulled by the drop in phase 3, which ends to break behind the drop, leading to the formation of a volume encapsulating the drop. The main mechanism responsible for the rupture of the entrained column is the competition between the inertial force that makes the drop rise in phase 3 and the gravity force acting on the lighter column phase that pulls it back in the opposite direction (towards the interface). In order to characterize the film entrainment and the encapsulating volume, these observations lead us to define two nondimensional parameters: The first is based on the ratio of these two forces (denoted by F * ) and the second corresponds to the first one rescaled by the criterion of crossing in static conditions (named inertia excess and denoted by I ex ).

(iv) The maximum column length at breakup increases with inertia excess I ex , the evolution being linear at low I ex but following two distinct slopes that depend upon the Reynolds number in the film column, independently of deformation. A higher growth rate corresponds to viscous flows in the column. The transition between these two regimes has not yet been elucidated.

(v) For solidlike droplets, the encapsulating film volume V * f is remarkably well described by the force ratio F * , exhibiting a minimum for increasing F * , which originates from the contribution of the film volume on the top part of the particle which is non-negligible at low inertia. For deformable droplets, V * f is also very well described by F * provided it is corrected by the viscosity ratio λ 12 between the drop and the film phase to the power 0.4, because λ 12 has a strong impact on the tangential velocity along the drop surface, then on the film drainage rate (in the range of parameters investigated with 0.05 λ 13 1; note that the viscosity ratio λ 12 has a greater influence on V * f than λ 13 ). At low inertia, the encapsulating volume is much smaller for deformable drops than for solidlike drops. This is promoted by two effects: First is the lower λ 12 and second is the oblate-prolate shape transition observed during the crossing between the interface level and the location of U min , both effects tending to accelerate the film drainage during the drop rising in phase 3, leading therefore to minimization of the film volume attached to the drop.

(vi) As the main inertia parameters (F * and I ex ) exhibited in this study are defined based on the minimal velocity U min , scaling laws have been proposed to predict the latter as a function of independent parameters of the problem. First the maximum impact velocity U max can be computed from Ar by using Fig. 17 or 18 and then the relationship of the decrease of kinetic energy during crossing as a function of I ex , displayed in Fig. 16, involves these two velocities and can be used to implicitly predict U min .

Together the results make possible the prediction of the critical condition and dynamics of crossing, maximum column length entrained, and encapsulating volume for solidlike and deformable droplets in a rather large range of flow parameters.

Future work based on this study should address the comparison of the numerical results with experimental data. Then, with the objective of developing an efficient encapsulation process based on interface crossing, surfactants are necessary to stabilize the encapsulated droplet in order to make it resistant to breakup while the drop is rising in phase 3. An analysis of the influence of surfactants, adsorbed at the liquid-liquid plane interface, on the crossing dynamics and encapsulation volume seems therefore to be highly relevant to continue this work.

FIG. 2 .

 2 FIG. 2. Reynolds-number evolution as a function of normalized distanced covered by the droplet z * = z d for different mesh resolutions for (a) case 1, Re T = 20, We 12 = 0.11, and simulation domain 6R × 12R and (b) case 2, Re T = 180, We 12 = 0.31, and simulation domain 8R × 16R.

FIG. 5 .

 5 FIG. 5. Diagram showing the crossing and no-crossing zones as a function of static theory predictions (solid lines), Eqs. (8)and (9). The squares and circles correspond to simulations of solidlike droplets and deformable droplets, respectively, in the dynamic case, with the points labeled by the Archimedes number Ar. Symbols are closed in the case of crossing in the simulation and open otherwise.

FIG. 6 .

 6 FIG. 6. Evolution of the Reynolds number Re = ρ 1 U d μ 1 (U is the instantaneous drop velocity) as a function of z * (drop center of mass position normalized by R), for cases S 1 and S 2 , along with screenshots of the phase indicator function field issued from the simulations. Here z * = 0 is the interface position. Figures appearing on the curves of the Reynolds number time signal are related to the image sequences; dimensionless time is indicated below the screenshots, using the characteristic timescale d ω 2 r i .

FIG. 7 .

 7 FIG. 7. Evolution of the Reynolds number Re and aspect ratio χ as a function of z * = z/R for cases D 1 and D 2 . Dashed lines correspond to the two-phase flow simulations (without the interface), used as a comparison basis for the three-phase flow case.

FIG. 8 .

 8 FIG. 8. Evolution of the Reynolds number as a function of z * = z/R for two cases of deformable droplets D 1 and D 2 and the case S 1 of a solidlike droplet at the same Ar and Bo 13 , along with screenshots of the simulations. The points indicated on the Reynolds number curve are relative to the images shown; dimensionless time is indicated below the screenshots, using the characteristic timescale d ω 2 r i .

FIG. 9 .

 9 FIG.9. From left to right, velocity field of cases S 1 , D 1 ,and D 2 right before the column detachment. The scaling velocity vector is indicated on the rightmost image. Here x * and z * are the coordinates of the calculation domain normalized by the radius of the droplet. In all cases, the viscosity of the two continuous phases is the same (λ 13 = 1).

FIG. 10 .

 10 FIG. 10. Illustration of the crossing time criterion.

FIG. 11 .

 11 FIG. 11. Crossing time t * cr as a function of 2d/(U max + U min ) for both solidlike and deformable droplets (dimensional physical parameters are given in Table I for all cases).

FIG. 12 .

 12 FIG.12. Mapping of nondimensional maximal column length L * max = L max /d as a function of the nondimensional inertia excess I ex for both nondeformable and deformable droplets (see Fig.6for the definition of L max ). According to Re col , two regimes of increase of L * max are observed. Based on mesh convergence tests, the uncertainty of the reported values of L * max after column breakup is about 5%.

FIG. 13 .

 13 FIG. 13. Dimensionless film volume V * f = V f /V drop as a function of the parameter F * in the case of solidlike droplets and as a function of F * λ 0.4 12 in the case of deformable droplets. Error bars have been computed according to the criterion discussed in Sec. II B. A series of screenshots of the phase indicator function at the instant of column detachment is presented for some cases.

FIG. 15 .

 15 FIG. 15. Evolution of the Reynolds number Re = ρ 1 U d μ 1 as a function of the centroid position for case D 13 . The droplet deforms into a prolate shape when decelerating and then accelerating in phase 3 (λ 13 = 0.1). Screenshots correspond to the different positions indicated on the red curve with corresponding values of the aspect ratio.

FIG. 16 .

 16 FIG. 16. Kinetic energy variation ratio E *k as a function of inertia excess I ex . Orange squares denote solidlike droplets, blue circles denote deformable droplets, and the dashed line is a double exponential fit of all points.

FIG. 18 .

 18 FIG.[START_REF] Dietrich | Passage of a bubble through a liquid-liquid interface[END_REF]. Maximum Reynolds number as a function of Archimedes number for deformable droplets. The dashed line represents a linear fit.

TABLE I .

 I Physical parameters of the simulations (Sim.) presented in this paper for solidlike and deformable droplets. Density is in kg/m 3 , dynamic viscosity in Pa s, and surface tension in N/m. The droplet size is d = 200 μm for all simulations.

	Sim.	ρ 1	ρ 2	ρ 3	μ 1	μ 2	μ 3	γ 12	γ 13	a c
						Solidlike
	S 1	995	1100	1000	0.001	0.05	0.001	1	0.01	2740g
	S 2	995	1100	1000	0.002	0.1	0.002	1	0.01	1551g
	S 3	995	1100	1000	0.001	0.05	0.001	1	0.01	1551g
	S 4	995	1100	1000	0.001	0.05	0.001	1	0.01	2660g
	S 5	995	1100	1000	0.001	0.05	0.001	1	0.01	3262g
	S 6	995	1100	1000	0.0015	0.075	0.0015	1	0.01	1551g
	S 7	995	1100	1000	0.0007	0.035	0.0007	1	0.01	1551g
	S 8	930	1080	1000	0.002	0.1	0.002	1	0.0045	1359g
	S 9	850	1100	1000	0.01	0.5	0.001	2	0.02	5097g
	S 10	980	1100	1000	0.001	0.05	0.001	1	0.01	2548g
	S 11	980	1100	1000	0.001	0.05	0.001	1	0.01	1733g
						Deformable
	D 1	995	1100	1000	0.001	0.05	0.001	0.01	0.01	2740g
	D 2	995	1100	1000	0.001	0.005	0.001	0.01	0.01	2740g
	D 3	995	1100	1000	0.001	0.02	0.001	0.01	0.01	2740g
	D 4	995	1100	1000	0.01	0.05	0.01	0.01	0.01	2740g
	D 5	995	1100	1000	0.01	0.05	0.01	0.0134	0.0134	2740g
	D 6	995	1100	1000	0.001	0.05	0.001	0.01	0.01	1529g
	D 7	960	1200	1000	0.0017	0.0085	0.0017	0.01	0.01	1359g
	D 8	960	1200	1000	0.0007	0.0035	0.0007	0.01	0.01	1359g
	D 9	960	1200	1000	0.0017	0.085	0.0017	0.01	0.01	1359g
	D 10	960	1200	1000	0.0007	0.035	0.0007	0.01	0.01	1359g
	D 11	980	1200	1000	0.001	0.0005	0.001	0.006	0.006	503g
	D 12	850	1100	1000	0.01	0.01	0.001	0.02	0.02	4077g
	D 13	850	1100	1000	0.01	0.05	0.001	0.02	0.02	5097g
	D 14	913	1081	997	0.005	0.0024	0.001	0.036	0.011	898g
	D 15	930	1081	997	0.01	0.0024	0.001	0.036	0.011	923g
	D 16	950	1081	997	0.02	0.0024	0.001	0.036	0.011	926g

TABLE II .

 II Nondimensional numbers of all simulations. Here Re max , Re min , and We max are based on the properties of phase 1. Note that drop breakup occurs during the interface crossing for simulations D 12 and D 14 . 12 λ 13 Bo 12 Bo 13 Ar Oh 12 Oh 13 Re max Re min We max I ex χ max χ min χ det

	Sim. ξ 12	ξ 13	ξ 12 ξ 13
	S 1 0.105 0.005 21 50 1 0.03 0.13 53 0.1 0.02 160 134 0.13	11	1	1	1
	S 2 0.105 0.005 21 50 1 0.02 0.076 20 0.21 0.04 48	16 0.046 1.56 1	1	1
	S 3 0.105 0.005 21 50 1 0.02 0.076 40 0.1 0.02 114 62 0.065 3.05 1	1	1
	S 4 0.105 0.005 21 50 1 0.03 0.13 52 0.1 0.02 158 130 0.125 10.5	1	1	1
	S 5 0.105 0.005 21 50 1 0.03 0.16 57 0.1 0.02 178 158 0.16 14.8	1	1	1
	S 6 0.105 0.005 21 50 1 0.02 0.076 27 0.16 0.03 69	25 0.054 1.76 1	1	1
	S 7 0.105 0.005 21 50 1 0.02 0.076 57 0.07 0.01 177 104 0.077 3.82 1	1	1
	S 8 0.161 0.075 2.14 50 1 0.02 2.07 22 0.22 0.07 52	17 0.058 1.61 1	1	1
	S 9 0.294 0.176 1.67 50 0.1 0.06 3.75 10 0.75 0.17 19	0 0.106 1.05 1	1	1
	S 10 0.12 0.02	6 50 1 0.03 0.5 54 0.1 0.02 94	85 0.045 5.07 1	1	1
	S 11 0.12 0.02	6 50 1 0.02 0.34 45 0.1 0.02 72	36 0.026 1.83 1	1	1
	D 1 0.105 0.005 21 50 1 2.82 0.134 53 1.06 0.02 124 95 7.72	6.4 1.77 1.93 1.63
	D 2 0.105 0.005 21 5	1 2.82 0.134 53 0.1 0.02 112 92 6.37	6.14 2.39 3.24 1.77
	D 3 0.105 0.005 21 20 1 2.82 0.134 53 0.43 0.03 116 95 6.8	6.38 2.21 2.12 1.808
	D 4 0.105 0.005 21 5	1 2.82 0.134 5 1.06 0.22 7	2.8 2.37	2.17 1.24 1.09 1.16
	D 5 0.105 0.005 21 5	1 2.1 0.1	5 0.92 0.19 7	1.3 1.77	1.4 1.21 0.998 1.11
	D 6 0.105 0.005 21 50 1 1.57 0.076 40 1.06 0.02 96	40 4.63	1.85 1.33 1.4 1.13
	D 7 0.25 0.042 6	5	1 2.43 0.4 28 0.17 0.03 60	38 5.4	3.67 2.1 0.95 1.26
	D 8 0.25 0.042 6	5	1 2.43 0.4 69 0.07 0.01 151 116 5.8	5 2.52 1.79 1.5
	D 9 0.25 0.042 6 50 1 2.43 0.4 28 1.74 0.03 64	30 6.2	2.83 1.4 1.16 1.16
	D 10 0.25 0.042 6 50 1 2.43 0.4 69 0.71 0.01 162 111 6.76	4.68 1.95 1.53 1.39
	D 11 0.224 0.02 11 0.5	1 1.8 0.16 33 0.01 0.02 77	58	5	4.37 2.23 1.48 1.42
	D 12 0.294 0.176 1.67 1 0.1	5	3	9 0.15 0.17 16	1.8 7.3	1 1.65 0.6 1.17
	D 13 0.294 0.176 1.67 5 0.1 6.25 3.75 10 0.75 0.17 17 -2.5 8.3	1.05 1.83 0.72 1.1
	D 14 0.184 0.092 2 0.48 0.2 1.37 2.25 18 0.02 0.08 48	37 4.76	4.7 1.82 1.39 1.31
	D 15 0.162 0.072 2.25 0.24 0.1 1.44 2.09 10 0.02 0.16 22	17 3.5	3.96 1.43 1.13 1.33
	D 16 0.138 0.049 2.78 0.12 0.05 1.36 1.6	5 0.02 0.3	8	6.5 1.97	5.14 1.21 0.92 1.76

λ

ACKNOWLEDGMENTS

The authors gratefully acknowledge financial support received from the French Research National Agency (ANR) through the program "Encapsulation Process by Droplet Interface Crossing" (project reference: ANR-15-CE08-0026).