Mass transfer from a spherical capsule in absence of flow

Résumé

Le transfert de masse à partir de particules présente de nombreuses applications biomédicales, telles que le relargage contrôlé de médicaments, la culture cellulaire, ou encore la conception d'organes bioartificiels. En conditions opérationnelles, les particules sont soumises à des écoulements de fluides de différentes natures, dont l'effet sur le transport du soluté est aujourd'hui encore mal connu. Dans cette thèse, nous nous intéressons à l'effet de l'écoulement sur le transfert de masse à partir de capsules et de fibres coeur-coque.

Le cas du relargage d'un soluté par une fibre coeur-coque confinée dans un canal plan et soumise à un écoulement de Poiseuille est étudié à partir de simulations numériques en deux dimensions basées sur la méthode de Boltzmann sur réseau. Les effets combinés de l'écoulement et de la présence de la coque sur le transfert de masse sont analysés sur une large plage de nombre de Reynolds couvrant des écoulements stationnaires et instationnaires. Une corrélation donnant le nombre de Sherwood (le coefficient de transfert de masse adimensionnel) en fonction de la perméabilité de la coque, du nombre de Reynolds, et du nombre de Schmidt est proposée.

La perméabilité de la coque est une propriété difficilement mesurable. Les méthodes actuelles de caractérisation ne permettent d'obtenir qu'une perméabilité effective pour l'ensemble de la particule (coeur et coque réunis), et non spécifiquement pour la coque. Une nouvelle méthode de caractérisation est proposée et validée sur des données numériques et expérimentales. Cette méthode consiste à extraire la perméabilité de la capsule par analyse inverse, en ajustant les courbes de relargage avec des solutions de la seconde loi de Fick calculées en une dimension par la méthode des différences finies.

Le relargage d'un soluté par une capsule placée dans un écoulement de cisaillement est également analysé par des simulations tridimensionnelles couplant la méthode de Boltzmann sur réseau et la méthode des frontières immergées pour les interactions fluidestructure. L'écoulement améliore significativement l'efficacité du transfert de masse par convection forcée. L'effet des conditions aux limites à la surface de la capsule est examiné. Considérer des conditions aux limites de Dirichlet, comme fait classiquement dans la littérature, conduit à des nombres de Sherwood bien plus élevés que dans le cas de conditions aux limites de continuité, plus adaptées pour modéliser le relargage d'un soluté.

La dynamique d'une suspension de capsules déformables en écoulement dans une constriction est aussi étudiée. Un diagramme d'états décrivant la transition entre le passage et le blocage de la constriction par les capsules est déterminé en fonction de leurs propriétés géométriques et mécaniques. Il peut être utilisé, par exemple, afin d'optimiser le design de systèmes microfluidiques dans le but d'améliorer l'efficacité du transfert de masse.

Mots clés : Transfert de masse, écoulement, convection forcée, capsules, fibres, membranes, perméabilité, relargage de médicament

Abstract

Mass transfer from particles is encountered in many biomedical applications, such as controlled drug delivery, cell culture, or in designing bioartificial organs. In operating conditions, particles are subjected to various flows, whose effect on solute transport is still not well understood and controlled. In this PhD thesis, we study the effect of the flow on mass transfer from core-shell capsules and fibers, using numerical simulations.

Solute release from a core-shell fiber confined in a channel and subjected to Poiseuille flow is studied using two-dimensional lattice Boltzmann simulations. The combined effects of the flow and the shell permeability on mass transfer are analyzed over a wide range of Reynolds number covering both steady and unsteady flows. A new correlation giving the Sherwood number (the dimensionless mass transfer coefficient) as a function of the shell permeability, the Reynolds number, and the Schmidt number is proposed.

The shell permeability is difficult to measure. Current characterization methods only allow the determination of an effective permeability for the whole particle (core and shell combined), and not specifically the shell permeability. A novel method to characterize capsules shell permeability is proposed and validated on both numerical and experimental data. This method consists in extracting the capsule permeability by fitting the release curves with numerical solutions of Fick's second law of diffusion computed in one dimension using the finite difference method, by assuming spherical symmetry.

Solute release from a capsule placed in a shear flow is also analyzed by the mean of three-dimensional simulations coupling the lattice Boltzmann method and the immersed boundary method for the fluid-structure interaction. The flow is found to enhance the mass transfer efficiency by forced convection. The effect of the boundary conditions set at the surface of the capsule is examined, as well. Considering Dirichlet boundary conditions at the particle surface, as classically done in literature, leads to significantly higher Sherwood numbers as compared to the case of continuity of both the concentration and the mass flux, which is more adapted to model solute release.

The suspension dynamics of soft capsules in a microfluidic constriction is also studied. A state diagram describing the transition from capsule passage to blockage is determined as a function of their geometrical and mechanical properties. It can be used, for example, to optimize the design of microfluidic devices in order to enhance the mass transfer efficiency. Keywords: Mass transfer, flow, forced convection, capsules, core-shell fibers, membranes, solute permeability, drug release Mass transfer is of fundamental and practical interest. It is encountered in many natural phenomena and industrial processes, spanning from transport industry to bioengineering. This PhD thesis focuses on biomedical applications of mass transfer from particles, such as spherical capsules and cylindrical core-shell fibers, under flow conditions. These particles are composed of an inner core, which is initially loaded with a solute (e.g. drug molecules or a contrast agent), and a coating polymeric thick shell or thin membrane. The role of the latter is, among others, to protect the encapsulated material from external mechanical and chemical damages. They have diameter ranging from few millimeters down to few nanometers depending on their applications. Figure 1 shows sets of spherical capsules and cylindrical fibers, where the core-shell structure of the particles is well visible. Capsules and fibers have known a growing interest in the last decades and are increasingly used in bioengineering.

Controlled drug delivery

Capsules and fibers are very promising for controlled drug delivery applications, where they are employed as carrier particles to deliver their encapsulated material, for example insulin to treat diabetes, at targeted sites and at desired rates. Drug carrier particles can be administrated either orally or by intravenous injection, or they can be implanted locally, for instance, near a tumor (see Fig. 2a). Depending on the way of administration, these devices may be faced with various physical obstacles, such as dense tissues, limiting drug penetration. They may also be subjected to chemical deteriorations, as those induced by extreme variations in pH in the gastrointestinal tract due to the food intake [START_REF] Tibbitt | Emerging frontiers in drug delivery[END_REF]. The main challenge in controlled drug delivery is to design systems that overcome these physical and chemical barriers, and release the encapsulated drug in a controlled manner. A poorly designed device may release drug rapidly causing the concentration to exceed the toxic limit (see Fig. 2b), and potentially leading to overdoses. Shortly after, the drug is being absorbed until reaching a concentration below the therapeutic level. This reduces the efficacy of the treatment and compels patients to repeat drug intakes regularly, which can be very constraining and uncomfortable. In contrast, a properly designed system should enable the control of the release rate in order to maintain the drug concentration within Figure 1: (a) Liquid-core alginate capsules with diameter of few millimeters. Image extracted from Ben Messaoud et al. [2]. (b) Micrometric core/shell poly(ethylene oxide)/Eudragit fibers used for site-specific drug release. Image taken from Jia et al. [3].

the therapeutic window for a predetermined period of time, from days to years. Capsules and fibers are excellent candidates for this purpose thanks to their shell (or membrane) that can be tailored to obtain the desired release kinetics.

Organs-on-a-chip

Fibers and capsules are also widely used in artificial systems, such as organs-on-a-chip devices, where their role is to moderate mass exchange between their encapsulated solute or cells and the surrounding flowing solvent, thanks to their membrane. These miniaturized systems are designed on microfluidic chips and aim at artificially reproducing the biological functions of healthy organs and tissues [4,5]. The bioartificial pancreas developed at BMBI laboratory (see Fig. 3) exploits this technique to mimic the biological function of a healthy pancreas [4]. To do so, islets of Langerhans are entrapped into microwells and secrete insulin as a response to glucose stimulations, in order to regulate the glucose level in the surrounding fluid. It was found that insulin secretion is enhanced when islets of Langerhans are cultivated in microfluidic biochips rather than when they are simply cultivated in Petri dishes. Furthermore, their viability is enhanced thanks to the surrounding flow that delivers all the necessary nutrients and oxygen to them. Blood oxygenators Blood oxygenators are another example of bioartificial systems developed to supplement biological organs, here, the lungs. Their role is to saturate the venous flow with oxygen, and remove carbon dioxide from it [6], as naturally done by healthy lungs. They are used, for instance, during open-heart surgery to replace the functions of the lungs for a short time (few hours), or to treat respiratory insufficiency for longer periods (few weeks to months). Most currently commercially-available blood oxygenators are extracorporeal and are based on fiber technology [7]. They are basically constituted of arrays of core-shell cylindrical fibers that are highly permeable to oxygen and carbon dioxide. Oxygen flows into the lumen of the fibers, and blood is pumped and flows transversely in the space outside the fiber bundles (see Fig. 4). Oxygen diffuses through the microporous walls of the fibers and is transferred into the bloodstream, while, oppositely, carbon dioxide is absorbed and captured by the fibers. Gas exchanges take place at the fibers surface and are controlled by their shell permeability. In all the above-mentioned applications, particles are systematically subjected to the flow of their surrounding fluid, which can be either a liquid (e.g. blood or water) or a gas (e.g. air), depending on the problem's nature. For example, drug carrier particles may be subjected to blood and air flows in the cardio-respiratory system, or to the flow of the interstitial fluid that surrounds the tissues. Particles in artificial devices undergo the flow of their surrounding solvent, which, for instance, can be triggered by syringe or peristaltic pumps.

The effects of the flow on mass transfer from core-shell particles are still not well characterized. On the one hand, theoretical models on mass transfer from core-shell particles [START_REF] Siepmann | Modeling of diffusion controlled drug delivery[END_REF][START_REF] Arifin | Mathematical modeling and simulation of drug release from microspheres: Implications to drug delivery systems[END_REF] are restricted to static conditions for simplicity, i.e. they assume that the surrounding fluid is at rest. On the other hand, most studies considering mass transfer under flow conditions are limited to uncoated particles [START_REF] Frankel | Heat and mass transfer from small spheres and cylinders freely suspended in shear flow[END_REF][START_REF] Acrivos | Heat transfer at high Péclet number from a small sphere freely rotating in a simple shear field[END_REF][START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF][START_REF] Morgan | The overall convective heat transfer from smooth circular cylinders[END_REF][START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Batchelor | Mass transfer from a particle suspended in fluid with a steady linear ambient velocity distribution[END_REF]. The effect of having a permeable coating shell on mass transfer under flow conditions has not been reported yet in the literature. Furthermore, these studies assume either constant concentration or constant mass flux at the particle surface. This assumption does not hold for applications such as drug delivery, for which the concentration within the particle is expected to decay over time as the encapsulated solute is released into the surrounding medium. There is, thus, need to better understand the effects of the flow on mass transfer from core-shell particles in order to optimize the design of future drug carriers and microfluidic devices to improve their performance. Trial-and-error experimental studies are generally long and expensive. Their costs can be significantly reduced by mathematical modeling and numerical simulations, where the design of the studied device can be adjusted easily, without any additional manufacturing cost and time.

The main objective of this PhD thesis is to numerically study the interplay between the solvent flow and the shell on mass transfer, while considering continuous boundary conditions at the particles surface. Several open questions will be addressed in this thesis, including:

1. How does the transition from steady to unsteady flow alter mass transfer?

2.
What is the effect of the coating shell on the release kinetics?

3. How to characterize the shell permeability to solute?

In what extent do the boundary conditions set at the particles surface affect the mass transfer efficiency?

To answer these questions, numerical simulations will be performed in both two and three dimensions, while varying the flow strength and the shell solute permeability. Most simulations are carried out using the lattice Boltzmann method, which allows to deal with complex geometries, offers large flexibility for the boundary conditions, and is highly parallelizable [START_REF] Succi | The lattice Boltzmann equation for fluid dynamics and beyond[END_REF][START_REF] Sukop | Lattice Boltzmann modeling[END_REF][START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF].

Several flow conditions are considered: Poiseuille flow, shear flow and flow in a microfluidic constriction. The aim of this PhD thesis is not to focus on a specific biomedical application of mass transfer from particles, with a specific encapsulated solute and surrounding solvent. On the contrary, the aim is to provide a better fundamental understanding of the underlying mechanisms that govern mass transfer from core-shell particles, considering several flow conditions that are commonly encountered in various biomedical applications. The emphasis is put on the impact of the flow and the contribution of the membrane permeability on the mass transfer efficiency. This PhD thesis is organized as follows. The first chapter is dedicated to mass transfer from a core-shell capsule in the absence of flow, in order to analyze solely the effect of the shell membrane permeability to solute, independently of the contribution of an external applied flow. The contribution of the flow to mass transfer is examined in the second chapter, in two dimensions first, for the case of a stationary core-shell cylindrical reservoir in cross-flow. The third chapter focuses on solute release from a spherical capsule freely rotating in shear flow, in three dimensions, while considering various boundary conditions at the particle surface. The fourth chapter deals with the collective motion of multiple soft capsules in a microfluidic constriction, which is exploited, for instance, to mechanically trigger solute release. Here, the emphasis is put on the clogging events that may dramatically alter the mass transfer efficiency.

Main contributions of this PhD thesis

The main contributions of this PhD thesis are listed in what follows, as a summary of the work achieved during these past three years.

• First contribution:

The first contribution concerns mass transfer from a core-shell cylindrical reservoir in cross-flow. It is known that mass (or heat) transfer is enhanced by forced convection, and this enhancement is quantified by various correlations for the Sherwood (or Nusselt) number [START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF][START_REF] Morgan | The overall convective heat transfer from smooth circular cylinders[END_REF][START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF][START_REF] Yousefi | Regimes of heat transfer in finitesize particle suspensions[END_REF]. However, no study accounting for the contribution of a coating shell could be found in the literature, despite core-shell particles being widely used in biomedical applications. It is found that decreasing the shell permeability shrinks the concentration boundary layer and, as a consequence, increases the Sherwood number. A new correlation for the Sherwood number combining both the contribution of the flow and the shell permeability is proposed. These original results were published as a regular article in Physical Review Fluids (Ref. [START_REF] Bielinski | Unsteady mass transfer from a core-shell cylinder in crossflow[END_REF]).

• Second contribution:

The shell permeability is a key parameter in designing capsules or fibers. However, this is an intrinsic property of the capsule and it is unknown prior to experiments. Current characterization methods [2,[START_REF] Dembczynski | Characterization of small molecules diffusion in hydrogel-membrane liquid-core capsules[END_REF][START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF][START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF] only allow the estimation of an effective diffusion coefficient for the whole capsule, core and membrane combined. Here, an efficient and accurate method is proposed to characterize diffusion coefficients in both the core and the membrane of capsules. The present method was successfully tested on experimental data taken from literature [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF][START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF], for the case of glucose release and absorption by capsules. The method also recovered accurately the estimated effective diffusion coefficients of these experimental studies.

• Third contribution:

The equations that are classically used to model solute release from controlled drug delivery devices are found to be not well adapted to core-shell particles, whose release kinetics is much more complex than their homogeneous counterparts. A practical equation modeling accurately the release kinetics of capsules as a function of their physicochemical properties is proposed. This equation is user-friendly and can be used in designing future core-shell capsules for pharmaceutical applications.

• Fourth contribution:

The fourth contribution of this PhD thesis concerns the boundary conditions set at the surface of the particle. In most studies [START_REF] Acrivos | Heat transfer at high Péclet number from a small sphere freely rotating in a simple shear field[END_REF][START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Batchelor | Mass transfer from a particle suspended in fluid with a steady linear ambient velocity distribution[END_REF][START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF], the particle is assumed to be sustained at constant concentration to simplify the problem and to limit the computations only to the outer fluid domain. However, this assumption is not appropriate to model solute release, during which the concentration within the particle decays over time. Here, we consider continuity of both the concentration and the mass flux at the surface of the particles, which naturally emerges in absence of interfacial resistance. It is found that imposing Dirichlet boundary conditions leads to higher Sherwood numbers, as opposed to when considering continuous boundary conditions, due to the increased concentration gradient at the particle surface.

• Fifth contribution: Many biomedical applications involve the flow of multiple particles in microfluidic constrictions [START_REF] Mcfaul | Cell separation based on size and deformability using microfluidic funnel ratchets[END_REF][START_REF] Sharei | A vector-free microfluidic platform for intracellular delivery[END_REF][START_REF] Myrand-Lapierre | Multiplexed fluidic plunger mechanism for the measurement of red blood cell deformability[END_REF][START_REF] Goff | Squeezing bio-capsules into a constriction: deformation till break-up[END_REF][START_REF] Luo | Solute release from an elastic capsule flowing through a microfluidic channel constriction[END_REF]. In such geometries, particles may clog the constriction entrance depending on their physical properties, and that may dramatically hinder the device efficacy. The clogging phenomenon is well characterized for rigid particles (see e.g. Ref. [START_REF] Marin | Clogging in constricted suspension flows[END_REF]). However, only few works deal with suspensions of soft particles in microfluidic constrictions. Here, the suspension dynamics of soft capsules flowing into a microfluidic constriction is studied using fully three-dimensional numerical simulations. The effects of both the capsules deformability and the constriction width on the formation of clogs are carefully analyzed. The results are reported in a state diagram indicating whether or not the particles clog the channel depending on their deformability and the neck-to-particle size ratio.

Introduction

The last decades have known an increasing interest in using capsules to encapsulate active agents, which has led to a huge amount of publications and patents. Capsules are composed of an inner liquid core, where active molecules are encapsulated, and an outer protective thick shell or thin membrane. They are fabricated through diverse techniques, including emulsification in flow-focusing microfluidic devices [START_REF] Seiffert | Microfluidic fabrication of smart microgels from macromolecular precursors[END_REF][START_REF] Chu | Fabrication and in situ characterization of microcapsules in a microfluidic system[END_REF][START_REF] Boskovic | Synthesis of polymer particles and capsules employing microfluidic techniques[END_REF][START_REF] Galle | On-chip encapsulation via chaotic mixing[END_REF][START_REF] Li | Microfluidic fabrication of microparticles for biomedical applications[END_REF] (see Fig. 1.1a). Capsules with more complex structures are being developed by adding multiple extra layers or chemically reinforced shells to meet specific performance requirements [START_REF] Dähne | Fabrication of micro reaction cages with tailored properties[END_REF][START_REF] Qiu | Studies on the drug release properties of polysaccharide multilayers encapsulated ibuprofen microparticles[END_REF][START_REF] Antipov | Sustained release properties of polyelectrolyte multilayer capsules[END_REF][START_REF] Zarket | Onion-like multilayered polymer capsules synthesized by a bioinspired inside-out technique[END_REF], as illustrated in Fig. 1.1b. The main advantage of capsules over their monolithic counterparts is their enhanced ability to monitor the release of their cargo in a controlled manner. Moreover, their membrane improves the mechanical robustness, and protects the encapsulated material against undesired external chemical reactions and mechanical damages. Capsules are also increasingly used as miniaturized bioreactors, where living cells are enclosed to carry biochemical reactions as a response to specific applied stimuli [4,[START_REF] Ma | Core-Shell hydrogel microcapsules for improved islets encapsulation[END_REF][START_REF] Fischer | Triggered release of loads from microcapsule-in-microcapsule hydrogel microcarriers: En-route to an "artificial pancreas[END_REF].

However, the design and the optimization of capsules are very challenging due to the multiple physicochemical and geometrical parameters that are involved. The release kinetics of the encapsulated material is markedly influenced by the capsule properties, which have to be optimized in order to deliver or absorb the solute at the desired rate. For example, a poorly designed capsule may release the encapsulated material too rapidly. If the aim is to deliver a drug to a patient, the drug concentration may exceed the toxicity threshold that might provoke overdose. On the contrary, if the release rate from this capsule is too slow, the efficacy of the treatment may be dramatically altered. The optimization of capsules design is thus crucial, and often requires extensive and expensive experimental studies. Mathematical modeling and numerical simulations are valuable alternative tools to better understand the physical mechanisms that govern and control mass transfer. They enable us to reduce significantly the number of expensive trial-and-error experiments and greatly help designing capsules and improving their performance.

In the first chapter of this PhD thesis, mass transfer from and to capsules is numerically studied under static conditions, i.e. in absence of flow, in order to highlight the effect of the shell (or membrane) permeability on mass transfer, independently of the contribution of an external applied flow. The emphasis is put on the release kinetics of the encapsulated solute, and especially on how the properties of the capsule's membrane affect the mass transfer efficiency. A practical equation modeling solute release from capsules as a function of their physicochemical properties is derived. Then, a novel and efficient method to characterize capsules membrane permeability from experimental data is proposed and compared with the currently used approach. The proposed empirical equation along with the permeability characterization method are user-friendly tools, and may help optimizing the design of future capsules for pharmaceutical and bioengineering applications, such as controlled drug delivery and cell culture. This chapter is organized as follows. Section 1.2 introduces the model and the problem governing equations. The numerical methods used to solve the governing equations are presented and validated in Sec. 1.3. Section 1.4 discusses the effect of the membrane permeability on the release kinetics. The proposed empirical equation describing the release kinetics as a function of the capsule physicochemical properties and the characterization method of the capsules membrane permeability are presented in Secs. 1.5 and 1.6, respectively. Finally, the conclusions are drawn in Sec. 1.7.

Mathematical formulation of the model

The capsule is modeled as a rigid composite sphere composed of a liquid core of radius R c coated by an outer polymeric thick shell or thin membrane of thickness δ, as shown in Fig. 1.2. It is immersed in a surrounding bulk medium of radius R b , which is much larger than the external radius of the capsule R = R c + δ. Both solute diffusion from the capsule to the bulk (release), and solute absorption from the bulk to the capsule (uptake) are considered in this chapter, with the corresponding mass transfer directions indicated in Fig. 1 In absence of any flow and chemical reactions, diffusion of solute species in a composite medium is unsteady and driven by Fick's second law with a space-dependent diffusion coefficient,

∂c ∂t = ∇ • (D∇c) , (1.1)
where t is the time, and c and D are, respectively, the local solute concentration and diffusion coefficient. Here, the local diffusion coefficient takes different values in the core, the membrane and the bulk, and is defined as,

D(r) =          D c , 0 ≤ r ≤ R c D m , R c < r ≤ R D b , R < r ≤ R b , (1.2)
where r is the distance from the capsule's center. As the shell is a cross-linked gel or a porous rigid wall, its diffusivity to the solute D m is lower compared to the diffusivities in the core and the bulk, if the latter are in liquid state. The problem presents a spherical symmetry, and, thus, the diffusion is assumed to be purely radial with c = c (r, t). Equation (1.1) is consequently expressed only as a function of its radial derivative terms,

∂c ∂t = 1 r 2 ∂ ∂r r 2 D ∂c ∂r , (1.3)
which can be rewritten as,

∂c ∂t = 2D r ∂c ∂r + ∂D ∂r ∂c ∂r + D ∂ 2 c ∂r 2
.

(1.4)

Initial condition

Diffusion both from and into capsules is considered in this chapter. For the release case, i.e. diffusion from the capsule to the bulk, the core is initially loaded with a uniform solute concentration C 0 , and the concentration is zero elsewhere,

c(r, t = 0) =    C 0 , 0 ≤ r ≤ R c 0, r > R c . ( 1.5)
This establishes an initial concentration gradient triggering diffusion from the core to the outer surrounding bulk. When absorption into the capsule is considered, the latter is initially free of solute, while the concentration is uniform and equal to C 0 in the bulk,

c(r, t = 0) =    0, r ≤ R C 0 , r > R . (1.6)
This condition rather triggers diffusion from the bulk to the capsule.

Boundary conditions

The resulting spatial distribution of the concentration depends also on the boundary conditions. Here, unsteady continuous concentration and mass flux boundary conditions naturally emerge in absence of any interfacial mass resistance at the core/membrane interface

r = R c , c(R - c , t) = c(R + c , t) and D c ∂c ∂r r=R - c = D m ∂c ∂r r=R + c , ( 1.7)
and at the membrane/bulk interface r = R,

c(R -, t) = c(R + , t) and D m ∂c ∂r r=R - = D b ∂c ∂r r=R + , ( 1.8)
where the superscripts (-) and (+) respectively refer to the inner and outer sides of the interfaces. There exist other boundary conditions to model mass transfer from and to capsules, depending on the operational conditions. The following Dirichlet boundary condition, c(R, t) = 0, (

is widely used to model situations where the capsule is immersed into a well-stirred solution of "infinite" volume (R b /R 1), and is referred as "perfect sink" condition. If the bulk solution is well-stirred, but has a finite volume, the condition at the membrane/bulk interface rather takes the form,

∂c ∂t = - 3R 2 R 3 b -R 3 D m ∂c ∂r r=R , ( 1.10)
modeling situations where the concentration is kept uniform in the bulk due to stirring, but still evolves in time. The derivation of this equation is given in Appendix A. These boundary conditions greatly simplify the problem by reducing it to only two domains: the inner core, and the membrane of the capsule. However, if they could in principle be fulfilled in an experiment, for example by using a magnetic stirrer to homogenize the concentration in the bulk, they are hardly realized under operational conditions.

No mass flux is imposed at r = R b when the three layers, i.e. core, shell and bulk, are considered,

∂c ∂r r=R b = 0. (1.11)
The boundary condition at the center of the capsule is,

∂c ∂r r=0 = 0, (1.12)
to fulfill radial symmetry. The governing equations presented in Sec. 1.2 are numerically solved using the finite difference (FD) method, while considering continuity in both the concentration and the mass flux at the core/membrane and membrane/bulk interfaces. To do so, the finite difference scheme proposed by Hickson et al. [START_REF] Hickson | Finite difference schemes for multilayer diffusion[END_REF] for diffusion in multilayer slabs in Cartesian coordinates is employed. It is adapted here for spherical coordinates, which introduces additional terms. The diffusion equation, Eq. (1.4), is discretized at the inner points of each subdomain (represented as closed circles in Fig. 1.3) by a central finite difference approximation in space. Considering the capsule is composed of homogeneous layers, i.e. layers of constant diffusivity, leads to,

Numerical methods

Finite difference method

dc k dt = ξ i k -1 k c k-1 -2ξ i c k + ξ i k + 1 k c k+1 , (1.13)
for the discretization at the inner points of layer i. Here, c k is the local concentration at a given inner point k, ξ i = D i ∆r 2 , and ∆r is the space step, which is uniform in the whole domain. The discretization of the domain proposed in Ref. [START_REF] Hickson | Finite difference schemes for multilayer diffusion[END_REF] assumes one grid-point j lies on the interface χ i between the layers i and i + 1, as illustrated in Fig. 1.3. Taking a central finite difference scheme to discretize Eq. ( 1.3) at the interface χ i , and then taking respectively first-order backward and forward differences for the spatial derivatives of the concentration in layers i and i + 1 gives,

dc j dt = ξ i 2 j -1 j 2 c j-1 -   ξ i 2 j -1 j 2 + ξ i+1 2 j + 1 j 2   c j + ξ i+1 2 j + 1 j 2 c j+1 .
(1.

                    . . . c j-3 c j-2 c j-1 c j c j+1 c j+2 c j+3 . . .                     =                       . . . . . . • • • ξ i j-3 j-2 -2ξ i ξ i j-1 j-2 0 0 0 0 • • • • • • 0 ξ i j-2 j-1 -2ξ i ξ i j j-1 0 0 0 • • • • • • 0 0 a -(a + b) b 0 0 • • • • • • 0 0 0 ξ i+1 j j+1 -2ξ i+1 ξ i+1 j+2 j+1 0 • • • • • • 0 0 0 0 ξ i+1 j+1 j+2 -2ξ i+1 ξ i+1 j+3 j+2 • • • . . . . . .                                           . . . c j-3 c j-2 c j-1 c j c j+1 c j+2 c j+3 . . .                     ,
(1.15) where c j is the concentration at the interface χ i located at node j (see Fig. 1

.3), ξ i = D i ∆r 2 , ξ i+1 = D i+1 ∆r 2 , a = ξ i 2 j-1 j

2

, and b = ξ i+1 2 j+1 j

2

. The concentration at iteration n + 1 is computed using a standard Euler scheme,

C n+1 = C n + ∆tMC n , (1.16)
where C n is the column vector containing the values of the concentration at each grid-point at iteration n, M is the tridiagonal matrix in Eq. (1.15) and ∆t the timestep.

The finite difference scheme is validated by comparing the obtained concentration profiles with those computed using the semi-analytical solution provided in Ref. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF], for the following parameters,

C 0 = 1 g/L, R c = 1.5 mm, δ = 0.2 mm, R b = 30 mm, D c = 3 × 10 -10 m 2 /s, D m = 0.5 × 10 -10 m 2 /s, D b = 3 × 10 -10 m 2 /s.
The resulting concentration profiles are presented in Fig. 1.4 at times t = 5, 50 and 150 min.

The figure is truncated at r/R = 3 for clarity reasons, since the concentration is almost zero beyond this point. The thin membrane is depicted by the gray area. As can be seen, excellent agreement is obtained between the numerical and semi-analytical concentration profiles, validating thus the FD scheme.

Figure 1.4: Concentration profiles taken at times t = 5, 50 and 150 min. The red solid line curves represent the numerical solution obtained with the finite difference (FD) method, while green dashed lines represent the semi-analytical solution taken from Ref. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF]. The gray solid line gives the local solute diffusion coefficient D(r). Excellent agreement is achieved between the FD and semi-analytical solutions.

The lattice Boltzmann (LB) method

Some simulations in this chapter are performed in three dimensions using the lattice Boltzmann (LB) method [START_REF] Succi | The lattice Boltzmann equation for fluid dynamics and beyond[END_REF][START_REF] Sukop | Lattice Boltzmann modeling[END_REF][START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF]. The basis of the LB method is the LB equation,

g i (r + e i , t + 1) -g i (r, t) = Ω i (r, t), (1.17)
describing the evolution in time and space of a distribution function g i (r, t), which is solved instead of solving directly the diffusion equation. Here, r = (x, y, z) is the position vector, and e i the ith discrete velocity vector direction, as illustrated in Fig. 1.5 for the D3Q19 lattice used in this chapter. Ω i is the collision operator, which is approximated by the Bhatnagar-Gross-Krook (BGK) [START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-Component systems[END_REF] operator,

Ω i (r, t) = - 1 τ D (g i (r, t) -g eq i (r, t)) , (1.18)
with τ D the relaxation time towards the equilibrium g eq i (r, t). It is related to the local solute diffusion coefficient D,

D(r) = c 2 s τ D (r) - 1 2 , ( 1.19)
where c s = 1/ √ 3 is the lattice speed of sound. As different solute diffusion coefficients are set in the core, the shell and the bulk, the relaxation time is a space-dependent quantity taking different values in each subdomain. A similar approach has been used by Kaoui et al. [START_REF] Ferrari | Analytical modeling for the heat transfer in sheared flows of nanofluids[END_REF][START_REF] Kaoui | How does confinement affect the dynamics of viscous vesicles and red blood cells?[END_REF][START_REF] Kaoui | Interplay between microdynamics and macrorheology in vesicle suspensions[END_REF][START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF] to include viscosity or diffusivity contrast in particles, and it has demonstrated its ability to capture the correct physics. The equilibrium distribution function g eq i is defined, in absence of convection, as,

g eq i = ω i c, (1.20)
where c is the local concentration computed as the zeroth moment of the distribution function,

c(r, t) = 18 i=0 g i (r, t). (1.21)
The weight factors of the D3Q19 lattice are: ω 0 = 1/3, ω i = 1/18 for i = 1 -6 and ω i = 1/36 for i = 7 -18. The LB solver is validated by comparing the obtained concentration profiles with the ones computed by the FD method. Figure 1.6 shows comparison between the numerical solutions obtained by both methods, with the following physical parameters,

C 0 = 1 g/L, R c = 1 mm, δ = 0.5 mm, R b = 12 mm, D c = 2 × 10 -10 m 2 /s, D m = 0.6 × 10 -10 m 2 /s, D b = 2 × 10 -10 m 2 /s.
For the LB simulation, the computation domain is a box of size 240 × 240 × 240, where the capsule is placed at its center. The conversion of the above physical parameters from SI units to lattice units gives for the parameters of the LB simulation: R c = 10, δ = 5, C 0 = 1, τ c = τ b = 1 and τ m = 0.65. As can be seen in Fig. 1.6, excellent matching is achieved between the LB and the FD methods.

Effect of the permeability on the solute release kinetics

The main characteristic of a capsule is its permeability to solute, which strongly depends on the internal structure of its membrane. It is dictated by multiple factors, such as the porosity and the tortuosity of the membrane, which result from the quality of the polymer cross-linking of which the membrane is made. From a macroscopic point of view, the membrane permeability is described by p = D m /δ, with D m the solute diffusion coefficient within the membrane, and δ the membrane thickness.

The effect of the membrane permeability is illustrated in Fig. 1.7, which shows the concentration profiles computed at various times for a homogeneous sphere and two capsules with permeabilities of 16 × 10 -8 m/s and 4 × 10 -8 m/s. The permeability is set to the desired value by adjusting solely D m , while keeping δ = 0.5 mm. The core radius is R c = 1 mm for all the particles and the bulk radius is R b = 12 mm. The solute diffusion coefficient in the core and in the bulk is 2 × 10 -10 m 2 /s, corresponding to the diffusion of insulin in an aqueous solution [START_REF] Morvan | Unsteady diffusion mass transfer in a microencapsulated islet of Langerhans for a bioartificial pancreas[END_REF]. The concentration profiles are smooth for the homogeneous sphere, in contrast to the case of capsules for which the concentration derivative with respect to the radial coordinate is discontinuous at the core/membrane and the membrane/bulk interfaces. This discontinuity is due to the conservation of the mass flux across the interfaces, which emerges naturally in absence of interfacial resistance. The leakage of the encapsulated solute from the core to the bulk is markedly slowed down by the presence of a membrane. The lower the membrane permeability, the slower the release. The permeability also controls the concentration boundary layer thickness around the particle. The boundary layer is thinner at low membrane permeability, meaning the solute remains close to the particle, and does not diffuse much in the bulk solution. The membrane acts, thus, as a barrier for the solute. The release kinetics is quantified by the cumulative amount of solute released,

R(t) = 1 - M t M 0 , ( 1.22)
where M t is the mass of solute remaining inside the particle at time t computed as,

M t = 4π R 0 c(r, t)r 2 dr, (1.23)
and M 0 is the initial mass of solute inside the capsule. This quantity is commonly used in pharmacokinetics to characterize the performance of drug carrier particles.

Figure 1.8 presents the release curves computed for the three particles considered in Fig. 1.7. The release kinetics is largely affected by the presence of the membrane and its permeability. For the homogeneous sphere, the solute is released instantaneously and very rapidly into the bulk solution. Almost 90% of the encapsulated material has already been released after only one hour. Homogeneous spheres do not allow the control of the rate at which the solute is released. The release process is delayed and considerably slowed down when a coating membrane is added. The initial delay is further increased when decreasing the permeability. Furthermore, the release rate dR/dt is much lower at smaller membrane permeabilities because solutes hardly diffuse in dense porous media. Capsules membranes, thus, allow us to control the moment at which the encapsulated solute starts to be released into the bulk, and they also enable to monitor the release rate. These two features make core-shell capsules very promising for targeted drug delivery applications, for which the drug must be released only once the device has reached its target, and the release rate must be controlled to maintain the drug concentration within the therapeutic window. Figure 1.8: Release kinetics computed for a homogeneous sphere and two capsules with membrane permeabilities p = 16 × 10 -8 m/s and p = 4 × 10 -8 m/s. The solute is released instantaneously by the homogeneous sphere, while the capsules enable the delay and slow down of the release process. The initial delay and the slow down are further increased when decreasing the membrane permeability.

Empirical equation for the release kinetics

A practical equation describing the release kinetics of capsules depending on their physicochemical properties would be very valuable in tailoring capsules that deliver their cargo at the desired rate.

Classical equations used in pharmacokinetics

Many theoretical and empirical equations already exist to predict the cumulative amount of drug release for a wide variety of controlled drug delivery systems [START_REF] Siepmann | Modeling of diffusion controlled drug delivery[END_REF][START_REF] Arifin | Mathematical modeling and simulation of drug release from microspheres: Implications to drug delivery systems[END_REF][START_REF] Higuchi | Rate of release of medicaments from ointment bases containing drugs in suspension[END_REF][START_REF] Ritger | A simple equation for description of solute release I. Fickian and non-Fickian release from non-swellable devices in the form of slabs, spheres, cylinders or discs[END_REF]. A review of the most commonly used equations in pharmacokinetics is given in what follows.

Higuchi equation

Takeru Higuchi is considered as the father of mathematical modeling in drug delivery. Higuchi's equation is one of the most famous equations in pharmacokinetics. It predicts that the amount of drug released from a thin film evolves as the square root of time [START_REF] Higuchi | Rate of release of medicaments from ointment bases containing drugs in suspension[END_REF],

Q = Dt (2A -C s ) C s .
(1.24)

Here, the same notation as in the original publication is used. Q is the amount of drug released per unit area at time t, D is the drug diffusion coefficient, A is the drug concentration and C s the drug solubility. However, this equation is analytically derived for a thin film and does not hold for other geometries. A more general equation is used for nonplanar systems,

R(t) = k √ t, (1.25)
where k is a constant which mainly depends on the geometry of the device and the solute diffusion coefficient. It is generally unknown prior experiments, and it has to be determined by fitting experimental data.

The Ritger-Peppas equation

The Ritger-Peppas equation [START_REF] Ritger | A simple equation for description of solute release I. Fickian and non-Fickian release from non-swellable devices in the form of slabs, spheres, cylinders or discs[END_REF] is another well-known equation in pharmacokinetics,

R(t) = kt n , (1.26)
with k a constant and n the diffusional exponent that characterizes the release mechanism. The latter depends on the geometry of the device and lies in the range 0.43 ≤ n ≤ 1. These two parameters also have to be determined after experiments by a fitting procedure. The case n = 1/2 recovers Higuchi's equation.

Reservoir with non-constant activity sources

The model of reservoir with non-constant activity sources (RNCAS) [START_REF] Siepmann | Modeling of diffusion controlled drug delivery[END_REF] is specifically derived for particles having core-shell structures. It corresponds to the solution of Fick's law of diffusion, assuming the concentration is constant in the bulk solution. This model predicts an exponential decay of the mass remaining inside the particle over time, independently of its geometry (slab, cylinder or sphere). For a spherical reservoir, which is the geometry of interest in this chapter, the release kinetics is described by,

R(t) = 1 -exp - 3RD m Kt R 2 c R -R 3 c , (1.27)
where K is the partition coefficient.

How accurate are these equations to model solute release from capsules?

The release curves predicted by these equations are compared with those computed by FD simulations in the following. The physical parameters set in the simulations are,

R c = 1 mm, R = 1.5 mm, R b = 12 mm, D c = 2 × 10 -10 m 2 /s, D m = 0.6 × 10 -10 m 2 /s, D b = 2 × 10 -10 m 2 /s.
They correspond to the release of insulin from a core-shell capsule to an aqueous solution [START_REF] Morvan | Unsteady diffusion mass transfer in a microencapsulated islet of Langerhans for a bioartificial pancreas[END_REF]. The bulk radius R b is chosen large enough with respect to the radius of the capsule to prevent saturation of the external solution. Insulin is then continuously released into the bulk until complete exhaustion of the capsule. The initial concentration is arbitrarily set to C 0 = 1 g/L in the core, whereas it is zero elsewhere. Continuous boundary conditions are considered at the core/membrane and membrane/bulk interfaces.

The resulting release profile is shown in Fig. 1.9 in log-log scale and as a function of the square root of time. The release kinetics computed for a homogeneous sphere is also plotted for comparison purposes. The release curve for the homogeneous sphere is a straight line of slope n = 0.5 in the log-log scale up to 60% of solute release. It follows the Higuchi (dashed line in Fig. 1.9b), and, a fortiori, the Ritger-Peppas models. However, these models predict nonphysical infinite amount of drug released as t → ∞. They are short-term approximations only. It is admitted that they are valid up to 60% of release [START_REF] Ritger | A simple equation for description of solute release I. Fickian and non-Fickian release from non-swellable devices in the form of slabs, spheres, cylinders or discs[END_REF], which is confirmed by the simulation. Mass transfer from capsules is much more complex. The solute first diffuses slowly through the polymeric membrane, which acts as a porous barrier and delays the release into the bulk. This initial delay is clearly visible in Fig. 1.9b. The release curve of the capsule does not follow a linear trend in log-log scale, and cannot thus be modeled by Higuchi or Ritger-Peppas models, as in the case of the homogeneous sphere. The release curve predicted by the model of reservoir with non-constant activity source, Eq. (1.27), is plotted as a function of the square root of time in Fig. 1.10 (dashed line), together with the FD numerical solution (solid line). The partition coefficient is K = 1. Equation (1.27) clearly overestimates the cumulative amount of drug released. This is due to the perfect sink conditions which increases the concentration gradient at the surface of the capsule, and thus speeds up the release process. Similar conclusions were recently reported in Ref. [START_REF] Maroto-Centeno | Coarse-grained simulations of diffusion controlled release of drugs from neutral nanogels: Effect of excluded volume interactions[END_REF]. Moreover, the initial delay induced by the membrane of the capsule is not well captured by the RNCAS model. Figure 1.10: Release kinetics of insulin from a capsule computed by the FD method (solid line) and the model of reservoir with non-constant activity source (RNCAS), Eq. (1.27), with K = 1 (dashed line). The RNCAS model clearly overestimates the amount of drug released and does not capture well the initial lag effect induced by the membrane.

Predicting solute release from capsules

The classical equations used in pharmacokinetics are not well adapted to predict solute release from capsules, as discussed in the previous section. They all fail in capturing the correct release kinetics of capsules, which is more complex than the one of homogeneous spheres because of the initial delay induced by the polymeric membrane. The aim of this section is to derive a practical equation describing the release kinetics of capsules as a function of their physicochemical properties.

The equation is derived while considering continuity of both mass flux and concentration at the core/membrane and membrane/bulk interfaces. This implies the partition coefficient K is unity, and the external solution is not subjected to stirring. Under these assumptions, mass transfer from capsules is expected to depend mainly on 6 physical parameters: the core radius R c , the capsule external radius R, the bulk radius R b , and the solute diffusion coefficients in the core D c , the membrane D m and the bulk D b . However, in most practical cases, the core of the capsule is filled with the same fluid as the bulk solution, and therefore it has the same solute diffusion coefficient (D c = D b ). Furthermore, the volume of the capsule is generally much smaller than the one of the bulk solution (R R b ). These assumptions greatly simplify the problem by reducing the number of parameters from 6 to 4 (R c , R, D m , D b ). This gives two key dimensionless control parameters: the ratio of the capsule's external radius to its core radius R/R c , and the membrane-to-bulk diffusivity ratio D m /D b . These parameters are varied between 1.1 ≤ R/R c ≤ 2 and 0.1 ≤ D m /D b ≤ 1, which are the common ranges encountered in literature. The dimensionless time is defined as

T = D b t/R 2 .
As the classical equations of pharmacokinetics are not appropriate to model the release kinetics of capsules, it is necessary to introduce other equations as a basis for the derivation. The equation describing the first-order kinetic release, as well as the Weibull and the loglogistic cumulative distribution functions (CDF) are commonly used to model solute release from various drug carrier particles [START_REF] Langenbucher | Linearization of dissolution rate curves by the Weibull distribution[END_REF][START_REF] Kosmidis | A reappraisal of drug release laws using Monte Carlo simulations: The prevalence of the Weibull function[END_REF][START_REF] Koester | Mathematical evaluation of in vitro release profiles of hydroxypropylmethylcellulose matrix tablets containing carbamazepine associated to β-cyclodextrin[END_REF][START_REF] Papadopoulou | On the use of the Weibull function for the discernment of drug release mechanisms[END_REF][START_REF] Ghosal | Mathematical modeling of drug release profiles for modified hydrophobic HPMC based gels[END_REF]. They are respectively defined by Eqs. (1.28), (1.29) and (1.30).

• First-order release:

R(T) = 1 -exp - T a F , (1.28)
• Weibull:

R(T) = 1 -exp   - T a W b W   , (1.29)
• Log-logistic:

R(T) = 1 1 + T a LL -b LL . (1.30)
The log-logistic CDF, Eq. (1.30), can be rewritten as,

R(T) = 1 - 1 1 + T a LL b LL , (1.31)
so as to have the same convenient form R(T) = 1 -f (T) as the first-order kinetics and Weibull models. The parameters a F , a W and a LL are related to the timescale of the release process. More specifically, a F and a W correspond to the time T at which 1 -1 e ≈ 63% of the solute has been released into the bulk, whereas a LL corresponds to the time required to release 50% of the solute. The parameters b W and b LL of the Weibull and log-logistic CDF are shape factors. They are strictly positive and characterize the shape of the release curve as exponential (b = 1), sigmoids (b > 1), or exponential with a steep initial slope (b < 1). These three equations are compared here in order to determine which one is the most appropriate to describe solute release from capsules. Figure 1.11 presents the release kinetics computed with the finite difference method for R/R c = 2 and D m /D b = 0.2, along with the fits obtained by the aforementioned equations. The goodness of the fits is measured by both the adjusted determination coefficient, R 2 adj , and the RMSE. The latter are listed in Tab. 1.1 for the three equations, with their corresponding fitting parameters. The log-logistic CDF produces the best fit of the numerical data, and it also captures correctly the delay induced by the membrane, in contrast to the first-order and the Weibull equations. It has the highest value of R 2 adj , which is very close to 1, and the lowest RMSE. For this reason, the log-logistic CDF is chosen as the basis equation to model solute release from capsules.

The question that arises is how to correlate the timescale parameter a LL and the shape factor b LL of the log-logistic CDF, Eq. (1.31), with the physicochemical properties of the capsule. To determine an accurate correlation of these parameters, 400 numerical simulations in total were performed over the ranges 0.1

≤ D m /D b ≤ 1 and 1.1 ≤ R/R c ≤ 2.
The values of a LL and b LL were then computed for each simulation by a fitting procedure. Their averaged values are represented as black dots in Figs. 1.12a

a LL (x, y) = α 1 x α 2 + α 3 y α 4 + α 5 x α 6 y α 7 + α 8 , (1.32a) b LL (x, y) = β 1 x β 2 + β 3 y β 4 + β 5 x β 6 y β 7 + β 8 , (1.32b)
where

x = D m /D b and y = R/R c . The values of α i and β i , i = 1, ..., 8 are listed in Tab. 1.2.
The timescale parameter a LL and the shape factor b LL have a very similar dependence on the ratios R/R c and D m /D b . They both increase when increasing R/R c and decreasing D m /D b , namely when the membrane becomes less permeable to the solute, further delaying and slowing down the release process. They are also more sensitive to variations in the diffusivities ratio than in the radii ratio. The timescale parameter is almost constant and particularly low when D m /D b > 0.5. For such high diffusivities ratios, the solute diffuses easily through the membrane, and, thus, the release process is fast, even for the thickest membrane. Oppositely, a LL sharply increases when D m /D b → 0, i.e. when the size of the membrane pores reduces. In this case, the solute hardly diffuses within the membrane and the release process is very slow. As expected, the shape factor b LL is always greater than 1, independently of the capsules design. This means the release curves are Sshaped (sigmoids), with an initial plateau that is characteristic of capsules release kinetics. Increasing the radii ratio R/R c while lowering the diffusivities ratio D m /D b increases the shape factor, and, thus, amplifies further the S-shape character of the release curve. The proposed final equation describing the release kinetics from capsules then reads,

R(T) = 1 - 1 1 + T a b , (1.33)
where the timescale parameter a and the shape factor b are functions of the diffusivities and radii ratios and are given by,

a = 0.193 D m D b -0.937 + 0.097 R R c -0.517 -0.192 D m D b -0.938 R R c -1.562 , (1.34a) b = 3.260 -1.614 D m D b 0.392 -1.346 R R c -1.789 + 0.821 D m D b 0.913 R R c -1.721 . (1.34b)
The release curves predicted by the proposed equation are compared with those computed by the semi-analytical solution of Ref. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF], and a lattice Boltzmann (LB) simulation, for two different capsule designs in Fig. 1.13. The present model is in excellent agreement with both the semi-analytical solution and the LB simulation. It is more practical and user-friendly than the semi-analytical solution, and it produces instantaneous results, in contrast to numerical simulations. It enables us to compute the release kinetics of capsules on the fly with excellent accuracy, which is very convenient for optimizing the design of capsules in order to obtain the desired release kinetics and improve their performance.

Characterization of capsules membrane permeability

The release kinetics of a capsule is well described by the equation proposed in the preceding section. It is found to strongly depend on the capsule's membrane permeability. However, the membrane permeability is an intrinsic property of the capsule and is hardly measurable. The solute diffusion coefficient D m required in the proposed equation, is, thus, a priori unknown. The question that arises is how to characterize the membrane permeability, which depends on the material and the conditions under which a capsule is fabricated. The permeability may differ from a population of capsules to another due to slight accidental changes in the fabrication process. It cannot be measured directly because of the impossibility to access the protected inner core of the capsule in order, for example, to place a probe apparatus. The closed topology of the capsules leaves only the option to access physically their external surrounding environment.

Classical method

The state-of-the-art approach to determine capsules permeability consists in measuring the evolution in time of the cumulative released mass of solute in the external surrounding medium, which is then fitted with a mathematical equation modeling the release or the uptake kinetics. This approach has been employed in many studies [2,[START_REF] Dembczynski | Characterization of small molecules diffusion in hydrogel-membrane liquid-core capsules[END_REF][START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF][START_REF] Yao | Mass transfer behavior of solutes in NaCS-PDMDAAC capsules[END_REF][START_REF] Ameloot | Interfacial synthesis of hollow metal-organic framework capsules demonstrating selective permeability[END_REF][START_REF] Henning | Sealing liquid-filled pectinate capsules with a shellac coating[END_REF][START_REF] Leick | Alginate/poly-L-lysine capsules: mechanical properties and drug release characteristics[END_REF][START_REF] Shahravan | Encapsulation and controlled release from core-shell nanoparticles fabricated by plasma polymerization[END_REF] to characterize permeability of capsules made of various materials. The fitting equation depends on the experimental conditions, such as the direction of mass transfer (release or uptake) and the volume of the bulk solution.

For the release case, the most widely considered equation is the analytical solution of Fick's second law for diffusion from a sphere to a well-stirred solution of limited volume [START_REF] Crank | The Mathematics of Diffusion[END_REF],

C t C eq = 1 - ∞ n=1 6α(1 + α) 9 + 9α + q 2 n α 2 exp - Dq 2 n R 2 t , (1.35)
where C t is the solute concentration in the bulk at time t, which is assumed to be uniform due to stirring, and C eq is the expected concentration at equilibrium, i.e. as t → ∞. R is the radius of the sphere, D is the solute diffusion coefficient in the sphere, and α = V b /V is the volumes ratio of the bulk and the sphere, assuming the partition coefficient is equal to unity. The q n 's are the nonzero positive roots of,

tan q n = 3q n 3 + αq 2 n , (1.36)
which is a transcendental equation that needs to be solved numerically. The expected equilibrium concentration C eq depends on the initial concentration C 0 within the sphere, and can be deduced from mass conservation,

C 0 V = C eq (V b + V ), (1.37)
which leads to,

C eq = C 0 /(α + 1). (1.38)
If the volume of the bulk solution is very large compared with the volume of the capsule (α → ∞), Eq. (1.35) simply reduces to,

C t C eq = 1 - 6 π 2 ∞ n=1 1 n 2 exp - Dn 2 π 2 R 2 t , (1.39)
which is more convenient since it does not imply computing the roots of a transcendental equation.

For the absorption case, namely for diffusion from the well-stirred bulk solution to the capsule, the analytical solution of Fick's second law is rather [START_REF] Crank | The Mathematics of Diffusion[END_REF],

C t C 0 = α 1 + α 1 + ∞ n=1 6 (1 + α) 9 + 9α + q 2 n α 2 exp - Dq 2 n t R 2 , ( 1.40)
where C 0 denotes now the initial concentration in the bulk solution. The derivation of Eqs. (1.35) and (1.40) is provided in Appendix A.

However, these equations have not been derived originally for a capsule with a composite structure, but simply for a homogeneous sphere. The diffusion coefficient D is then considered by most authors as the effective diffusion coefficient of the overall capsule (core and membrane combined). In this way, the fit only enables the determination of an average diffusion coefficient for the whole capsule, and not specifically for the capsule membrane. As the membrane permeability is a key parameter in designing, for example, capsules as drug carrier particles or bioreactors to encapsulate cells [START_REF] Ma | Core-Shell hydrogel microcapsules for improved islets encapsulation[END_REF], it is essential to be able to measure its value.

Only few authors have attempted to fit their experimental data with mathematical models for membranes. Kondo [START_REF] Kondo | Preparation and permeability characteristics of microcapsule membranes[END_REF] has used a model for a planar membrane by assuming negligible effects of the capsule curvature when the aspect ratio of the membrane thickness δ to the capsule radius tends to zero (δ/R → 0). Henning et al. [START_REF] Henning | Characterizing permeability and stability of microcapsules for controlled drug delivery by dynamic NMR microscopy[END_REF] have measured concentration profiles inside capsules using nuclear magnetic resonance (NMR) microscopy and have fitted their experimental data with numerical solutions computed with the Matlab pdepe solver, instead of analytical solutions of Fick's second law of diffusion. However, the pdepe solver does not allow large flexibility, especially for the boundary conditions, and does not enable to model the well-stirred condition, Eq. (1.10).

Proposed method

This section proposes a novel and accurate method to characterize capsules membrane permeability by an inverse approach. The present method is based on fitting experimental data with numerical solutions of the diffusion equation Eq. (1.1) computed by the finite difference method. It accounts for several boundary conditions (see Sec. 1.2) in order to meet various experimental and practical conditions, in contrast to the Matlab pdepe solver.

The inputs required for the characterization are the experimental data that describe either the release or the uptake kinetics of a solute, and the geometrical characteristics of the studied capsules. A wide range of the fitting parameters space (D m , D c ) is scanned, and the resulting output release (or uptake) curves are then compared against the experimental data to quantify the root mean square error (RMSE),

RMSE = 1 N N i=1 (C exp (i) -C num (i)) 2 , (1.41)
where N is the number of the input experimental data points, C exp , and C num is the computed numerical solution. Only cases where D m ≤ D c are considered because the diffusion in the membrane is expected to be slower than in the core, since the role of the membrane is to moderate mass transfer. The optimal values of the diffusivities (D m , D c ) are those leading to the minimal error, and the estimated capsule membrane permeability p is then evaluated as D m /δ.

The proposed method is first validated considering numerical data obtained by the lattice Boltzmann method. The efficiency of the proposed method is then demonstrated by applying it to available experimental data in literature, and, in particular in the cases of glucose release and uptake by capsules.

Examples of applications Application based on numerical data

First, the method is tested on numerical data obtained using a fully three-dimensional LB simulation for which the diffusion coefficients are set and known in advance release kinetics is reported in Fig. 1.15a as square symbols. The characterization method consists in taking these data, as if they were experimental data C exp , and fit them with FD numerical solutions using the solute diffusivity in both the membrane and the core (D m , D c ) as fitting parameters, while holding constant the values of R c , δ and R b . This means the release process is mainly controlled by the membrane permeability.

The release curve obtained for D c = 2.1 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s by the FD method fits almost perfectly the input LB release data, as shown in Fig. 1.15a, in contrast to the one computed for a homogeneous sphere, which leads to a poorer fit. This example has allowed us to validate the proposed method and demonstrate its efficiency to characterize the membrane permeability of capsules.

Application to experimental data of glucose release

The proposed method is now applied to determine the membrane permeability to glucose of hydrogel alginate capsules, whose release kinetics has been measured experimentally by Rolland et al. [START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF]. In this study, the size of the capsules is R = 1.73 mm, and the average membrane thickness is δ = 50 µm. The capsules are suspended in an aqueous solution whose glucose concentration is kept uniform thanks to stirring.

For this case, the computational domain in the FD simulation consists of only two layers: the core and the membrane, excluding the bulk, where the concentration C t is forced to be uniform but yet varies over time. The stirring condition is modeled by Eq. (1.10). The initial concentration of glucose is C 0 = 0.1 g/mL in the whole capsule, i.e. both in the core and the membrane. The experimental data of Rolland et al. are depicted in Fig. 1.16a as square symbols. They show the evolution in time of the glucose concentration in the liquid bulk surrounding the capsule. The fit of these experimental data while scanning the values of the diffusion coefficients in both the core and the membrane leads to the RMSE reported in Fig. 1.16b. The minimum located on the line D m = D c (symbol ⊕) corresponds to assume the capsule is a homogeneous sphere, as done in Ref. [START_REF] Rolland | Physicochemical properties of aqueous core hydrogel capsules[END_REF], whereas the minimum located within the area D m < D c (symbol ⊗) considers the composite nature m/s. The corresponding release curve fits the experimental data better than the one obtained when assuming the capsule is a homogeneous sphere, as shown in Fig. 1.16a. Moreover, Rolland et al. have pointed out that their measured effective diffusion coefficient is about two times larger than the expected one in an infinitely dilute solution, that is 6.75 × 10 -10 m 2 /s [START_REF] Gladden | Diffusion in supersaturated solutions: Glucose solutions[END_REF]. They attributed their observed increase of D to the stirring and motion of the capsule that speeds up the release. Here, the obtained value of D m is close to 6.75 × 10 -10 m 2 /s. The flow around the capsule induced by stirring is expected to enhance the release rate due to convection [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF][START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF], but not to modify and significantly increase the intrinsic permeability of the capsule membrane.

Application to experimental data of glucose absorption

The present characterization method can also be used in the case of solute absorption. Here, it is applied to the case of glucose absorption by capsules made of alginate and polyethylene glycol, using the experimental data of Koyama and Seki [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF]. In this study, the inner core radius of the capsules is R c = 2.17 mm and the membrane thickness is δ = 160 µm. The capsules were immersed into a well-stirred solution of glucose which diffuses to the inner core of the capsules through their membrane. The authors have measured the variation of glucose concentration in the surrounding bulk over time, which is reported in Fig. 1.17a as square symbols. The RMSE between these experimental data and the numerical solution obtained by the FD method while varying the fitting parameters D c and D m is shown in Fig. 1.17b. As in the case of solute release, the RMSE is very sensitive to variations of the membrane diffusivity D m , while it is not affected considerably by changes in the core diffusivity D c , especially at low membrane diffusivities. The absorption process is then also mainly controlled by the membrane permeability /s, which is very close to D = 7.9 × 10 -10 m 2 /s estimated by Koyama and Seki [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF]. The value of D m obtained by the present method is much smaller than the effective diffusion coefficient determined by Koyama and Seki, which is expected since the membrane is a porous medium where diffusion is naturally slower than in an aqueous solution. The corresponding absorption curves are shown in Fig. 1.17a. The curve for the capsule model considering an internal composite structure fits almost perfectly the experimental data, while the curve obtained assuming the capsule as a homogeneous sphere leads to a poorer fit.

The proposed method has been applied to determine the permeability of capsules in the cases of solute release and absorption, considering various boundary conditions at the surface of the particle. It allows the determination of the diffusion coefficient both in the core and the membrane, as opposed to the classically used approach, allowing only the estimation of an effective diffusion coefficient for the whole capsule. The obtained diffusion coefficients in the membrane are much smaller than the effective diffusion coefficients determined by the classical method, which is naturally expected since the membrane is a porous medium whose role is to control and slow down mass transfer. Furthermore, the present method can also recover accurately the effective diffusion coefficients estimated by the classical method. Its ease-of-use makes the proposed method an efficient tool for determining accurately the capsule membrane permeability in drug delivery applications. It can also be used as a computer-aided tool in designing future capsules with the desired release or absorption kinetics. Squares: experimental data taken from Ref. [START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF] showing the glucose absorption kinetics for capsules made of alginate and polyethylene glycol, with core radius R c = 2.17 mm and membrane thickness δ = 160 µm. Solid line curves: fits obtained while considering the capsule as a homogeneous sphere (red), and with its real composite structure (green).

Conclusions

This chapter was devoted to mass transfer from capsules under static conditions and has introduced the key elements for the next chapters, where various flow conditions will be considered.

A robust finite difference scheme is proposed to solve the unsteady diffusion equation with space-dependent diffusion coefficient in spherical coordinates. The FD solver is used to study the effect of capsules membrane permeability on mass transfer, and highlight the main differences between the release kinetics of capsules and homogeneous spheres. Capsules enable the delay and slow down of the release process, as opposed to homogeneous spheres whose solute release is instantaneous and fast. Capsules are therefore particularly interesting for controlled drug delivery applications thanks to their membrane that can be tailored to obtain the desired release kinetics.

After discussing the key role of the membrane with regard to mass transfer, two useful tools are proposed to help designing future capsules.

The first tool is a practical and user-friendly equation predicting solute release from capsules, depending on their physicochemical properties. The proposed equation is based on the log-logistic cumulative distribution function, whose timescale and shape parameters are related to the ratio of the capsule external radius to the core radius R/R c , and the membrane-to-bulk diffusivity ratio D m /D b . However, the proposed equation is derived considering continuous boundary conditions at the core/membrane and membrane/bulk interfaces, and in absence of stirring. The accuracy may degrade when the partition coefficient is not unity or the bulk solution is subjected to stirring. In such cases, the log-logistic model may still be appropriate to describe the release kinetics, but the timescale and shape parameters should be adjusted depending on the experimental conditions.

The second tool is a novel and accurate method to characterize capsules membrane permeability from experimental release or absorption curves. The proposed method is able to estimate the diffusivities in both the core and the membrane of capsules, in contrast to the currently used approach, which is limited to the estimation of an effective diffusion coefficient for the whole capsule, treated as a homogeneous sphere. It has been successfully applied in the cases of solute release and absorption with different boundary conditions, both on numerical and experimental data sets.

Chapter 2

Unsteady mass transfer from a core-shell cylinder in cross-flow

Introduction

The first chapter was devoted to mass transfer from spherical capsules in absence of flow and has highlighted the effect of the membrane permeability on the release kinetics. In this chapter, we study the impact of an external flow on mass transfer from a core-shell cylindrical fiber in cross-flow. This situation is encountered, for instance, in blood oxygenators as presented in the general introduction. Core-shell fibers are also commonly used in controlled drug release and tissue engineering [START_REF] Abdullah | Core-Shell Fibers: Design, Roles, and Controllable Release Strategies in Tissue Engineering and Drug Delivery[END_REF].

The presence of a surrounding fluid flow is known to enhance mass transfer from particles due to forced convection [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF][START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF][START_REF] Karabelas | Use of asymptotic relations to correlate mass transfer data in packed beds[END_REF]. Classical studies [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF][START_REF] Whitaker | Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles[END_REF][START_REF] Michaelides | Particles, bubbles and drops, their motion, heat and mass transfer[END_REF] have quantified the effect of the flow strength, the particles shape, and the transport properties of both the solute and the solvent on mass transfer. They propose empirical correlations for the Sherwood number Sh (the dimensionless mass transfer coefficient) as a function of the Reynolds number Re, the Schmidt number Sc, and some of them also on other parameters, Sh = f (Re, Sc, ...).

However, all these studies deal with uncoated particles (devoid of any shell) and they are limited to cases of constant and uniform concentration or mass flux at the surface of the particles. The boundary conditions considered in these studies are not appropriate to model solute release from particles, whose concentration decreases over time as the encapsulated material is released into the surrounding medium. Mass transfer boundary conditions for a core-shell reservoir are rather expected to be unsteady and nonuniform. Most existing theoretical models describing mass transfer from core-shell reservoirs with appropriate boundary conditions are limited to static conditions [START_REF] Siepmann | Modeling of diffusion controlled drug delivery[END_REF][START_REF] Arifin | Mathematical modeling and simulation of drug release from microspheres: Implications to drug delivery systems[END_REF], despite external flows are omnipresent under practical conditions. There is, thus, no systematic study on the interplay between the flow and the presence of a coating shell, whereas coating shells have numerous benefits for biomedical applications, as mentioned in the first chapter. The main advantage of coating shells is their ability to monitor and control the rate at which the solute is released [START_REF] Zarket | Onion-like multilayered polymer capsules synthesized by a bioinspired inside-out technique[END_REF][START_REF] Leick | Alginate/poly-L-lysine capsules: mechanical properties and drug release characteristics[END_REF][START_REF] Henning | Characterizing permeability and stability of microcapsules for controlled drug delivery by dynamic NMR microscopy[END_REF]. Moreover, shells significantly increase the mechanical robustness of drug carrier particles, protect their encapsulated agents from external mechanical and chemical deteriorations, and improve biocompatibility.

Modeling mass transport from composite particles such as core-shell reservoirs considering unsteady boundary conditions presents some mathematical and numerical challenges [START_REF] Siepmann | Modeling of diffusion controlled drug delivery[END_REF][START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF]. For this, the lattice Boltzmann method is used and the computational domain is decomposed into subdomains as done in Ref. [START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF]. The aim of this chapter is to systematically study the effect of the flow and the coating shell on mass transfer, and to answer the following open questions:

1. How does the transition from a flow regime to another alter mass transfer from a core-shell cylinder?

2. How does the shell solute permeability affect mass transfer?

What is the correlation for the mass transfer coefficient in the presence of a shell?

This chapter is organized as follows. The problem setup and the numerical method are presented and explained in Sec. 2.2. The transition from steady to unsteady flows and their corresponding solute spatial distribution are described in Secs. 2.3 and 2.4, respectively. The effects of the flow and the shell solute permeability are analyzed in Secs. 2.5 and 2.6. The proposed correlation for the Sherwood number in the case of a cylindrical core-shell fiber is derived in Sec. 2.7 and discussed within the context of existing classical correlations obtained for cylinders without any shell [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF]. Conclusions are given in Sec. 2.8. The original results reported in this chapter highlight for the first time the contribution of having a coating shell on mass transfer, and have been published as a regular article in the journal Physical Review Fluids, see Ref. [START_REF] Bielinski | Unsteady mass transfer from a core-shell cylinder in crossflow[END_REF].

Problem setup and numerical method

In absence of turbulence, at small and moderate Reynolds numbers Re < 200, the problem of flow perpendicular to an infinite circular cylinder can be solved in two dimensions [START_REF] Haeri | A correlation for the calculation of the local Nusselt number around circular cylinders in the range 10 ≤ Re ≤ 250 and 0.1 ≤ P r ≤ 40[END_REF]. The numerical setup is shown in Fig. 2.1. The core-shell cylindrical reservoir has an external radius R. It is composed of an inner core initially loaded in solute which is coated by a permeable thick shell of uniform thickness δ = 0.3R. The shell is considered as a porous medium that slows down and moderates diffusion of the encapsulated solute into the surrounding flowing fluid. The intrinsic permeability of the shell remains unchanged if its porous microstructure is not altered by hydrodynamical stresses. The particle is placed in a channel of length L = 20R and width W = 4R, at a distance x c = L/4 from the channel inlet and at y c = 0. The solvent enters the channel from the left with a uniform flow velocity U ∞ and with zero solute concentration c ∞ = 0. The reservoir is stationary and the hydrodynamic part of the problem corresponds to the classical case of a flow past a stationary circular cylinder. The mass transfer part is, however, less studied. The domain is decomposed into three adjacent subdomains that correspond to the core, the shell and the fluid, having respectively solute diffusion coefficients D c , D s and D f . The core of the reservoir is filled with the same fluid as the external solvent, and thus D c = D f all along this chapter. The ratio of the solute to solvent diffusivity is varied in the range 0.1 ≤ D s /D f ≤ 1, as generally encountered in literature [START_REF] Morvan | Unsteady diffusion mass transfer in a microencapsulated islet of Langerhans for a bioartificial pancreas[END_REF].

!"#$%&'"() !"#$%&'$(()*+(,-.$# θ δ d =2 R !'/--$()0/(( W L/4 U ∞ x y L r x c y c W 2 !"#$)1))))))2 3'$(()1))))))2 4(5,.)1)))))))))))))))))2 D c D s D f , ν, ρ Ω c Ω s Ω f 0 Figure 2.1:
Geometry of the problem. The solvent enters the channel from the left with uniform velocity profile, passes the stationary core-shell cylinder, and exits the channel with constant pressure. The domain is decomposed into three adjacent subdomains: the core, the shell and the fluid. The core is initially loaded with a solute that leaks out through the shell toward the outer passing fluid.

Mathematical formulation

Governing equations

The problem requires solving the Navier-Stokes equations (2.1) and (2.2) for the fluid, which is considered to be incompressible and Newtonian, along with the advection-diffusion equation (2.3) for the solute transport,

∂u ∂t + u • ∇u = - ∇p ρ + ν∇ 2 u, (2.1)
∇ • u = 0, (2.2) ∂c ∂t + u • ∇c = ∇ • (D∇c) , (2.3)
where u(x, y, t), p(x, y, t) and c(x, y, t) are respectively the local velocity of the solvent, the local pressure, and the local solute concentration at time t. ρ and ν are the mass density and the kinematic viscosity of the fluid, respectively. The solute diffusion coefficient is a space-dependent quantity and is defined as,

D(x, y) =          D c if (x, y) ∈ Ω c , D s if (x, y) ∈ Ω s , D f if (x, y) ∈ Ω f .
(2.4)

Equation (2.
3) is solved without the advection term in the cylindrical reservoir since it is stationary and there u(x, y, t) = 0.

Initial condition

Equations (2.1) to (2.3) are solved for the initial condition,

c(x, y, t = 0) =    c 0 if (x, y) ∈ Ω c , 0 elsewhere, (2.5)
that models a reservoir whose core is initially loaded with a uniform solute concentration c 0 . This establishes an initial concentration gradient that induces solute transport from the core to the bulk fluid and leads to concentration decay over time inside the reservoir until exhaustion of the solute (equilibrium state). There is no sustained supply of solute that would maintain constant concentration at the reservoir surface, as done in most classical studies. The initial flow field is a fully developed flow.

Boundary conditions

For the flow part, non-slip velocity boundary condition u = (0, 0) is set on the channel walls y = ±W/2, and on the outer surface of the cylinder r = R, r being the distance from the center of the reservoir. The flow is uniform and parallel to the walls at the inlet u(-L/2, y) = (U ∞ , 0), and it exits the channel with constant pressure. where the subscripts stand for quantities evaluated on the sides facing the core, the shell and the fluid. These unsteady continuous boundary conditions lead to different mass transfer scenario compared to the largely used constant boundary conditions, and to nonuniform boundary conditions examined recently in Ref. [START_REF] Vandadi | Reciprocal theorem for convective heat and mass transfer from a particle in Stokes and potential flows[END_REF]. They are more appropriate to model solute release from particles for which the concentration decays over time, and are then well adapted to the problem of drug delivery.

Key physical quantities

Mass transfer from a cylindrical core-shell reservoir is expected to depend on the Reynolds number Re and the Schmidt number Sc, as is the case for other types of particles, and in particular, on its shell solute permeability. The flow pattern that develops around the cylindrical reservoir is dictated by the imposed Reynolds number that measures the importance of inertia with respect to viscous forces,

Re = d U max ν , ( 2.8)
where d = 2R is the reservoir diameter, and U max = 3U ∞ /2 the maximum velocity at the channel centerline in absence of any obstacle. Because the cylinder is confined between two walls, U max is taken as the characteristic velocity, instead of U ∞ , which is usually used for unbounded geometries. In this study, it is varied in the range 0.01 ≤ Re ≤ 180 that covers Stokes flows encountered in microfluidic devices, and unsteady flows that take place in bioreactors and large blood vessels. The flow pattern is also influenced by the position of the cylinder with respect to the channel inlet (x c = L/4 here) and on the blockage degree that is defined as the ratio of the cylinder outer diameter to the channel width B = d/W [START_REF] Sahin | A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder[END_REF]. It is set to B = 0.5 (except where otherwise specified) all along this chapter to model situations with moderate effect of confinement. Moreover, mass transfer depends also on the Schmidt number, defined as the ratio between the rates of momentum and mass diffusion,

Sc = ν D f . (2.9)
It is typically of order unity for diffusion in gases, and is larger for diffusion in liquids.

Here, 1 ≤ Sc ≤ 25. For the core-shell reservoir case, mass transfer is expected to depend also on the shell permeability to solute that is characterized by the ratio of the solute diffusivity in the shell D s to the thickness of the shell δ, p s = D s /δ. It measures how easily the solute diffuses through the shell. Hereafter, it is scaled as,

P = d D f p s , ( 2.10)
and varied within the range 0.66 ≤ P ≤ 6.67 that spans from very little permeable shells made with dense polymers to highly permeable shells with large pores. Five observable quantities are measured and analyzed:

1. The solute concentration c s (θ, T) at the surface of the reservoir (r = R), which is computed using fine bilinear interpolation of the known concentration at the onlattice points. The angular position θ is defined in Fig. 2.1, with θ = 0 at the rear of the particle (downstream) and θ = π at the front (upstream). T = (D f /d 2 )t is the dimensionless time.

2. The surface mass flux,

ϕ(θ, T) = -D f ∂c ∂r r=R , ( 2.11)
that here, is unsteady and nonuniform. It evolves freely in contrast to some classical studies [START_REF] Krall | Local heat transfer around a cylinder at low Reynolds number[END_REF][START_REF] Karniadakis | Numerical simulation of forced convection heat transfer from a cylinder in crossflow[END_REF] where it is rather set as a control parameter.

3. The local Sherwood number quantifying the ratio of the convective mass transfer to the rate of diffusive mass transport at the surface of the reservoir,

Sh(θ, T) = d D f h(θ, T), (2.12)
where h(θ, T) = ϕ(θ, T)/ [c s (θ, T) -c ∞ ] is the convective mass transfer coefficient, and c ∞ the concentration at the channel inlet. It is inversely proportional to the concentration boundary layer thickness [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF].

4. The average Sherwood number, which is computed as the integration of the local Sherwood number Eq. (2.12) over the surface of the reservoir [START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF],

Sh(T) = 1 2πR 2π θ=0 Sh(θ, T)ds(θ), (2.13)
and where s is the curvilinear coordinate along the cylinder circumference.

5. The steady average Sherwood number, denoted by Sh (without T between parentheses), is used hereafter to derive a mass transfer correlation.

The lattice Boltzmann method

The lattice Boltzmann (LB) method is opted to compute both the flow and the mass transfer [START_REF] Succi | The lattice Boltzmann equation for fluid dynamics and beyond[END_REF][START_REF] Sukop | Lattice Boltzmann modeling[END_REF][START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF][START_REF] Wolf-Gladrow | Lattice-gas cellular automata and lattice Boltzmann models: An introduction[END_REF][START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF], as done in Refs. [START_REF] Ferrari | Analytical modeling for the heat transfer in sheared flows of nanofluids[END_REF][START_REF] Kaoui | Flow and mass transfer around a core-shell reservoir[END_REF][START_REF] Kaoui | Computer simulations of drug release from a liposome into the bloodstream[END_REF][START_REF] Kaoui | Algorithm to implement unsteady jump boundary conditions within the lattice Boltzmann method[END_REF], instead of solving Eqs. (2.1) to (2.3).

In the limit of small Mach numbers (ratio of the speed of a fluid particle in a medium to the speed of sound in that medium) and Knudsen numbers (ratio of the molecular mean free path to the macroscopic characteristic length scale) the LB method computes with good approximation the flow of an incompressible Newtonian fluid.

Two distribution functions f i and g i are needed with their respective lattice Boltzmann equations. The distribution f i is used for the flow part. It gives the probability to find a population of fluid particles at position r = (x, y) at time t, with the ith discrete velocity vector e i , where i depends on the lattice type. Here, the D2Q9 lattice is used, with i = 0-8. It is represented in Fig. 2.2 with its 9 discrete velocity vectors e i .

f i (r + e i , t + 1) -f i (r, t) = - 1 τ ν [f i (r, t) -f eq i (r, t)] , (2.14)
where f eq i is the equilibrium distribution function given by,

f eq i (r, t) = ω i ρ 1 + u • e i c 2 s + (u • e i ) 2 2c 4 s - u • u 2c 2 s , ( 2.15)
where c s = 1/ √ 3 is the lattice speed of sound and ω i are weight factors with ω i = 4 9 for i = 0, ω i = 1 9 for i = 1 -4, and ω i = 1 36 for i = 5 -8. The relaxation time τ ν is related to the fluid kinematic viscosity ν via the relationship,

ν = 1 3 τ ν - 1 2 .
(2.16)

The local density ρ and the local macroscopic velocity u are computed respectively by the zeroth and the first order moments of the distribution function f i ,

ρ(x, y, t) = 8 i=0 f i (x, y, t) and u(x, y, t) = 1 ρ 8 i=0 f i (x, y, t)e i .
(2.17)

The local pressure is computed as p = ρc 2 s . The zero non-slip conditions are set on the surface of the cylindrical reservoir and on the channel walls using the bounce-back boundary condition. Zou and He boundary conditions [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] are used to set uniform incoming flow u = (U ∞ , 0) at the channel inlet, and constant fluid density at the outlet in order to recover constant pressure outflow boundary condition.

For the inlet boundary condition, the wall velocity u w = (u w , v w ) is imposed but the fluid density ρ w is unknown. The distribution functions f 1 , f 5 and f 8 are also unknowns. The strategy proposed by Zou and He to solve for these 4 unknowns relies on the expressions of the moments of the distribution function,

ρ w = 8 i=0 f i = f 0 + f 1 + f 2 + f 3 + f 4 + f 5 + f 6 + f 7 + f 8 , (2.18a
)

ρ w u w = 8 i=0 f i e i,x = f 1 -f 3 + f 5 -f 6 -f 7 + f 8 , (2.18b
)

ρ w v w = 8 i=0 f i e i,y = f 2 -f 4 + f 5 + f 6 -f 7 -f 8 . (2.18c)
Because there are four unknowns, a fourth equation is required to close the problem. This equation arises from equilibrium conditions normal to the boundary,

f 1 -f eq 1 = f 3 -f eq 3 . ( 2
f 1 = f 3 + 2 3 ρ w u w , ( 2.20a
)

f 5 = f 7 - 1 2 (f 2 -f 4 ) + 1 6 ρ w u w + 1 2 ρ w v w , ( 2.20b
)

f 8 = f 6 + 1 2 (f 2 -f 4 ) + 1 6 ρ w u w - 1 2 ρ w v w , ( 2.20c
)

ρ w = f 0 + f 2 + f 4 + 2 (f 3 + f 6 + f 7 ) 1 -u w . ( 2

.20d)

A similar approach is used to set the desired mass density ρ w at the channel outlet. For the right boundary, the unknown distribution functions are f 3 , f 6 and f 7 . Solving Eqs. (2.18)(a-c) and (2.19) and assuming v w = 0 give,

f 3 = f 1 - 2 3 ρ w u w , ( 2.21a
)

f 6 = f 8 - 1 2 (f 2 -f 4 ) - 1 6 ρ w u w , ( 2.21b
)

f 7 = f 5 + 1 2 (f 2 -f 4 ) - 1 6 ρ w u w , ( 2.21c
)

u w = -1 + 1 ρ w [f 0 + f 2 + f 4 + 2 (f 1 + f 5 + f 8 )] . (2.21d)
The distribution function g i is used for the mass transfer solver with its corresponding discrete lattice Boltzmann equation,

g i (r + e i , t + 1) -g i (r, t) = - 1 τ D [g i (r, t) -g eq i (r, t)] , (2.22)
which is solved on the same D2Q9 lattice as f i . The equilibrium distribution function g eq i is defined as,

g eq i (r, t) = ω i c 1 + u • e i c 2 s + (u • e i ) 2 2c 4 s - u • u 2c 2 s , (2.23)
where c s and ω i are the same as those of the flow solver, and the relaxation time τ D is related to the solute diffusion coefficient D by, Zero mass flux is set on the channel walls using the standard bounce-back boundary condition [START_REF] Sukop | Lattice Boltzmann modeling[END_REF]. Zero concentration and zero mass flux are respectively set at the inlet and the outlet of the channel using the boundary conditions developed by Inamuro et al. [START_REF] Inamuro | A lattice Boltzmann method for a binary miscible fluid mixture and its application to a heat transfer problem[END_REF]. These boundary conditions have originally been derived for the temperature, but they can also be used for the concentration thanks to the heat/mass transfer analogy. The idea proposed by Inamuro et al. consists in assuming that the pre-collision distribution function can be represented by the equilibrium distribution function with an unknown (at the moment) concentration C ,

D = 1 3 τ D - 1 2 . ( 2
g i = g eq i (C , u w ), (2.26)
where u w is the wall velocity. The desired wall concentration c w is computed as the zerothorder moment of g i and decomposed as the sum of the known (e i • n ≤ 0) and unknown (e i • n > 0) populations [START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF],

c w = i g i = e i •n≤0 g i known + e i •n>0 g i unknown = e i •n≤0 g i + e i •n>0 g eq i (C , u w ), (2.27)
where n stands for the normal to the boundary wall. From this, one can deduce the unknown concentration C ,

C = c w -e i •n≤0 g i e i •n>0 ω i 1 + u•e i c 2 s + (u•e i ) 2 2c 4 s -u•u 2c 2 s , ( 2.28)
and insert it into Eq. (2.26) to compute the unknown populations. This technique can also be used to set the desired mass flux j on a boundary by setting,

c w = c b - j • n |x b -x w | D , (2.29)
where c b is the known concentration at node x b located along the normal direction and D is the solute diffusion coefficient.

Additional simulations are performed in the case of an uncoated cylinder whose surface is maintained at constant concentration to allow comparison with the results obtained for the core-shell cylinder, and, thus, to appreciate the effects of the shell and the continuous boundary conditions. The desired wall concentration c w is imposed on the cylinder surface using the scheme proposed by Zhang et al. [START_REF] Zhang | General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method[END_REF], which is similar to Ladd's scheme for momentum [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF],

g -i (r, t + 1) = -g * i (r, t) + 2ω i c w , (2.30)
where e -i = -e i and the superscript * refers to the post-collision state.

Validation of the numerical method

The present problem of mass transfer from a core-shell cylindrical reservoir in cross-flow does not have any analytical solution that can be used for the validation of the numerical method. The flow and mass transfer solvers have then been validated with the classical benchmark tests of Poiseuille flow and diffusion of a Gaussian hill, for which analytical solutions exist. The analytical solution for the Poiseuille flow in a rectangular channel is given by,

u x (y) = u max 1 - y H 2 , ( 2.31)
where u max is the velocity at the centerline and H the height of the channel. The analytical solution for the diffusion of a Gaussian hill is [START_REF] Crank | The Mathematics of Diffusion[END_REF], A study of the Sherwood number dependency on the grid resolution has been conducted to determine the optimal mesh resolution maximizing accuracy while keeping reasonable computational times. For large grid resolutions, the Zou and He boundary conditions [START_REF] Zou | On pressure and velocity boundary conditions for the lattice Boltzmann BGK model[END_REF] compute correctly the expected flow pattern. The stairs-like approximation of the cylinder surface inherited from the bounce-back boundary condition is smooth enough, and the wake as well as the vortices at the rear of the cylinder are well captured. Furthermore, there are enough grid points within the gaps between the cylinder surface and the channel walls for the flow to be accurately resolved. Similarly, there are enough grid points within the reservoir to resolve the diffusion precisely and to compute the bilinear interpolation of the surface concentration with good accuracy.

c(x, y, t) = M 4πDt exp - x 2 + y 2 4Dt , ( 2
The relative deviation of the computed steady average Sherwood number Sh, with respect to the finest considered resolution 2401 × 481, when coarsing the mesh is reported in Tab. 2.1 for the largest used Reynolds number (Re = 180). The measured deviations are always below 5%, even when dividing the original referential resolution by half, which is satisfactory enough to trust the computed results. The grid resolution 2001 × 401 leads to a deviation of only 0.72% with respect to the finest grid resolution. It is chosen to carry out all the simulations, with a reservoir radius R = 100 and a shell thickness δ = 30. The introduction of a coating shell has added numerical challenges to conduct stable and accurate simulations with the LB method. Achieving desired large Reynolds numbers, low shell solute permeabilities and large Schmidt numbers requires using relaxation times close to 1/2, leading to numerical instabilities. Simulations are then performed with both the Bhatnagar-Gross-Krook (BGK) and the multiple relaxation time (MRT) collision operators [START_REF] Krüger | The lattice Boltzmann method -Principles and practice[END_REF][START_REF] Bhatnagar | A model for collision processes in gases. I. Small amplitude processes in charged and neutral one-Component systems[END_REF][START_REF] Mohamad | Lattice Boltzmann method: Fundamentals and engineering applications with computer codes[END_REF].

Grid resolution

Steady and unsteady flow regimes

The transition from a flow regime to another is triggered by tuning solely the speed of the flow U ∞ to set the desired Reynolds number, while holding constant all the other parameters. Figures 2.5 to 2.7 present the flow field and the streamlines around the coreshell reservoir obtained at Re = 1, 80 and 180 corresponding to three typical flow regimes. At Re = 1, the flow is steady, laminar and unseparated. Its streamlines exhibit almost a fore-aft symmetry and are parallel to the channel walls far away from the cylinder. The velocity magnitude is maximal in the confined regions between the reservoir and the walls. For Re = 80, the hydrodynamic boundary layer detaches from the cylinder surface to form closed steady wake recirculations at the rear of the reservoir. The streamlines remain parallel to the channel walls far away downstream. At Re = 180, the flow becomes unsteady and develops periodically shed vortices, known as von Kármán vortex street. The curve computed for Re = 30 in Fig. 2.8a has a short portion of negative horizontal velocity at the rear of the cylinder. This indicates the presence of a recirculation bubble, whose size increases with the Reynolds number. Fig. 2.8b shows a zoom at the rear of the reservoir and reveals that the recirculation bubble emerges for 20 ≤ Re ≤ 22. The critical Reynolds for the first transition is then estimated to be Re c 1 = 21. According to Fig. 2.8a, the transition from steady flow with recirculations to unsteady flow with vortex shedding occurs for 125 ≤ Re ≤ 130. The critical Reynolds corresponding to the second transition is then estimated to be Re c 2 = 128. The transitions depend on the blockage degree and the distance of the cylinder from the channel inlet [START_REF] Sahin | A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder[END_REF][START_REF] Shair | The effect of confining walls on the stability of the steady wake behind a circular cylinder[END_REF][START_REF] Chen | Bifurcation for flow past a cylinder between parallel planes[END_REF]. Figure 2.9 reports the first and second transitions obtained by Chen et al. [START_REF] Chen | Bifurcation for flow past a cylinder between parallel planes[END_REF] and Sahin et al. [START_REF] Sahin | A numerical investigation of wall effects up to high blockage ratios on two-dimensional flow past a confined circular cylinder[END_REF] at various blockage ratios. The critical Reynolds of the first transition monotonically increases when further narrowing the confining walls, while the second transition towards unsteady flows is much more complex and is non-monotonic.

Solute spatial distribution at different flow regimes

The concentration is computed outside as well as inside the cylindrical core-shell reservoir, in contrast to most classical studies that are limited to solve mass transfer only outside the particle, see for example Refs. [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Michaelides | Particles, bubbles and drops, their motion, heat and mass transfer[END_REF]. Figures 2.10, 2.12 and 2.14 show the solute spatial distribution obtained at Re = 1, 80 and 180, corresponding to the three typical flow patterns previously highlighted. The two black solid circles represent the inner and the outer interfaces of the shell. For these simulations, the Schmidt number is Sc = 5 and the shell permeability to solute is P = 3.70.

For low Reynolds numbers, namely for steady laminar and unseparated flows, mass transfer is dominated by diffusion, as shown in Fig. 2.10 at Re = 1, where the solute diffuses rapidly around the cylinder, and is barely advected downstream by the flow. The concentration boundary layer is unsteady, and continues to expand until it touches the channel walls.

For 21 < Re < 128, the concentration boundary layer is thin at the front of the reservoir, and a long plume develops downstream, as illustrated in Fig. 2.12 for Re = 80. In this situation, mass transfer in the fluid is dominated by advection that transports the solute downstream until the channel outlet. The released solute at the front of the reservoir is instantaneously washed away by the newly upcoming solvent flow, leading to the emergence of a steady thin concentration boundary layer there and up to the separation points. Most of the solute is trapped in the wake vortices, where it accumulates for a while before it skips out by diffusion. Here, the isoconcentration contours around the reservoir are unsteady because neither the concentration nor the mass flux are maintained at constant values on the reservoir surface. They keep evolving until complete exhaustion of the reservoir.

For Re > 128, the solute is transported by the von Kármán vortex street, as shown in Fig. 2.14 for Re = 180. The concentration boundary layer is now even thinner at the front of the cylinder because the flow strength is larger, and so does advection. The solute never reaches the channel walls when Re > 21, and, thus, the concentration boundary layer thickness is not affected by the presence of the walls at large Re, for which the mass transfer correlation is extracted later on. The reported observations of the effect of steady and unsteady flows on the spatial distribution of the released solute by a core-shell reservoir with unsteady and non-uniform boundary conditions has not been reported elsewhere. The previously reported results computed in the limit of Stokes flows [START_REF] Clift | Bubbles, drops, and particles[END_REF], and for unsteady flows [START_REF] Haeri | A correlation for the calculation of the local Nusselt number around circular cylinders in the range 10 ≤ Re ≤ 250 and 0.1 ≤ P r ≤ 40[END_REF][START_REF] Karniadakis | Numerical simulation of forced convection heat transfer from a cylinder in crossflow[END_REF][START_REF] Bouhairi | Two-dimensional simulation of unsteady heat transfer from a circular cylinder in crossflow[END_REF] around a circular cylinder are obtained for constant and uniform surface boundary conditions. This simplifies the calculations by limiting the computation only to the outer domain.

For comparison purpose, simulations with constant concentration c w = 1 at the surface of an uncoated cylinder are also performed at the same Schmidt number Sc = 5 and Reynolds numbers Re = 1, 80 and 180. Snapshots taken at similar dimensionless times as in Figs. 2.10, 2.12 and 2.14 are shown in Figs. 2.11, 2.13 and 2.15. They allow us to appreciate the differences between having the classically imposed constant concentration boundary condition at the cylinder surface, and having the continuous boundary conditions involved by the reservoir, Eqs. (2.6) and (2.7). Mass transfer is, again, dominated by diffusion at Re = 1, while it is dominated by advection at Re = 80 and 180. However, Figs. 2.11, 2.13 and 2.15 exhibit differences with the case of the core-shell reservoir with unsteady and non-uniform boundary conditions. For example, the solute concentration in the vicinity of the cylinder and downstream is higher for the isoconcentration boundary conditions. The maximal concentration in all the snapshots is the same because the concentration inside the cylinder does not decrease over time. It is maintained constant at c 0 = c w = 1. The system evolves to a steady state at Re = 80 (Fig. 2.13), in contrast to the case of the reservoir (Fig. 2.12). For Re = 180, the released solute develops sustained periodical moving patterns downstream (Fig. 2.15), while the overall concentration decays over time in Fig. 2.14. The background applied flow is the same as in Fig. 2.12. The mass transfer reaches a steady regime, in contrast to the case of the core-shell reservoir reported in Fig. 2.12 where the concentration decays over time. The background applied flow is the same as in Fig. 2.14. The resulting solute concentration is higher and undergoes sustained oscillations compared to the case of a cylindrical core-shell reservoir reported in Fig. 2.14.

Effect of the flow on mass transfer quantities

The essential of mass transfer occurs at the surface of the core-shell reservoir, which is in contact with the ambient fluid, and where the local concentration c s (θ, T), the local mass flux ϕ(θ, T) and the resulting local Sherwood number Sh(θ, T) are found to evolve differently depending on the flow regime. At Re = 1, mass transfer is dominated by diffusion and the resulting concentration at the surface of the reservoir c s (θ, T) is almost uniform. However, at Re = 80 and Re = 180, c s (θ, T) strongly varies over the surface of the reservoir, adopting two maximal values around θ = ± π 4 , which correspond to the separation points of the concentration boundary layer. A common feature observed at all Reynolds numbers is that c s (θ, T) is always uniform at the front of the cylinder, for -3π 4 ≤ θ ≤ 3π 4 . There, the thickness of the concentration boundary layer is uniform because the solvent flow pushes the released solute against the reservoir surface. In all cases, the overall concentration in the reservoir tends to vanish at long-term as the encapsulated solute is completely released out (equilibrium state). At any given time, the local concentration is lower at higher Reynolds numbers since the solute release is enhanced by the flow. This corresponds to speeding up the exhaustion of the reservoir by increasing the flow velocity U ∞ . The reported non-uniform and unsteady evolution in time of c s (θ, T) is a characteristic feature of having a reservoir with an initial non-sustained solute load, Eq. (2.5), and continuity of mass transfer boundary conditions at its surface, Eqs. (2.6) and (2.7). Non-uniform, but steady, character of c s (θ, T) can be expected for cylinders without a shell when their surface is hold at constant mass flux. Unfortunately, these quantities are not usually reported in literature. In the present work, both c s (θ, T) and ϕ(θ, T) are unsteady and non-uniform, as opposed to the aforementioned classical studies, where only one of them is set constant, as a control parameter, and the other evolves accordingly into a steady non-uniform distribution. Here, both c s (θ, T) and ϕ(θ, T) are observable quantities and continue to evolve in time unsteadily until complete exhaustion of the reservoir from its initial solute load. Figures 2.16 and 2.17 show how mass transfer tends to equilibrium without adopting a steady regime. The surface concentration c s (θ, T) and the surface mass flux ϕ(θ, T) are used in Eq. (2.12) to compute the local Sherwood number Sh(θ, T), which is reported in Fig. 2.18. It is lower, uniform and unsteady at Re = 1 because diffusion is the dominant mass transfer mechanism. For Re = 80, Sh(θ, T) is larger at the front of the reservoir where advection is important, and lower at the rear where solute transport is operated mostly by diffusion within the recirculation zone. Sh(θ, T) does not vary in time on the upstream front of the reservoir because the released solute is constantly washed away by the flow, and, thus, the concentration boundary layer thickness is maintained constant in time. Sh(θ, T) gets minima at the separation points due to the local switch of the mass transfer mechanism from advection to diffusion. For Re = 180, Sh(θ, T) is further increased because of the amplified contribution of advection, and the further reduction in the concentration boundary layer thickness at the front of the reservoir. Sh(θ, T) is steady everywhere, even beyond the separation points, since the released solute at the rear of the reservoir is evacuated by the periodically shed vortices and does not accumulate in contrast to the case of Re = 80. Interestingly, the Sherwood number of the core-shell reservoir is steady despite the unsteadiness of both the local concentration and the local mass flux, whose instantaneous ratio is, however, a constant. The resulting local Sherwood number has a similar qualitative behavior as observed by Karniadakis [START_REF] Karniadakis | Numerical simulation of forced convection heat transfer from a cylinder in crossflow[END_REF] for a cylinder without a shell (see Fig. 2.19). However, here, the measured values are larger due to the contribution of the shell that limits solute diffusion to the reservoir surface, and, consequently, decreases the concentration boundary layer thickness.

Instantaneous average Sherwood number

The instantaneous average Sherwood number Sh(T), evaluated with Eq. (2.13), is reported in Fig. 2.20 for different flow regimes. Sh(T) is transiently infinite at the early stage of all the simulations, when the solute still diffuses through the shell and has not yet come into contact with the fluid. Later on, the three curves detach from each other when the solute reaches the surface of the reservoir. For Re = 1, Sh(T) decays continuously towards zero since the concentration boundary layer expands unsteadily as the solute is released. However, at larger Reynolds numbers, the released solute is quickly moved away from the reservoir surface by the flow, causing the concentration boundary layer to shrink into a steady thin thickness. This is why the average Sherwood number adopts a steady plateau, whose value Sh increases with Re, i.e., when the mass transfer is further enhanced by forced convection. The unsteady and oscillatory character of the von Kármán vortex street has no visible qualitative signature on the overall temporal evolution of the average Sherwood number. Fluctuations may show up if Re is further increased to trigger turbulent flows, as reported for a cylinder without a shell in Ref. [START_REF] Bouhairi | Two-dimensional simulation of unsteady heat transfer from a circular cylinder in crossflow[END_REF]. The dependency of the steady value of Sh(T) on the Schmidt number Sc and the dimensionless shell permeability to solute P is investigated in the following. The solute first diffuses from the core to the shell. During this pure diffusion process, the solute is still protected by the shell and the flow has then no impact on the release. At T ≈ 0.015, the three curves detach from each other because the solute reaches the surface of the reservoir and starts being advected by the flow. Afterwards, R(T) evolves approximately as √ T for a short time (see Fig. 2.21b). This follows the Higuchi model [START_REF] Higuchi | Rate of release of medicaments from ointment bases containing drugs in suspension[END_REF], which has originally been derived to describe drug release from particles in absence of external flow. However, this model approximates well the computed data at short-term only. At long-term, the release rate gradually slows down as the reservoir exhausts its content and R(T) deviates from the square root model. Higuchi's model cannot then be used to describe the whole release curves, as also discussed in the first chapter for a spherical capsule in absence of flow. The release is faster at larger Re due to the mass transfer enhancement by forced convection.

Solute release kinetics

Effect of the shell permeability

The main property of the shell is its permeability p s = D s /δ that quantifies how easily the solute diffuses through the pores of the shell. Understanding and characterizing the effect of the shell solute permeability on mass transfer is essential in designing efficient controlled delivery systems and devices that use core-shell or hollow fibers to control the mass transfer rate. The concentration profiles computed for Sc = 1, Re = 80 and P = 2.22 are plotted at various times in Fig. 2.22. They show similar evolution as observed in the first chapter for a spherical capsule (see Fig. This leads to discontinuities of the concentration radial derivative at the shell interfaces, and to the slowdown of the release process.

The influence of the dimensionless shell solute permeability P on the steady average Sherwood number Sh is reported in Fig. 2.23 with both linear and log-log scales. The results are also shown as a function of the diffusivity ratio D s /D f , with D f = D c . Each curve corresponds to a specific flow regime: steady laminar (Re = 10), steady laminar with recirculations (Re = 80) and unsteady with vortex shedding (Re = 180). Sh decreases in both linear and log-log scales for all the flow regimes. Figure 2.23a shows deep decline of Sh at lower values of P. These results reveal for the first time the effect of the shell solute permeability on mass transfer from a circular cylinder, coated with a semi-permeable shell, under cross-flow. The present study highlights how the internal structure of a coreshell cylinder and the transport properties of its shell alter substantially the mass transfer performance. While one would expect enhancement of Sh by increasing P because the latter enhances the mass transfer rate R(T) across the shell, data in Fig. 2.23a show the opposite scenario. For a given flow pattern, increasing P while holding all other parameters constant enhances the solute release R(T), which further expands the thickness of the concentration boundary layer to which Sh is inversely proportional. This is why Sh decreases when increasing P.

Mass transfer coefficient correlation

The aim of this section is to derive a correlation for the mass transfer coefficient accounting for the effect of the shell solute permeability. Mass transfer correlations are usually obtained by scaling experimental data of the steady average Sherwood number Sh as a function of mainly the Reynolds number Re and the Schmidt number Sc. A general form has been proposed by Hilpert for a single circular cylinder in cross-flow [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF],

Sh = cRe m Sc n . (2.33)
This correlation is adapted from the heat transfer problem to the actual mass transfer problem by replacing the Nusselt number (Nu) by the Sherwood number (Sh), and the Prandtl number (Pr) by the Schmidt number (Sc). The prefactor c depends on Re and on other factors, such as the shape of the cylinders. The value of the power exponent m depends also on Re, and it is generally within the range of 0.3 ≤ m ≤ 0.5 in absence of turbulence [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Bergman | Incropera's Principles of Heat and Mass Transfer[END_REF]. Most studies adopt n = 1 3 , as predicted by the boundary layer theory. However, in some studies, deviation from 1 3 is measured [START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF]. There are many other correlations coming with different degrees of accuracy [START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF][START_REF] Morgan | The overall convective heat transfer from smooth circular cylinders[END_REF][START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF][START_REF] Whitaker | Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles[END_REF]. The most popular one is the correlation of Churchill and Bernstein [START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF], (2.34)

Sh = 0.3 + 0.62Re 1/2 Sc 1/3 1 + (0.4/Sc) 2/3 1/4   1 + Re 282, 000

CHAPTER 2. UNSTEADY MASS TRANSFER FROM A CORE-SHELL CYLIN[...] 81

Despite its complexity compared to Hilpert's correlation, the Churchill-Bernstein correlation has the advantage of being valid for a wide range of the Péclet number Pe = Re Sc ≥ 0.2, and without involving exponents that depend on either Re or Sc. This correlation, and all others in literature, consider Dirichlet or Neumann boundary conditions corresponding to either constant concentration or constant mass flux at the surface of the cylinder. Thus, they are not adapted to the case of solute release from a cylindrical reservoir that implies unsteady continuous boundary conditions, Eqs. (2.6) and (2.7). Moreover, they do not consider the presence of a coating semi-permeable shell that introduces an additional parameter: the shell solute permeability. In this study, the general form of Hilpert, Eq. (2.33), is used and the dependency of its power exponents and prefactor on Re, Sc, and more importantly on P is examined. Two main behaviors are obtained depending on the value of Re. For Re ≤ 1, Sh is independent of Re and it is relatively low since mass transfer is dominated by diffusion. It declines very slowly toward zero as the solute is released from the reservoir and the concentration boundary layer continues to expand without adopting a steady thickness, as shown in Fig. 2.10. At long term, Sh is expected to vanish (Sh → 0 as Re → 0) as indicated in Ref. [START_REF] Whitaker | Forced convection heat transfer correlations for flow in pipes, past flat plates, single cylinders, single spheres, and for flow in packed beds and tube bundles[END_REF]. For Re 1, advection is the dominant mass transfer mechanism, and, consequently, Sh increases with the flow strength. In the range 1 ≤ Re ≤ 10, mass transfer is operated by both advection and diffusion, with none of them being predominant. The obtained Sh has similar main qualitative behaviors and features as reported for other particles and obstacles [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Michaelides | Particles, bubbles and drops, their motion, heat and mass transfer[END_REF][START_REF] Bergman | Incropera's Principles of Heat and Mass Transfer[END_REF]], but with different numerical values due to the contribution of the shell and the unsteady boundary conditions.

The contribution of the Reynolds number, the Schmidt number and the shell permeability to the mass transfer correlation is analyzed within the range 10 ≤ Re ≤ 180, in which Sh strongly depends on Re. Sh is plotted as a function of Re and Sc in log-log scale for different dimensionless shell solute permeabilities P in Fig. 2.25. Sc is varied from 1 to 25. The numerical data are represented as symbol points, and their respective fits as solid lines. For comparison purposes, numerical data obtained for an uncoated cylinder sustained at constant surface concentration are also reported (down-pointing triangles). For these data, the Reynolds number is defined as Re = U ∞ d/ν to allow quantitative comparison with the classical correlations of Churchill-Bernstein [START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF] (dashed line) and Hilpert [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF] (dotted line), which are also plotted in Fig. 2.25. For comparison purposes, data from simulations with constant surface concentration on a cylinder without a shell are also shown (down-pointing triangles), along with the correlations of Churchill-Bernstein [START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF] and Hilpert [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Bergman | Incropera's Principles of Heat and Mass Transfer[END_REF]. The presence of the shell shifts Sh upwards to larger values.

For the core-shell cylinder, the fits of the numerical data in both subfigures are linear and adopt the same slope in each subfigure regardless of the value of P. The measured slopes do not depend on either Re or Sc. Sh increases with decreasing P, as explained in the previous section. The few data points that deviate from the solid lines at either low Re or Sc, or large P are discarded in the analysis. For these parameters, the advection is relatively weak and the concentration boundary layer is altered by the channel confining walls. All the data computed for the core-shell reservoir are larger than those computed for a cylinder without a shell and with constant surface concentration. Even the largest used value of P, which corresponds to the case of a reservoir with an internal homogeneous structure (D s = D c ), leads to larger Sh. The correlations of Hilpert [START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF] and Churchill-Bernstein [START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF] both fall below all other data since they are obtained for unconfined cylinders (B → 0). The increase of the blockage ratio B from 0 to 0.5 increases Sh, but not as much as the presence of a shell does. It is remarkable how decreasing P shifts Sh upwards far from the values of the classical correlations. where the prefactor c is certainly a function of the shell solute permeability P. Other couple of exponents than (0.43, 0.40) have been tested to fit the data, but all of them lead to larger errors. The Reynolds number exponent m = 0.43 is within the range

1 3 ≤ m ≤ 1 2
that is consistent with Hilpert's correlation for the studied range of the Reynolds numbers (Re < 4000). The deviation of the Schmidt number exponent n from its classical value of 1/3 could be explained by the unsteady and continuous boundary conditions on the surface of the reservoir. The deviation may also be attributed to the effect of the wall confinement. Indeed, even when setting constant concentration boundary conditions on the surface of the cylinder with B = 0.5, the measured Schmidt exponent is n = 0.35, which is slightly larger than 1/3.

It remains to determine how the prefactor c in Eq. (2.35) depends on P in order to quantify the effect of the shell on mass transfer. In the following, this is determined for the blockage ratio B = 0.5, which is set along the present study. The slight deviation observed in Fig. 2.25 for the case of constant and uniform boundary conditions with respect to the correlations of Hilpert and Churchill-Bernstein, suggests that the effect of confinement is not as such significant at B = 0.5, as also discussed in Ref. [START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF]. Figure 2.26 gives the scaled correlation prefactor c = Sh/(Re 0.43 Sc 0.40 ) versus the shell solute permeability P, and the corresponding ratio of the diffusivity in the shell to the one in the fluid D s /D f . The represented data are obtained within the Reynolds number range of 10 ≤ Re ≤ 180 and the Schmidt number range of 1 ≤ Sc ≤ 25. All data collapse into a single line of slope l = -0.40, except few that correspond to low Péclet numbers. For these points, the contribution of diffusion to mass transfer is important and the solute reaches the channel walls, which distorts the concentration boundary layer and decreases the Sherwood number. (2.36)

This equation correlates well all the numerical data with an error of 4%, which is very satisfactory. The proposed correlation demonstrates the non-negligible contribution of the presence of a coating shell on the mass transfer coefficient. However, further systematic analysis of the dependency of c, m and n on the blockage degree is essential before proposing a universal correlation. This is left for a future investigation. Here, the effect of the blockage on the mass transfer coefficient is briefly explored for a given value of the dimensionless shell solute permeability P = 3.7, and at Sc = 5, see Fig. Sh increases with increasing B, i.e., by confining further the cylindrical reservoir between the channel walls. This is consistent with the behavior of uncoated cylinders reported in Refs. [START_REF] Zukauskas | Heat transfer from tubes in crossflow[END_REF][START_REF] Morgan | The overall convective heat transfer from smooth circular cylinders[END_REF][START_REF] Khan | Fluid flow and heat transfer from a cylinder between parallel planes[END_REF]. The confinement speeds up the fluid flow within the gaps between the cylinder surface and the walls, amplifying the contribution of convection, and, consequently, increasing Sh. The data points are perfectly fitted with a mathematical function of the form α exp (βB γ ), where the fitting parameters α, β and γ depend on Re. Estimation of Sh at B = 0 is obtained by extrapolation. Figure 2.27 shows that Sh varies slowly within the range 0 ≤ B ≤ 0.5, and it varies very much when B > 0.5 that corresponds to the strong confinement limit. Moreover, the estimated values at B = 0 are larger than the values predicted by the correlation of Hilpert, Eq. (2.33), due to the presence of the coating shell.

Conclusions

Mass transfer from a single cylindrical core-shell reservoir in cross-flow is studied numerically using two-dimensional and two-component lattice Boltzmann simulations. The main outcome is that the reservoir character leads to unsteady and non-uniform distribution of both the concentration and the mass flux at the cylinder surface. This mechanism differs from the classical studies, where either the concentration or the mass flux is set constant (as a control parameter) and the other quantity evolves to a steady non-uniform distribution. However, the resulting Sherwood number manifests the same qualitative, but not quantitative, spatial distribution as if the surface concentration or the mass flux has been set constant. Having a coating shell with low permeability is found to shrink the concentration boundary layer around the reservoir, and, thus, to increase the overall Sherwood number. The non-negligible contribution of the shell permeability to mass transfer is highlighted for the first time by the proposed correlation, Sh = 1.99P -0.40 Re 0.43 Sc 0.40 , which is extracted from a large data set computed within the ranges 10 ≤ Re ≤ 180, 1 ≤ Sc ≤ 25, and 0.66 ≤ P ≤ 6.67 for a given blockage ratio B = 0.5. This correlation complements the list of existing mass transfer correlations with the presently studied case of a cylindrical reservoir coated with a semi-permeable shell.

Chapter 3

Mass transfer from a spherical capsule in shear flow

Introduction

The second chapter dealt with mass transfer from a stationary core-shell cylindrical reservoir in cross-flow, in two dimensions. Here, the system of interest is rather a spherical capsule immersed and freely rotating in a shear flow, in three dimensions. This situation is encountered, for example, for drug release from a capsule in the gastrointestinal tract, where the contractions of the intestine walls trigger a shear flow [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF]. It is also largely encountered in microfluidic devices, such as bioreactors or organ-on-a-chip systems, where particles are used to release solute (e.g. insulin for the bioartificial pancreas developed at BMBI laboratory [4]) in the external flowing solvent. For these applications, it is necessary to understand the effect of the shear flow on mass transfer in order to optimize the design of both the particles and the microfluidic device.

Transport phenomena from solid spheres suspended in a simple shear flow has first been theoretically studied by Acrivos and co-workers [START_REF] Frankel | Heat and mass transfer from small spheres and cylinders freely suspended in shear flow[END_REF][START_REF] Acrivos | Heat transfer at high Péclet number from a small sphere freely rotating in a simple shear field[END_REF][START_REF] Poe | Closed streamline flows past small rotating particles: Heat transfer at high Péclet numbers[END_REF] in the limit of small Reynolds numbers (Re → 0), and at either low Péclet numbers (Pe → 0) or high Péclet numbers (Pe → ∞). They derived asymptotic relations for the Nusselt number (Nu) or the Sherwood number (Sh), which are respectively the heat and mass transfer dimensionless coefficients. Later, Batchelor [START_REF] Batchelor | Mass transfer from a particle suspended in fluid with a steady linear ambient velocity distribution[END_REF] used a completely different approach and obtained correlations that are similar to those of Acrivos. Polyanin and Dil'man [START_REF] Polyanin | man, New methods of the mass and heat transfer theory-II. The methods of asymptotic interpolation and extrapolation[END_REF] used an advanced fitting procedure to bridge the gap between the correlations of Acrivos in the limits of Pe → 0 and Pe → ∞. Subramanian and Koch [START_REF] Subramanian | Inertial effects on the transfer of heat or mass transfer from neutrally buoyant spheres in a steady linear velocity field[END_REF] extended the works of Acrivos and Batchelor to small, but finite Reynolds numbers, and reported the effects of inertia on heat and mass transfer. Longest and Kleinstreuer [START_REF] Longest | Interacting effects of uniform flow, plane shear, and near-wall proximity on the heat and mass transfer of respiratory aerosols[END_REF] investigated the effect of having bounding walls on the Sherwood number using numerical simulations, and proposed other correlations that take into account the contribution of the blockage ratio. More recently, Wang and Brasseur [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF] studied numerically mass transfer from a freely rotating sphere in simple shear flow, with application to drug dissolution in the intestines. They proposed correlations for the Sherwood number, including the effect of local shear, and which are consistent with those theoretically derived by Batchelor and Acrivos.

These studies consider steady and uniform boundary conditions at the surface of the particles, most of the time for simplification reasons. These boundary conditions are not appropriate to model solute release from particles, whose surface concentration is expected to be nonuniform and to decrease in time, as the encapsulated solute is released into the bulk. While there is a huge amount of studies on transport phenomena considering simple steady and uniform boundary conditions, there is only few works dealing with more realistic boundary conditions. This chapter focuses on mass transfer from a freely rotating capsule under shear flow in three dimensions, while considering unsteady and nonuniform surface boundary conditions that are relevant to drug release applications. As the contribution of a coating shell to mass transfer has already been characterized in the previous chapters, here, we consider a capsule with zero-thickness membrane and rather put the emphasis on the effect of the flow and the boundary conditions. The numerical setup of the problem is shown in Fig. 3.1. The capsule (Ω in ) is assumed to be rigid and is modeled as a liquid-filled sphere of radius R, which is suspended in the middle of a channel of length L = 30R, width W = 10R, and height H = 10R. The encapsulated fluid has density ρ in , kinematic viscosity ν in and solute diffusion coefficient D in . These physical quantities are labeled by the subscript "out" in the surrounding fluid (Ω out ). The upper and lower walls translate in opposite directions along the x-axis at a constant velocity magnitude U . This generates a simple shear flow inside the channel that induces the rotation of the capsule and its inner fluid at an angular velocity ω. The capsule is initially loaded at constant concentration c 0 , while the surrounding medium is initially free from solute. This establishes an initial concentration gradient that triggers unsteady solute diffusion from the capsule to the bulk. Continuity of both concentration and mass flux is considered at the capsule surface, in contrast to most classical studies that rather assume constant concentration within the particle, for simplicity. As the bulk volume is very large compared to the one of the capsule, the concentration far from the capsule is set to c ∞ = 0 (sink condition), as done in Refs. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF][START_REF] Longest | Interacting effects of uniform flow, plane shear, and near-wall proximity on the heat and mass transfer of respiratory aerosols[END_REF].

Problem setup and mathematical formulation

Governing equations

The encapsulated and the external fluids are assumed to be incompressible Newtonian fluids, and their flow is governed by the Navier-Stokes equations,

∂u ∂t + u • ∇u = - ∇p ρ + f + ν∇ 2 u, (3.1)
∇ • u = 0, (3.2)
where u = u(r, t) is the local fluid velocity, p = p(r, t) is the local pressure, and f = f (r, t) is the body force acting on the fluid at time t and at position r = (x, y, z). ρ and ν are the mass density and the kinematic viscosity of the fluid, respectively. Solute release from the capsule to the surrounding medium is described by the advection-diffusion equation,

∂c ∂t + u • ∇c = ∇ • (D∇c) , (3.3)
where c = c(r, t) is the local concentration, and D = D(r) is the local diffusion coefficient. Equations (3.1) to (3.3) are solved in the external medium and inside the capsule, as well.

Initial condition

The above governing equations are solved for the initial condition,

c(r, t = 0) =    c 0 , if r ∈ Ω in 0, if r ∈ Ω out . ( 3.4)
This initial condition models a capsule uniformly loaded in solute at concentration c 0 , and that is immersed in a surrounding fluid free from solute. It establishes an initial concentration gradient which triggers mass transfer from the capsule surface to the external fluid. The initial condition for the velocity field is a fully developed shear flow.

Boundary conditions

The capsule is immersed into a surrounding fluid in shear flow, which is induced by the translation of two parallel plates moving in opposite directions along the x-axis, at constant velocity magnitude U (see Fig. 3.1), u(x, y, ±H/2) = (±U, 0, 0). (

No-slip boundary conditions are applied on the surface of the capsule,

u in = u out on ∂Ω in , (3.6)
where u in and u out are respectively the velocity of the capsule's membrane and the velocity of its surrounding fluid. As the capsule center of mass is located precisely in the middle of the computational domain, the no-slip boundary condition induces the rotation of the particle around the y-axis due to local shear, without any translation. Periodic velocity boundary conditions are set at the domain vertical edges (x = ±L/2 and y = ±W/2).

For the mass transfer part, the boundaries of the domain are maintained at constant concentration c ∞ = 0, as done in Refs. [START_REF] Kaoui | Mechanistic modelling of drug release from multi-layer capsules[END_REF][START_REF] Longest | Interacting effects of uniform flow, plane shear, and near-wall proximity on the heat and mass transfer of respiratory aerosols[END_REF], to model sink conditions. Furthermore, continuity of both the concentration and the mass flux naturally emerges at the capsule's surface in absence of any interfacial mass transfer resistance,

c(R -, t) = c(R + , t) and D in ∂c ∂n R - = D out ∂c ∂n R + , (3.7)
where R -and R + respectively refer to the inner and outer sides of the interface. These unsteady continuous boundary conditions lead to different mass transfer scenario compared to the largely used constant and uniform surface concentration or mass flux. They enable to model properly the release of the encapsulated solute, in contrast to classical studies that assume the particle is sustained at constant concentration, and where mass transfer is solved only in the external medium.

Key physical quantities

The present problem is expected to be governed by two key dimensionless parameters: the shear Reynolds number Re and the Schmidt number Sc. In the following, the inner and the outer fluids are assumed to have the same physical properties, and thus, ν in = ν out = ν and D in = D out = D. The flow pattern around the capsule is dictated by the shear Reynolds number,

Re = γR 2 ν , ( 3.8)
where γ = 2U/H is the flow shear rate and R the capsule radius. It is adjusted to the desired value by solely varying the wall velocity U , while keeping all other parameters constant. The explored range of Reynolds number is 0.01 ≤ Re ≤ 1, which is relevant to most microfluidic applications, as well as to drug release, for example, in the gastrointestinal tract. Mass transfer also depends on the Schmidt number, which is defined as the ratio between the rates of momentum and mass diffusion,

Sc = ν D . ( 3.9)
The Schmidt number is set to Sc = 10 all along this study. Another important dimensionless parameter is the Péclet number,

Pe = ReSc = γR 2 D (3.10)
that measures the relative importance of advection and diffusion. The time t is normalized and expressed as,

T = Dt R 2 . ( 3.11)
Three observable local quantities are measured at the surface of the capsule to quantify mass transfer. The first one is the instantaneous surface concentration c s (θ, φ, t), which is evaluated using a fine trilinear interpolation of the known concentration at the on-lattice points. The variables θ and φ respectively refer to the colatitude and the longitude of the spherical coordinates system. From the surface concentration, one can compute the mass flux at the surface of the capsule,

ϕ(θ, φ, t) = -D ∂c ∂n r=R + . ( 3.12)
The third observable local quantity is the local Sherwood number, which is the dimensionless mass transfer coefficient,

Sh(θ, φ, t) = 2R D ϕ(θ, φ, t) c s (θ, φ, t) -c ∞ , ( 3.13)
where c ∞ is the concentration far from the capsule, which is set to zero all along this chapter. The instantaneous average Sherwood number is computed by integration of the local Sherwood number over the surface of the capsule,

Sh(t) = 1 4πR 2 ∂Ω in Sh(θ, φ, t)dS. ( 3.14)
The mathematical formulation still holds for the problem of heat transfer from a sphere in shear flow, by replacing the concentration c by the temperature T , the solute diffusion coefficient D by the thermal diffusivity α, the Schmidt number Sc = ν D by the Prandtl number Pr = ν α , and the Sherwood number Sh by the Nusselt number Nu.

Numerical methods

The present problem is a multiphysics problem. It involves fluid flow, mass transfer and structural mechanics. The basis of the numerical method used to solve this problem is schematized in Fig. 3.2. The flow field and the transfer part are both computed on an Eulerian grid using three-dimensional lattice Boltzmann simulations. There is a one-way coupling between the flow and the mass transfer solver. The flow is indeed expected to transport the solute, but the molecules of the solute are supposed to be small enough so that they do not affect the solvent flow. The membrane of the capsule is modeled by a set of massless nodes distributed on a Lagrangian mesh and connected together by a spring network. The fluid flow advects and deforms the capsule's membrane, which, in turn, exerts a force back on the fluid by relaxing to its lowest energy configuration. There is then a two-way coupling between the fluid flow and the structure, which is achieved using the immersed boundary method. The numerical method is explained in further details in what follows.

FLUID FLOW (LBM)

MASS TRANSFER

The lattice Boltzmann method

Both the flow and the concentration fields are computed in three dimensions using the lattice Boltzmann (LB) method. As in the second chapter, two distribution functions are used: f i for the flow part, and g i for the mass transfer part. However, the present problem is three-dimensional and the lattice Boltzmann equations are solved on the D3Q19 lattice rather than on the D2Q9 lattice. As the LB method has already been presented, only the outlines and the main differences with the previous chapters are mentioned in this section. Here, the capsule is expected to exert a force on the fluid. This force is taken into account by introducing a forcing term f b i in the lattice Boltzmann equation associated to f i ,

f i (r + e i , t + 1) -f i (r, t) = -(f i -f eq i ) /τ ν + f b i . (3.15)
This forcing term will be clarified in the next section describing the immersed boundary method. For the mass transfer part, there is no source term and the LB equation is the same as in the preceding chapters,

g i (r + e i , t + 1) -g i (r, t) = -(g i -g eq i ) /τ D . ( 3.16)
The equilibrium distribution functions and the expressions of the relaxation times are the same as in the second chapter. Again, the local density ρ and the local macroscopic velocity u are respectively computed as the zeroth and the first order moments of the distribution

function f i , ρ(r, t) = 18 i=0 f i (r, t), (3.17) u(r, t) = 1 ρ 18 i=0 f i (r, t)e i , (3.18)
and the local concentration is computed as the zeroth order moment of the distribution function

g i , c(r, t) = 18 i=0 g i (r, t). (3.19)
The translation of the upper and lower walls is achieved using the scheme of Ladd [START_REF] Ladd | Numerical simulations of particulate suspensions via a discretized Boltzmann equation. Part 1. Theoretical foundation[END_REF],

f -i (r, t + 1) = f * i (r, t) -6ω i ρ e i • u w c s , ( 3.20)
where e -i = -e i , u w is the desired wall velocity, and the superscript * refers to the postcollision state. The concentration is set to zero at the simulation box edges by using the scheme proposed by Zhang et al. [START_REF] Zhang | General bounce-back scheme for concentration boundary condition in the lattice-Boltzmann method[END_REF], which is similar to Ladd's scheme,

g -i (r, t + 1) = -g * i (r, t) + 2ω i c w , (3.21)
where c w is the desired wall concentration. This scheme has also been used for some simulations performed while imposing constant concentration within the capsule, for comparison with the studies mentioned in the introduction.

Immersed boundary method

The surrounding fluid exerts a force on the surface of the capsule, which, in turn, exerts a force back on the surrounding fluid. The two-way coupling between the fluid flow and the particle is achieved with the immersed boundary method (IBM) [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. The surface of the capsule is approximated by a set of marker points {r j (t)} that are distributed on a Lagrangian mesh. The Lagrangian mesh is generated by performing successive subdivisions of a regular icosahedron [START_REF] Wil | Icosphere, MATLAB Central File Exchange[END_REF], as illustrated in Fig. 3.3a. At each iteration, the edges of the polyhedron are split in their midpoint, and the newly created vertices are radially projected on the circumscribed sphere of the desired radius R. This process ensures the symmetry of the mesh. The Lagrangian nodes are interconnected by a spring network (see Fig. 3.3b) that enables to model the deformation of the capsule. The velocity of the Lagrangian markers ṙj (t) is interpolated from the known Eulerian velocity field u(r, t) computed by the LB method, ṙj (t) = r ∆x 3 u(r, t)∆(r j (t), r),

where ∆ is a smooth approximation of the Dirac distribution given by,

∆(r 1 , r 2 ) = δ(x 1 , x 2 )δ(y 1 , y 2 )δ(z 1 , z 2 ), (3.23)
with,

δ (a, b) =    1 4 1 + cos π(a-b) 2 if |a -b| ≤ 2 0 else . ( 3.24)
The new position of the Lagrangian nodes are then computed at time t + 1 using the explicit Euler scheme, r j (t + 1) = r j (t) + ṙj (t).

(3.25)

After advection, the Lagrangian nodes exert a force back on the fluid, which is extrapolated to the Eulerian mesh using the same approximation of the Dirac distribution ∆,

F(r, t) = j f j (t)∆(r j (t), r), (3.26)
f j being the total force acting on the Lagrangian marker j. It is computed as the sum of the spring forces acting on node j,

f j = i =j -e s r ij -r 0 ij r ij r ij , ( 3.27)
with e s the spring constant and r ij = r ir j . The summation is over the nodes i that are connected to node j by a spring. Superscript "0" refers to the initial nondeformed state.

The force density F is then injected into the forcing term of the lattice Boltzmann equation (3.15) such as,

f b i = ω i c 2 s (F • e i ) , (3.28)
to compute the new flow field, modified by the advection of the capsule.

Validation of the numerical method

The numerical method has been carefully validated by comparing the solutions it produces to known solutions of benchmark tests. The validation is done in three steps. First, the flow solver is validated by comparing the simulated velocity field in absence of the sphere to the analytical solution of a simple shear flow,

u = (γz, 0 , 0) , (3.29)
where γ = 2U/H is the flow shear rate. Figure 3.4 shows perfect matching between the numerical solution and the analytical one, which validates the implementation of the flow solver.

The next step of validation concerns the fluid-structure coupling, which is achieved by the immersed boundary method. The capsule is inserted in the simple shear flow and its angular velocity ω is computed at different shear rates. The shear rate is varied by solely varying the wall velocity U from 0.01 to 0.05, while keeping all other parameters constant. For an unbounded domain, and in the Stokes flow limit (Re 1), the angular velocity ω of a rigid sphere in a simple shear flow is proportional to the shear rate [START_REF] Clift | Bubbles, drops, and particles[END_REF],

ω = 1 2 γ. ( 3.30)
Figure 3.5 shows the position over time of the point P 0 that is initially located on the surface of the capsule at r = (-R, 0, 0). For these simulations, the spring constant is e s = 7, which is large enough to prevent deformation of the capsule (see Fig. 3.6b) and thus to allow comparison with the analytical solution for rigid spheres. The wall velocity is U = 0.05 l.u.

The x and z coordinates exhibit regular oscillations, while the y coordinate stays at 0, indicating the capsule does not translate along the y-axis. This confirms the mesh remains Figure 3.4: Numerically computed velocity profile (symbols) compared to the expected analytical solution (solid line). Excellent agreement is achieved between the LB and the analytical solutions. symmetric, and the code is robust and numerically stable, even after many time iterations. Similar results, but with different oscillations frequencies were obtained at all other wall velocities. The instantaneous angular position of the point P 0 is defined as,

θ = arctan z x , (3.31)
and is reported in The computed angular velocities ω are reported in Tab. 3.1 at different flow shear rates γ, and the resulting ratios ω/γ are plotted in Fig. 3.6a together with the analytical solution Eq. (3.30). The computed ratios ω/γ are all very close to the expected theoretical value 1/2, which confirms that the fluid-structure interaction is well implemented. The slight deviation from 1/2 when increasing the Reynolds number is due to the fact that the theory is valid only in the Stokes flow limit (Re 1) and for an unbounded domain.

U

Re (10 -1 ) γ (10 -4 ) ω (10 The angular velocity of the capsule while varying the spring stiffness e s is shown in Fig. 3.6b. For small stiffnesses (e s < 2), the capsule's membrane deforms under the shear forces exerted by the flow, and it adopts an ellipsoidal shape at equilibrium. This results in higher drag forces because of the less hydrodynamical shape of the capsule, and that leads to a slower angular velocity ω. At higher values of e s , the capsule can be considered as rigid and remains spherical, despite the forces exerted by the shear flow.

Finally, the mass transfer solver is validated by comparing the concentration field computed by the LB method to the one obtained with the finite difference (FD) method, in absence of any flow (Re = 0). In this case, the diffusion is purely radial and the solution has spherical symmetry. Figure 3.7 shows the LB and FD normalized concentration profiles at different dimensionless times. Both solutions are in excellent agreement, which validates the mass transfer solver.

Fully developed flow pattern

The flow patterns are very similar over the explored range of Reynolds number 0.01 ≤ Re ≤ 1. Figure 3.8 represents the typical topology of the streamlines computed at Re = 1, although the same characteristic structures are observed at all other Reynolds numbers considered in this study. Figure 3.8b shows the streamlines in the x-z plane passing by the center of the capsule, which is made transparent to appreciate the inner flow. The encapsulated fluid is driven by the rotation of the sphere due to shear, and its streamlines are closed and ellipsoidal. There is a region of closed streamlines around the particle that tends to elongate with distance. As closed streamlines hinder the transport by advection, mass transfer is expected to be diffusion-dominated in this region. The region of closed streamlines is surrounded by a pair of recirculating wakes, which originates from inertial forces [START_REF] Subramanian | Inertial effects on the transfer of heat or mass transfer from neutrally buoyant spheres in a steady linear velocity field[END_REF]. These recirculations are expected to transport the solute far away downstream and to enhance mass transfer by forced convection. The streamlines remain parallel to the channel walls far away from the particle, as in a simple shear flow in absence of any obstacle. The rotation of the capsule generates a spiraling flow in its neighborhood as illustrated in Fig. 3.8c. The streamlines emerge from the surface of the sphere and carry away the solute following spirals up to the far field flow, where it is advected downstream.

Effect of the shear flow on mass transfer

Effect of the flow on the solute spatial distribution

The impact of the flow on the solute spatial distribution is depicted in Fig. 3.9, which represents the concentration isocontours inside as well as outside the capsule at Re = 0.01, 0. The concentration, the mass flux and the Sherwood number at the surface of the capsule are plotted in Fig. 3.10 at Re = 0.01 (left), 0.4 (middle) and 1 (right). The surface concentration c s is evaluated by a refined trilinear interpolation of the on-lattice known concentration, and the surface mass flux ϕ and the local Sherwood number Sh are computed from Eqs. (3.12) and (3.13), respectively. At Re = 0.01, these quantities are uniform all over the surface of the capsule, since mass transfer is mainly due to diffusion. For higher Reynolds numbers, the contribution of advection to mass transfer becomes important and the resulting local quantities show strong variations all over the surface of the capsule. This contrasts with the classical studies where the surface concentration or the mass flux are assumed to be uniform. The qualitative behavior of the surface quantities is very similar for Re = 0.4 and 1. The surface concentration is maximal in the two dark red opposite areas that are located downstream, where the solvent velocity is relatively low. There, diffusion is the dominant mass transfer mechanism and the released solute is not transported by the flow, it remains and accumulates in the region of closed streamlines (see Fig. 3.8) that surrounds the capsule. The mass flux is minimal in these areas because they are rapidly saturated in solute. The areas of maximal concentration are shifted towards the channel centerline (z = 0) when the Reynolds number is increased, as also depicted in Fig. 3.9.

Surface mass transfer quantities

On the contrary, the concentration is minimal on the two diametrically opposed surfaces that are the most exposed to the flow (dark blue). They are located at approximately 90 degrees from the regions of maximal concentration. In these regions, advection prevails and the newly released solute is efficiently transported downstream by the upcoming flow. The mass flux is then particularly high there. Consequently, the resulting Sherwood number is maximal on the most exposed surfaces, where advection is dominant, and it is minimal in the two regions where diffusion prevails and the solute accumulates.

Instantaneous Sherwood number

The instantaneous average Sherwood number is computed from the integration of the local Sherwood number over the surface of the capsule (see Eq. (3.14)). It is plotted in Fig. 3.11 at Re = 0.01, 0.4 and 1. The Sherwood number is initially particularly high because the initial condition Eq. (3.4) establishes a steep concentration gradient at the particle surface. It then dramatically drops as the solute diffuses from the particle to the surrounding fluid that rapidly increases the concentration boundary layer thickness around the capsule. Later on, the boundary layer reaches a steady thickness and consequently the resulting Sherwood number adopts a steady value. This steady value increases with the Reynolds number, because mass transfer is enhanced by forced convection and local shear in the neighborhood of the particle. The instantaneous Sherwood number shows similar evolution in time as reported for a core-shell cylinder in crossflow in the second chapter.

Steady Sherwood number

The steady values of the average Sherwood number are reported in Fig. 3.12 as a function of the Reynolds and the Péclet numbers in logarithmic scale, at Sc = 10. For comparison purposes, the values of Sh computed when considering constant concentration at the surface of the capsule are also represented, along with the correlations of Wang and Brasseur [START_REF] Wang | Enhancement of mass transfer from particles by local shearrate and correlations with application to drug dissolution[END_REF],

Polyanin and Dil'man [START_REF] Polyanin | man, New methods of the mass and heat transfer theory-II. The methods of asymptotic interpolation and extrapolation[END_REF], and Longest and Kleinstreuer [START_REF] Longest | Interacting effects of uniform flow, plane shear, and near-wall proximity on the heat and mass transfer of respiratory aerosols[END_REF], which are defined below.

• Correlation of Wang and Brasseur: Sh = Sh 0 + 0.03748Pe 0.674 Re 0.583-0.032 ln Pe , (

where Sh 0 is the Sherwood number in the limit of Re → 0, given by,

Sh 0 =          2
Sh = 2 + 0.52Pe 1/2 1 + 0.057Pe 1/2 , Re → 0, 0 < Pe < ∞.
(3.34)

• Correlation of Longest and Kleinstreuer: where Re d = γd 2 ν , h is the half distance between the moving walls, and d the particle diameter. This correlation is valid for 0 ≤ Re d ≤ 32, and 1 ≤ h/d < 10.

Sh = C 1 + C 2 Re C 3 d Sc 0.
For both boundary conditions, two distinct behaviors are observed depending on the Péclet number. When Pe < 1, mass transfer is dominated by diffusion and the Sherwood number is barely affected by the flow, it is almost constant and relatively low. For Pe > 1, advection becomes the dominant mass transfer mechanism and the resulting Sh increases linearly in the logarithmic scale. These two characteristic behaviors have also been reported for other types of particles subjected to various flow conditions in Refs. [START_REF] Clift | Bubbles, drops, and particles[END_REF][START_REF] Bielinski | Unsteady mass transfer from a core-shell cylinder in crossflow[END_REF][START_REF] Karabelas | Use of asymptotic relations to correlate mass transfer data in packed beds[END_REF]. The Sherwood number strongly depends on the boundary conditions set at the surface of the capsule. Continuous boundary conditions lead to significantly lower Sherwood numbers than the classically considered constant boundary conditions, because the concentration gradient at the surface of the capsule is lower when the latter is not sustained at its initial concentration. This is illustrated in Fig. 3.13 which shows snapshots obtained at the same Reynolds numbers as in Fig. 3.9 when considering Dirichlet boundary conditions. The qualitative behavior is the same as reported for continuous boundary conditions in Fig. 3.9, but, here, the concentration around the capsule is higher since its surface is maintained at the concentration c 0 , resulting in a higher concentration gradient.

However, continuous boundary conditions lead to a stronger dependence of Sh on the Reynolds number in the advection-dominated regime (Pe > 1), where it scales as Re 0.59 . The Sherwood number rather scales as Re 0.19 in this regime if the capsule is maintained at constant concentration, which is consistent with the correlation of Wang and Brasseur that gives a Reynolds exponent of 0.187. For these boundary conditions, it tends to Sh 0 = 2 as Re → 0, which is the expected value for pure diffusion from a sphere whose surface is sustained at constant concentration [START_REF] Clift | Bubbles, drops, and particles[END_REF] (see derivation in Appendix B).

The dependency of Sh for the case of constant boundary conditions is consistent with the correlation of Wang and Brasseur. The present numerical data adopts the same shape as the curve representing the correlation of Wang and Brasseur over the whole range of Reynolds number. However, the present numerical values are slightly lower than the one predicted by Wang and Brasseur, especially in the advection-dominated regime, where they are in closer agreement with the correlations of Polyanin and Dil'man and Longest and Kleinstreuer, with a small deviation of about 10%. This deviation is not significant compared to the large difference between the Sherwood numbers computed when considering constant or continuous boundary conditions. The qualitative behavior is similar as the one reported for continuous boundary conditions in Fig. 3.9, however, here, the overall concentration is higher since the surface of the capsule is sustained at constant concentration c 0 , leading to larger surface gradients.

Conclusions

The proposed numerical approach is based on the two-component and three-dimensional lattice Boltzmann method to compute both the flow and the concentration fields, and the immersed boundary method to achieve the fluid-structure interaction. It enables to model unsteady and nonuniform boundary conditions on the particle surface. These boundary conditions differ from the constant and uniform surface concentration (Dirichlet boundary condition) or mass flux (Neumann boundary condition), which are classically used in literature. They are more adapted to model solute release from particles, during which the concentration of the encapsulated material decreases over time. The numerical method is carefully validated using benchmark tests for which analytical or numerical solutions are known. It is used to model mass transfer from a capsule freely suspended in a simple shear flow. The velocity and concentration fields are computed outside as well as inside the particle, and the impact of the flow strength on the solute spatial distribution is analyzed.

The effects of the flow on the surface concentration, mass flux and Sherwood number are also reported and discussed. The flow is found to enhance the mass transfer efficiency by forced convection. The effect of the boundary conditions set at the surface of the capsule on mass transfer is also investigated. Having continuity of both concentration and mass flux at the surface of the capsule leads to significantly lower Sherwood numbers than when the surface is maintained at constant concentration, as classically done in literature. Continuous boundary conditions also lead to a stronger dependency of the Sherwood number on the Reynolds number. In the advection-dominated regime, Sh ∝ Re 0.59 for continuous boundary conditions, while Sh ∝ Re 0.19 for Dirichlet boundary conditions that is consistent with the correlation of Wang and Brasseur. However, the computed Sherwood number for the case of constant boundary conditions were closer to the correlations of Polyanin and Dil'man and Longest and Kleinstreuer. The present method models accurately solute release from a capsule in simple shear flow, with realistic boundary conditions at the surface of the particle, and has then many potential applications in biotechnology, such as controlled drug delivery, cell culture, or for the optimization of mass transfer in bioreactors and organ-on-a-chip devices.

Chapter 4

Squeezing multiple capsules into a constriction: transition to clogging

The previous chapters considered solute transport from a single particle under various flow conditions, and has highlighted the contribution of both the flow strength and the membrane permeability on the mass transfer efficiency. Here, we aim at studying the situation of multiple soft capsules flowing into a microfluidic constriction, which is encountered in many biomedical applications. The emphasis is put on the collective motion of the capsules, and more specifically on the formation of clogs at the constriction entrance.

Introduction

The flow of particles in microfluidic constrictions is encountered in many biomedical applications. For example, narrow constrictions are used in microfluidic devices to sort living cells or particles [START_REF] Mcfaul | Cell separation based on size and deformability using microfluidic funnel ratchets[END_REF], to measure their deformability [START_REF] Myrand-Lapierre | Multiplexed fluidic plunger mechanism for the measurement of red blood cell deformability[END_REF], or to characterize their breakup [START_REF] Goff | Squeezing bio-capsules into a constriction: deformation till break-up[END_REF]. All these applications are represented in Figs. 4.1 to 4.3, where details on the devices and their utility are given in the captions. Constricted microfluidic channels have also been recently proposed to mechanically trigger and enhance drug release from particles [START_REF] Sharei | A vector-free microfluidic platform for intracellular delivery[END_REF]. This novel technique is illustrated in Fig. 4.4. It consists in exploiting a flow to push and force capsules to pass through a narrow constriction. The capsules are squeezed into the constriction, which enlarges their membrane pores size and, consequently, increases their permeability and promotes the release of the solute into the surrounding bulk. This delivery technology can be used, for example, to release locally a fat dissolving agent into a blood vessel that would be constricted by accumulation of grease along its walls. However, the performance of such devices can be significantly hindered by clogging events that take place at the constriction entrance, especially when the particle size is of the same order as the channel width, or when multiple particles arrive suddenly at the narrow passage [START_REF] Dressaire | Clogging of microfluidic systems[END_REF][START_REF] Hong | Clogging of soft particles in two-dimensional hoppers[END_REF][START_REF] Zhang | Particle squeezing in narrow confinements[END_REF][START_REF] Van Zwieten | From cooperative to uncorrelated clogging in cross-flow microfluidic membranes[END_REF]. For instance, Luo et al. [START_REF] Luo | Solute release from an elastic capsule flowing through a microfluidic channel constriction[END_REF] reported no solute release when the capsules remained trapped at the constriction inlet. They also observed an enhancement of the release when narrowing the constriction walls, while maintaining the aperture large enough to prevent clogging. The release is, thus, expected to be affected by the mechanical properties of the particles and the channel geometry, as also reported in Ref. [START_REF] Zhang | ATP release by red blood cells under flow: Model and simulations[END_REF] in the case of adenosine triphosphate release by red blood cells. There must be then an optimal constriction width maximizing, at a given particle deformability, the release of the solute while preventing the formation of clogs. Thus, it is essential to understand the clogging phenomenon in microfluidic constrictions in order to optimize the design of these devices and improve their efficiency. Furthermore, clogging events may also occur in blood vessels if those are constricted, by grease for example. The formation of clogs in a vessel may obstruct the flow of red blood cells and other particles, such as drug carrier particles, and that can have dramatic consequences. This is another example showing the importance of studying and being able to predict the formation of clogs by soft particles in constricted channels.

Clogging events are expected to depend mainly on the particles mechanical and physical properties, as well as on the channel geometry. Recent experimental studies have investigated the clogging phenomenon by rigid particles [START_REF] Marin | Clogging in constricted suspension flows[END_REF][START_REF] Souzy | Transition from clogging to continuous flow in constricted particle suspensions[END_REF], which are found to form a stable arch at the constriction entrance, and thus, block the channel. The neck-toparticle size ratio is the leading parameter that determines the transition threshold to clogging, with permanent clogs systematically formed for a neck-to-particle size ratio below 3 (see Fig. 4.5). For sufficiently large neck-to-particle size ratios, particles flow either continuously or intermittently depending on the solid volume fraction. The corresponding dynamics have been characterized by stochastic mathematical models. Jäger et al. [START_REF] Jäger | Clogging at pore scale and pressure-induced erosion[END_REF] also studied the formation of clogs by multiple rigid particles in a micrometric pore. They determined relations between the hydraulic pressure that triggers erosive bursts and the particles diameter and the pore size. In the case of soft particles, the existing literature focuses solely on the flow and deformation of a single isolated particle passing through a constriction [START_REF] Goff | Squeezing bio-capsules into a constriction: deformation till break-up[END_REF][START_REF] Rorai | Motion of an elastic capsule in a constricted microchannel[END_REF][START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF][START_REF] Luo | Off-center motion of a trapped elastic capsule in a microfluidic channel with a narrow constriction[END_REF][START_REF] Lei | Transport mechanism of deformable microgel particle through micropores with mechanical properties characterized by AFM[END_REF]. Kusters et al. [START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF] determined phase diagrams indicating the passage and blockage of a single particle in a narrow constriction (smaller than the capsule diameter at rest) as a function of the confinement degree, the length of the neck and the intensity of the applied body force that triggers the flow (see Fig. 4.6). However, these phase diagrams have been established for a single soft particle and may not be applicable to the case of multiple soft particles. While there is an increasing need to improve the capability of microfluidic devices to handle high particle throughput, this study is carried out to address the scenarios that emerge when multiple soft capsules are pushed to pass through an abrupt constriction. This chapter extends the previous works and gets one step toward mimicking living cells in microfluidic chips by considering a suspension of soft capsules. The collective motion of the capsules is studied using fully three-dimensional simulations, as done in Ref. [START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF] for a single soft particle. Over 200 highly accurate numerical simulations were performed, while varying two control parameters: the particle deformability, and the constriction width. The effect of these two parameters on the dynamics of the particle passage is reported and analyzed. In addition to that, a state-diagram indicating whether or not particles clog the channel depending on their deformability and the ratio of their diameter to the constriction width is determined. This chapter is organized as follows. The model and the numerical method are described in Sec. 4.2. Then, the transition from the particle passage to blockage, and the dynamics of particle passage are analyzed in Secs. 4.3 and 4.4, respectively. Conclusions are given in Sec. 4.5.

Model and numerical method

The numerical setup is shown in Fig. 4.7. The dimensions of the microfluidic channel are L = 37.5d (length), W = 6.25d (width) and H = 1.875d (height), where d is the diameter of the capsules at rest, when they adopt a spherical equilibrium shape. The particle diameter is set to d = 8 lattice units (l.u.) all along this study. Both the suspending and the encapsulated fluids are considered to have the same kinematic viscosity ν and mass density ρ, as done in chapter 3. The channel contains 38 particles that are initially nondeformed and located in the inlet compartment. All the particles have the same geometrical and mechanical properties (i.e. monodisperse suspension). The constricted passage has a length L obs = 6.25d and width D, which is a control parameter. It is placed at a distance L in = 12.5d from the channel entrance. The problem is expected to depend mainly on two dimensionless control parameters: the aperture D/d, and the particle deformability that is quantified by the capillary number,

Ca = ρν d 2 γ κ s , ( 4.1)
with γ = 4u max /W the measured shear rate at the wall and κ s the shear elastic modulus.

Here, u max = fzW 2 8ρν is the velocity at the channel centerline in absence of particles and obstacle, with f z the applied body force in the z-direction generating the flow. The dimensionless elastic modulus of the membrane K = 1/Ca is also used to express the particle rigidity. The aperture D/d is set to desired value by varying solely D, while holding d constant. The values of Ca are adjusted by varying only κ s , while keeping all the other parameters constant. The Reynolds number is defined as,

Re = u max d 2 ν . (4.2)
It is set to Re = 0.1 all along this study, as encountered in microfluidic devices. The simulations of this chapter are performed using the LB3D code [START_REF] Schmieschek | LB3D: A parallel implementation of the Lattice-Boltzmann method for simulation of interacting amphiphilic fluids[END_REF].

Fluid flow solver

Here again, the flow field u(r, t) is computed using three-dimensional lattice Boltzmann (LB) simulations performed with the D3Q19 lattice and the Bhatnagar-Gross-Krook collision operator. The dynamics of the distribution function f i is described by the same lattice Boltzmann equation as in chapter 3,

f i (r + e i , t + 1) -f i (r, t) = Ω i (r, t) + F i (r, t). (4.3)
Here, the only difference with the previous chapter is the forcing term F i , which is computed as,

F i (r, t) = 1 - 1 2τ ω i e i -u c 2 s + e i • u c 4 s e i • F(r, t), (4.4)
where F(r, t) accounts for either the body force or the membrane forces. The weights ω i , the discrete velocities e i and the lattice speed of sound c s are the same as those used in the third chapter. No-slip boundary conditions are imposed on the channel walls using the bounce-back boundary conditions.

Fluid-structure interaction

As in the third chapter, the two-way coupling between the fluid flow and the capsule dynamics is achieved by the immersed boundary method [START_REF] Peskin | Numerical analysis of blood flow in the heart[END_REF][START_REF] Peskin | The immersed boundary method[END_REF]. However, here, there are multiple capsules in the channel and they may be highly deformable depending on the value of Ca set in the simulations. Their membranes are then rather discretized into 1280 triangular elements to capture accurately large deformations. This differs from the preceding chapter where the capsule was assumed to be rigid, and its membrane was discretized into 320 triangular faces to reduce computation times. The advection and reaction steps are computed using the same smooth Dirac function ∆. The main difference with chapter 3 is the membrane force, which here, is not evaluated using a spring network, but it is rather computed using a more sophisticated mechanical model that is described in the next section.

Mechanics of particle deformation

Capsules are strain-softening with zero-thickness membranes that exhibit Neo-Hookean 2D hyperelastic mechanical behavior with the energy [START_REF] Barthes-Biesel | Motion and deformation of elastic capsules and vesicles in flow[END_REF],

E s = κ s 2 A I 1 -1 + 1 I 2 + 1 dA, (4.5)
where κ s is the shear elastic modulus,

I 1 = λ 2 1 + λ 2 2 -2 and I 2 = λ 2 1 λ 2 2 - 1
are the two deformation invariants, λ 1 and λ 2 are the principal stretching ratios, A is the surface of each capsule, and dA the surface element. In addition, we enforce the constraint of capsule volume conservation by using an energy that gives the cost of any deviation of the volume from its original value as,

E v = κ v 2 (V -V 0 ) 2 V 0 , ( 4.6)
where κ v is a numerical parameter whose value is taken to be large enough to fulfill the volume conservation constraint. V 0 is the volume of the spherical capsule at rest, and V is the actual volume of the capsule. Non-physical wrinkles may emerge at the surface of capsules, and which we avoid by applying a bending force F b that is derived as a functional derivative of the Helfrich free energy proposed originally for lipid membranes [START_REF] Helfrich | Elastic Properties of Lipid Bilayers: Theory and Possible Experiments[END_REF],

F b (x i ) = 2κ b [2H(H 2 -K) + ∆ s H]n, (4.7)
where κ b is the bending modulus, H = 1 2 2 i=1 c i is the mean curvature, K = 2 i=1 c i is the Gaussian curvature, c i is the principal curvature, ∆ s is the Laplace-Beltrami operator and n the normal vector pointing outward from the membrane. The bending modulus is chosen such as the dimensionless number B = κ b /(κ s r 2 ), quantifying the relative importance of the bending rigidity with respect to the shear elasticity is small. Thus, the bending force does not affect the global capsule dynamics and deformation. H, K and ∆ s are computed following a discrete differential geometry operators approach, as described in Ref. [START_REF] Guckenberger | Theory and algorithms to compute Helfrich bending forces: A review[END_REF].

The membrane of each capsule is discretized into triangular elements, and the force acting on each membrane node x i , with i referring to the index of the node, is evaluated following the principle of virtual work such as,

F α (x i ) = - ∂E α ∂x i . ( 4.8)
The subscript α denotes either "s" for the strain energy or "v" for the volume energy. The derivatives involved in evaluating the membrane forces are computed numerically using the finite element method as described in Ref. [START_REF] Krüger | Computer Simulation Study of Collective Phenomena in Dense Suspensions of Red Blood Cells under Shear[END_REF].

Transition from particle passage to blockage

The transition from particle passage to blockage depends on geometrical properties, such as the confinement degree and the length of the neck, as reported by Kusters et al. for a single particle [START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF], and the angle of the constriction with respect to the channel walls [START_REF] Lopez-Rodriguez | Effect of hopper angle on granular clogging[END_REF]. It is also expected to depend on the applied flow strength and the particle deformability [START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF].

The volume fraction is an additional parameter which also influences the particle passage [START_REF] Marin | Clogging in constricted suspension flows[END_REF]. The clogging phenomenon is then a complex problem involving a large number of parameters. Here, we focus the analysis on the effect of the confinement and the capillary number.

The typical flow pattern in the microfluidic constriction, in absence of the capsules, is shown in Fig. 4.8 at D/d = 1.25, which is the smallest aperture considered in this study. The colorbar indicates the normalized velocity magnitude and the white solid lines represent the flow streamlines. They are parallel to the channel walls far away from the narrow passage and exhibit symmetry with respect to the channel centerline. There is no recirculations in the corners because the Reynolds number is chosen small enough so that the boundary layer does not detach. The velocity magnitude is greatly increased in the narrow passage, where one would expect an acceleration of the capsules. In contrast, it is particularly low in the corners, where capsules may eventually remain trapped. Snapshots showing both the effect of varying the aperture D/d and the particle deformability Ca are given in Fig. 4.9. For the left-hand panel, the aperture degree is increased from 1.25 to 3.25, and the capillary number is set to Ca = 0.001 modeling almost rigid capsules. For D/d = 1.25, which is the smallest aperture considered in this study, two particles are found to be sufficient to permanently clog the constriction entrance. The clogging event generally occurs at the very beginning of the simulation, and few capsules are able to pass to the outlet compartment (only one here). As the capsules are rigid, they do not (or barely) deform under the pressure exerted by the flow and the other particles that accumulate at the constriction. At moderate apertures, particles still permanently clog the constriction by building a stable arch. The stability of the arch is mainly due to the pressure applied by the flow and the forces exerted from a capsule to another. However, the clogging event occurs later than for the case of D/d = 1.25 and more particles are able to pass through the constriction, as shown here for D/d = 2. For D/d = 3.25 (the widest constriction), the aperture is large enough to prevent the formation of an arch by the capsules. Such wide apertures allow to maintain a continuous flux of particles without observing any clogging event.

On the right-hand panel of Fig. 4.9, the aperture is set to D/d = 1.25, and the capillary number is varied from Ca = 0.001 to 0.1 that covers both rigid and highly deformable capsules. Again, for small values of Ca, namely for rigid, or very little deformable particles, two capsules are sufficient to form a permanent clog at the constriction entrance. However, the number of particles that can pass to the outlet compartment increases with the deformability. When Ca is further increased, capsules become deformable enough so that they can squeeze and pass easily through the constriction. No clogging event is observed for particles with large deformation capability under flow. At large Ca and narrow constriction, the capsules cross the constriction one by one by forming a regular train. They exit the constriction as a jet and recover their initial equilibrium spherical shape far away from the narrow passage. The reported results are observed independently of the initial position of the particles.

State diagram

The transition to clogging is explored here for a given number of particles (N = 38), and a given applied flow strength. The state diagram indicating whether or not particles clog the channel as a function of the aperture D/d and the particle rigidity K is shown in Fig. 4.10. Two regions are well distinguished: a region at the upper-left corner where particles systematically form permanent clogs (red), and a region at the bottom-right corner where capsules never clog the channel (green). Symbols represent different initial positions of the particles. The transition between the clog and no clog states is weakly sensitive to the initial position of the particle, although both states may coexist close to the boundary between both domains (dashed black line). However, far from this boundary, the clog and no clog states are uniquely determined with confidence. Capsules with small rigidity (K < 200) pass freely without clogging the constriction, even if they have almost the same size as the channel width (D/d → 1). For the less deformable capsules considered in this study (K = 3000), the transition occurs at D/d = 2.75, which is close to the critical value D/d = 3 measured experimentally by Marin et al. [START_REF] Marin | Clogging in constricted suspension flows[END_REF] for rigid spherical particles. The slight deviation may be attributed to the difference in the constriction geometry. Marin et al. have used a progressive constriction, which promotes the formation of clogs upstream from the neck, while, here, a sudden abrupt constriction is considered. Nevertheless, the threshold seems to tend to the asymptotical value of D/d = 3 determined by Marin et al. as K → ∞. This state diagram reports the effect of the deformability of multiple fluid-filled capsules, and shows that soft particles can pass easily through the constriction, and by decreasing their elasticity, they are more likely to form clogs, especially at small apertures.

Dynamics of the particle passage

The number of escapees N (t), defined as the number of particles that have already passed through the constriction at time t, is given in Fig. 4.11 for stiff particles with Ca = 0.001 (left panel), and for soft particles with Ca = 0.1 (right panel). Three aperture degrees are considered, D/d = 1.25, 2.25, and 3.25.

Stiff particles rapidly form permanent clogs at small apertures, and, consequently, the number of escapees rapidly reaches a plateau, whose value is generally low. This is shown by the red curve in Fig. 4.11a for Ca = 0.001 and D/d = 1.25, which plateaus at a value of 2. At either high aperture or high capsule deformability, N (t) evolves approximately linearly at early stage, meaning the flux of particles dN/dt is almost constant. It increases with both D/d and Ca. Indeed, on the one hand, when the constriction width is enlarged, more particles can pass at the same time, which enhances the particle flux across the constriction. On the other hand, when Ca is increased, particles are more likely to deform and flatten so that they can "cooperate" and pass side by side through the neck, even at small apertures. Furthermore, they can adopt a more hydrodynamic shape that minimizes drag and reduces the residence time within the constriction. In the second stage, the number of escapees raises relatively slowly, and in a nonlinear manner. This stage corresponds to the passage of the remaining particles that were trapped in the corners and pushed against the lateral walls by the other particles, and to the particles that were initially located near the channel walls whose velocity is lower. When Ca and D/d are large enough, all the particles can pass easily through the constriction, and there is a total evacuation. This is shown in Fig. 4.11b for Ca = 0.1 and D/d = 2.25 and 3.25, where the number of escapees N (t) reaches a plateau at 38, corresponding to the total number of capsules contained in the channel. The time at which the last particle enters the outlet compartment is called "total evacuation time" and is noted T evac in what follows. It is reported in Fig. 4.12 as a function of D/d and Ca in logarithmic scale. It is found to decay when increasing the aperture D/d and the capillary number Ca, i.e. as the particle flux is enhanced. The evacuation time depends nonlinearly on the aperture and the capillary number. By fitting the data of Fig. 4.12, it has been possible to derive the following scaling laws giving the dependency of T evac on both D/d and Ca,

T evac ∝ D d -0.827
and T evac ∝ Ca -0.108 . (4.9)

The absolute value of the exponent related to Ca is relatively low compared to the one obtained for D/d, meaning the constriction width is the leading parameter for the evacuation time. The exponent of d/D (0.827) is larger than 1/2 known for the flow of granular dry particles through a pore. This might be due to the contribution of the particle deformation and the presence of a suspending fluid in the present study.

Conclusions

Multiple deformable capsules passing through a microfluidic constriction exhibit mainly the same dynamics as observed for a single capsule with D/d < 1 reported in Ref. [START_REF] Kusters | Forced transport of deformable containers through narrow constrictions[END_REF], with slight differences due to the collective motion. The transition from the blockage state to the passage is expected to depend on the applied flow strength. However, for a given Reynolds number, it is also mainly controlled by the ratio of the constriction width to the particle diameter, and the membrane deformability. Capsules either pass the constriction or get stuck by building a stable arch at the neck entrance. The coexistence of both states, as well as transient clogging, were observed for some couples of (D/d, Ca) that are close to the transition between the "clog" and "no clog" states. Nevertheless, for values of D/d and Ca far from the transition, the state is uniquely determined, independently of the initial random position of the particles. Capsules with extremely weak deformability were found to have a threshold that is very close to the one reported by Marin et al. [START_REF] Marin | Clogging in constricted suspension flows[END_REF] for rigid particles (D/d = 3). The established state diagram may serve in designing microfluidic devices for handling high throughput of soft particle suspension while preventing the formation of clogs. For instance, it can be used for mechanically-triggered drug delivery applications in order to enhance solute release from soft particles by determining the optimal aperture degree preventing clogs.

The dynamics of capsules passage through the constriction has also been carefully analyzed. The flux of particles crossing the constriction is roughly constant at early stage, and then it gradually decreases before all the particles have been evacuated, if no clog is formed. At large Ca and narrow constrictions (D/d → 1), particles cross the neck one by one, by forming a regular train, which then diffuses as a jet in the outlet compartment. This train can be exploited, for example, to load capsules in solute or to trigger local drug release thanks to the increased membrane permeability induced by the deformation in the narrow passage. It can also be used for diagnoses purposes. The evacuation times are found to depend nonlinearly on both D/d and Ca, although the effect of the particles deformability is small compared to the one of the aperture degree.

This study has examined the role of both the constriction confinement and the particle deformability for a given number of particles, while further exploring the parameter space is left for a future work, by varying, for example, the number of particles and the flow strength. The small number of particles used in the present study does not lead to a strong stochastic behavior as observed experimentally for a large number of particles [START_REF] Marin | Clogging in constricted suspension flows[END_REF][START_REF] Harth | Intermittent flow and transient congestions of soft spheres passing narrow orifices[END_REF], whose passage dynamics is described by statistics and probabilistic models.

General conclusion

Along this PhD thesis, we studied numerically mass transfer from core-shell capsules and fibers under various flow conditions relevant to many biomedical applications. Until now, most studies dealing with mass transfer from particles under flow consider constant and uniform surface boundary conditions, most of the time, for simplification reasons. However, such boundary conditions are not well adapted to the case of solute release during which the concentration inside the particles is not uniform and decreases over time. Here, we complement the existing studies by considering unsteady and nonuniform boundary conditions that appropriately model solute release.

We studied the contribution of capsules membrane permeability to mass transfer, first in absence of the flow, using the finite difference method. We showed that the release kinetics mainly depends on the membrane permeability and we developed tools to help designing future capsules. We also examined the effect of the flow on mass transfer from core-shell cylindrical fibers and spherical capsules using lattice Boltzmann simulations. The flow is found to enhance mass transfer by forced convection, which is quantified by correlations for the Sherwood number. The differences between considering a reservoir with unsteady and nonuniform boundary conditions and assuming constant surface concentration or mass flux have also been highlighted.

The reported results give a better fundamental understanding of mass transfer from reservoir particles by characterizing the effects of both the flow and the membrane permeability on the mass transfer efficiency. They may help in improving the performance of future particles for pharmaceutical and bioengineering applications, such as controlled drug release and cell culture.

However, the major part of this PhD thesis considers mass transfer from a single isolated reservoir particle. Modeling solute release from multiple soft particles is an ambitious problem and still remains an open issue. Here, the case of multiple capsules flowing into a microfluidic constriction has been examined, but the mass transfer part has not been studied. The collective motion of the capsules is expected to affect the mass transfer efficiency by creating local variations of the velocity field, but further works are needed to better understand solute release from multiple soft particles in complex flows.

List of publications and conferences

Well-stirred boundary condition

The well-stirred boundary condition set at the surface of the capsule is modeled by,

∂c ∂t = - 3R 2 R 3 b -R 3 D m ∂c ∂r r=R , (A.15)
where R is the radius of the capsule, R b is the radius of the bulk, and D m the membrane solute diffusion coefficient. This equation results from mass balance, expressing the fact that the mass of solute that leaves the capsule is equal to the mass of solute that enters into the bulk solution, whose concentration is assumed to be uniform due to stirring. Conversely, the mass of solute entering into the capsule is equal to the mass of solute leaving the bulk solution. This condition reads, Simplifying and rearranging terms finally give the desired result, .20)

∂c ∂t = - 3R 2 R 3 b -R 3 D m ∂c ∂r r=R . ( A

Figure 2 :

 2 Figure 2: (a) Different routes of drug administration (image from Tibbitt et al. [1]). (b) Schematic curves showing the drug concentration evolution over time for both conventional (uncontrolled) and controlled drug delivery.

Flow

Figure 3 :

 3 Figure 3: Microwell-based pancreas-on-a-chip designed at BMBI laboratory [4]. (a) Microfluidic chip, where the flow is from left to right. (b) Zoom on the microwell structures.

Figure 4 :

 4 Figure 4: (a) Blood oxygenator Maquet Quadrox based on fiber technology. The arrows show the directions of blood and oxygen flows. (b) Zoom in the fiber bundles, where O 2 flows within the lumen, and blood flows in the spaces between the fibers. Reproduced from Potkay et al. [7].

Figure 1 . 1 :

 11 Figure 1.1: (a) Capsule fabrication using flow-focusing technique [35]. (b) Multilayer capsule composed of an alginate core, an inner layer of N-isopropylacrylamide (NIPA) and an outer layer of N,N -dimethylacrylamide (DMAA) [41].

Figure 1 . 2 :

 12 Figure 1.2: Schematic diagram of a core-shell capsule immersed in a bulk fluid, with R c the radius of the core, R the external radius of the capsule, and R b the radius of the bulk medium. The shell is supposed to be homogeneous with a uniform thickness δ.

Figure 1 . 3 :

 13 Figure 1.3: Scheme showing the spatial discretization in the vicinity of the interface χ i located between the two adjacent layers i and i + 1. The inner points are represented as closed circles, and the point at the interface is denoted by an open circle. The space step ∆r is uniform in the whole domain.

Figure 1 . 5 :

 15 Figure 1.5: The D3Q19 lattice with its 19 discrete velocity directions.

Figure 1 . 6 :

 16 Figure 1.6: Concentration profiles at various times computed by the FD (solid lines) and the LB (squares) methods. Both solutions are in excellent agreement. The gray rectangle represents the location of the shell.

Figure 1 . 7 :

 17 Figure 1.7: Concentration profiles at various times obtained for a homogeneous sphere (a) and two capsules with membrane permeabilities p = 16×10 -8 m/s (b) and p = 4×10 -8 m/s (c). The membrane of the capsules has a thickness δ = 0.5 mm and is represented by the gray area. The gray solid line gives the solute diffusion coefficient.

Figure 1 . 9 :

 19 Figure 1.9: Percentage of solute released from a capsule (green) and a homogeneous sphere (red) in log-log scale (a), and as function of √ t (b).

Figure 1 . 11 :

 111 Figure 1.11: Comparison of the log-logistic (solid line), Weibull (dashed line) and firstorder (dotted line) equations. The log-logistic model leads to the best fit of the release curve computed by the FD solver (square symbols). The radii ratio is R/R c = 2 and the diffusivities ratio is D m /D b = 0.2. The values of the fitting parameters are given in Tab. 1.1.

Figure 1 . 12 : 2 :

 1122 Figure 1.12: Timescale parameter a LL (a) and shape factor b LL (b) as a function of D m /D b and R/R c . The black dots are the averaged values extracted from numerical simulations, and the colored surfaces are their fits computed by Eqs. (1.32a) and (1.32b), respectively.

Figure 1 .

 1 Figure 1.13: (a) Comparison between the release curves computed by the proposed equation and the semi-analytical solution of Ref. [45], with R c = 1.5 mm, R = 1.7 mm, D b = 3×10 -10 m 2 /s and D m = 0.5×10 -10 m 2 /s. (b) Comparison of the release kinetics predicted by the present equation with the one obtained by a LB simulation, with R c = 1 mm, R = 1.5 mm, D b = 2 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s (b). The release curves predicted by the present model compare very well with the expected solutions.

Figure 1 . 14 :

 114 Figure 1.14: Snapshots of the concentration field at various times computed by the lattice Boltzmann method. The solid lines delimit the core and the shell of the capsule. The indicated times are rounded down to the nearest hour.

Figure 1 .

 1 15b reports the RMSE measuring how much the FD solutions deviate from the LB solution. The red solid line gives the value of the membrane diffusivity D m minimizing the RMSE at any given core diffusivity D c . The global minimum of the RMSE (symbol ⊗) is obtained at D c = 2.1 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s. This estimated value of the membrane diffusivity is exactly equal to the value set in the LB simulation. The proposed method, thus, recovers exactly the expected membrane permeability p = 1.2 × 10 -7 m/s. The line D m = D c corresponds to the case of homogeneous spheres, as classically assumed in literature. The local minimum of the RMSE on this line gives the effective diffusivity D for the whole capsule, which here, is located at D m = D c = D = 1.1 × 10 -10 m 2 /s (symbol ⊕). The RMSE is strongly sensitive to changes in the membrane diffusivity D m , and it is weakly affected by variations of the core diffusivity D c , especially at low values of D m .

Figure 1 .

 1 Figure 1.15: (a) Release kinetics of insulin from a core-shell capsule. Squares: release kinetics computed by the LB method and used as input data for the characterization. Green line: best fit of the LB data when considering a capsule; Red line: best fit of the LB data when assuming the capsule is a homogeneous sphere. (b) The RMSE between the input data points and the FD solutions over a wide range of the fitting parameters D c and D m . The red solid line gives the value of the membrane diffusivity D m which minimizes the RMSE at any given core diffusivity D c . The absolute minimum of the RMSE is obtained at D c = 2.1 × 10 -10 m 2 /s and D m = 0.6 × 10 -10 m 2 /s (⊗), which match the values used in the LB simulation. They fit better the data compared with the solution obtained when considering the capsule as a homogeneous sphere with D = D m = D c = 1.1 × 10 -10 m 2 /s (⊕). The core radius, the membrane thickness and the bulk size are R c = 1 mm, δ = 0.5 mm, and R b = 12 mm, respectively.

Figure 1 .

 1 Figure 1.16: (a) Release kinetics of glucose from a hydrogel alginate capsule of size R = 1.73 mm and membrane thickness δ = 50 µm. Square symbols correspond to the experimental data of Rolland et al. [25], whereas solid lines correspond to the fits obtained with the proposed method. (b) The RMSE error between the FD solution and the experimental data, while varying the solute diffusion coefficient in both the core and the membrane. The global minimum of RMSE (⊗) gives the diffusion coefficients D c = 17.11 × 10 -10 m 2 /s and D m = 7.54 × 10 -10 m 2 /s, which fit better the experimental data than the effective diffusion coefficient D = 13.32 × 10 -10 m 2 /s (⊕) obtained assuming the capsule as a homogeneous sphere. The membrane permeability is then estimated to be p = 1.51 × 10 -5 m/s.

  . The global minimum (⊗) is obtained at D c = 76.40 × 10 -10 m 2 /s and D m = 2.96 × 10 -10 m 2 /s. The estimated membrane permeability of the capsule is then p = 1.85 × 10 -6 m/s. The local minimum along the line D m = D c (⊕), corresponding to having homogeneous spheres, is located at D m = D c = D = 7.98 × 10 -10 m 2

Figure 1 .

 1 Figure1.17: (a) Evolution of the normalized glucose concentration in the bulk over time. Squares: experimental data taken from Ref.[START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF] showing the glucose absorption kinetics for capsules made of alginate and polyethylene glycol, with core radius R c = 2.17 mm and membrane thickness δ = 160 µm. Solid line curves: fits obtained while considering the capsule as a homogeneous sphere (red), and with its real composite structure (green). (b) RMSE between the input experimental data and the FD solutions obtained over a wide range of the fitting parameters D c and D m . The local minimum along the line D m = D c (denoted by ⊕) corresponding to the case of a homogeneous sphere is obtained at D = 7.98 × 10 -10 m 2 /s, which is close to the value estimated by Koyama and Seki (D = 7.9×10 -10 m 2 /s). The global minimum of RMSE is located at D c = 76.40×10 -10 m 2 /s and D m = 2.96 × 10 -10 m 2 /s (symbol ⊗), leading to a better fit of the experimental data than the homogeneous sphere. The estimated membrane permeability is p = 1.85 × 10 -6 m/s.

  Figure1.17: (a) Evolution of the normalized glucose concentration in the bulk over time. Squares: experimental data taken from Ref.[START_REF] Koyama | Evaluation of mass-transfer characteristics in alginate-membrane liquid-core capsules prepared using polyethylene glycol[END_REF] showing the glucose absorption kinetics for capsules made of alginate and polyethylene glycol, with core radius R c = 2.17 mm and membrane thickness δ = 160 µm. Solid line curves: fits obtained while considering the capsule as a homogeneous sphere (red), and with its real composite structure (green). (b) RMSE between the input experimental data and the FD solutions obtained over a wide range of the fitting parameters D c and D m . The local minimum along the line D m = D c (denoted by ⊕) corresponding to the case of a homogeneous sphere is obtained at D = 7.98 × 10 -10 m 2 /s, which is close to the value estimated by Koyama and Seki (D = 7.9×10 -10 m 2 /s). The global minimum of RMSE is located at D c = 76.40×10 -10 m 2 /s and D m = 2.96 × 10 -10 m 2 /s (symbol ⊗), leading to a better fit of the experimental data than the homogeneous sphere. The estimated membrane permeability is p = 1.85 × 10 -6 m/s.

  For the mass transfer part, zero mass flux (Neumann) boundary condition ∂c/∂y| y=±W/2 = 0 is set along the channel walls (i.e. impermeable walls), zero concentration Dirichlet boundary condition c ∞ = 0 is set at the channel inlet, and zero gradient of the concentration ∂c/∂x| x=+L/2 = 0 is set at the channel outlet. Neither the concentration, nor the mass flux are constant at the surface of the reservoir. They are rather time-varying quantities that evolve in time when the reservoir is leaking its initial solute load. In absence of any interfacial resistance, continuity of both the concentration and the mass flux emerges naturally at the core-shell interface (r = R -δ), c core = c shell and D c ∂c ∂r core = D s ∂c ∂r shell , (2.6) and at the shell-fluid interface (r = R), c shell = c fluid and D shell ∂c ∂r shell = D f ∂c ∂r fluid , (2.7)

Figure 2 . 2 :

 22 Figure 2.2: The D2Q9 lattice used to compute the velocity and the concentration fields, with its 9 discrete velocity vectors e i .

. 24 )

 24 Different relaxation times are used to set different diffusion coefficients in each subdomain: τ Dc in the core, τ Ds in the shell and τ D f in the surrounding fluid (see Fig.2.1). The local concentration is computed by the zeroth-order moment of the distribution function g i , c(x, y, t) = 8 i=0 g i (x, y, t).(2.25)

. 32 )Figure 2 . 3 :

 3223 Figure 2.3: (a) Velocity profile computed by the LB flow solver (square symbols) and the expected analytical solution for the Poiseuille flow (green solid line). (b) Concentration profiles computed at dimensionless times T = 0.01, 0.02 and 0.04 with the LB method (red solid line) and with the analytical solution (green dashed line). The LB solution matches perfectly with the analytical one for both benchmark tests.

Figure 2 .

 2 4 shows the instantaneous average Sherwood number computed for the smallest value of the LBM relaxation time used in this study, τ Ds = 0.509. Both schemes give almost the same curves. Because the MRT is time consuming, most of the simulations are conducted with the BGK scheme in this study.

Figure 2 . 4 :

 24 Figure 2.4: Instantaneous average Sherwood number Sh(T) computed at various Reynolds numbers with both the BGK (solid line) and the MRT (dashed line) collision operators. The dimensionless solute shell permeability is P = 0.66 corresponding to the smallest value of the LBM relaxation time used in this study, τ Ds = 0.509. The Schmidt number is Sc = 1.

Figure 2 . 5 :

 25 Figure 2.5: Steady laminar flow past a cylindrical core-shell reservoir computed at Re = 1.

Figure 2 . 6 :

 26 Figure 2.6: Velocity field and streamlines obtained at Re = 80. The flow is steady and exhibits recirculations at the rear of the particle.

Figure 2 . 7 :

 27 Figure 2.7: Flow over a cylindrical core-shell reservoir computed for Re = 180 at T = 0.03, 0.06 and 0.09. The flow is unsteady and exhibits a von Kármán vortex street.

  The flow pattern observed by Sahin et al. at Re = 200 and B = 0.5 is very similar to the one reported here in Fig. 2.7 for Re = 180. The presently measured thresholds, Re c 1 = 21 and Re c 2 = 128, are in close agreement with the values given in Refs. [75, 87] at a blockage degree B = 0.5.

Figure 2 . 8 :

 28 Figure 2.8: (a) Normalized velocity profiles along the channel centerline computed for various Reynolds numbers covering the three typical flow patterns mentioned previously. (b) Zoom at the rear of the cylinder. The figures exhibit the main features of these three flow regimes, with regions of negative axial velocity indicating recirculations, and oscillations that are characteristic of unsteady flows with vortex shedding. From these figures, one can estimate that the first and second transition thresholds are Re c 1 = 21 and Re c 2 = 128, respectively. The black dotted lines represent the location of the cylinder.

Figure 2 . 9 :

 29 Figure 2.9: (a) Critical Reynolds number corresponding to the first transition at which a recirculation zone emerges at the rear of the cylinder as a function of the blockage degree (image taken from Chen et al. [87]). (b) Effect of the blockage ratio on the second transition, with the corresponding typical flow patterns. The second transition is much more complex than the first one, it is non-monotonic and shows bifurcations. The flow pattern observed by Sahin et al. at Re = 200 and B = 0.5 (figure 1) is similar to the one reported here in Fig. 2.7. Figure extracted from Sahin et al. [75].

Figure 2 . 10 :

 210 Figure 2.10: Flow streamlines and solute concentration field computed for a steady laminar flow at Re = 1, and with B = 0.5, Sc = 5 and P = 3.70. Mass transfer is mainly driven by diffusion, although the solute is slightly advected downstream by the flow.

Figure 2 . 11 :

 211 Figure 2.11: Solute concentration field around a cylinder whose surface is maintained at uniform and constant concentration c w = 1, in absence of the coating shell, computed at Re = 1, B = 0.5 and Sc = 5. The flow is the same as in Fig. 2.10.Mass transfer is also dominated by diffusion, and the solute spreads faster compared to the core-shell reservoir.

Figure 2 . 12 :

 212 Figure 2.12: Flow streamlines and solute concentration field computed at Re = 80, B = 0.5, Sc = 5 and P = 3.70. Steady recirculations develop at the rear of the core-shell cylinder, where most of the released solute is trapped before it skips downstream by forming a long plume-shaped concentration front. The overall concentration front does not reach the channel walls due to the increased efficacy of advection.

Figure 2 . 13 :

 213 Figure 2.13: Solute concentration field around a cylinder whose surface is maintained at uniform and constant concentration c w = 1, in cross-flow and in absence of the coating shell computed at Re = 80, B = 0.5 and Sc = 5.The background applied flow is the same as in Fig.2.12. The mass transfer reaches a steady regime, in contrast to the case of the core-shell reservoir reported in Fig.2.12 where the concentration decays over time.

Figure 2 . 14 :

 214 Figure 2.14: Flow streamlines and solute concentration field computed at Re = 180, B = 0.5, P = 3.70 and Sc = 5. The flow is unsteady and exhibits a von Kármán vortex street that advects the solute to the channel outlet in contrarotary vortices. The vortices trap a small amount of solute whose concentration decays later downstream.

Figure 2 . 15 :

 215 Figure 2.15: Solute concentration field around a cylinder whose surface is maintained at uniform and constant concentration c w = 1 at Re = 180, B = 0.5 and Sc = 5. The background applied flow is the same as in Fig. 2.14. The resulting

Figures 2 .Figure 2 . 16 :

 2216 Figure 2.16: Solute concentration c s (θ, T) at the surface of the cylindrical core-shell reservoir computed at various Reynolds numbers Re, while holding the same Schmidt number Sc = 5 and the same dimensionless shell solute permeability P = 3.70. Distortion of the concentration uniformity is largely amplified by increasing the flow speed.

Figure 2 .

 2 Figure 2.17 shows how the local mass flux ϕ(θ, T) is uniform only at the front of the reservoir and drops down to a minimum around θ = ± π 4 (the separation points) at larger Re. It gets larger values where the flow can advect the released solute downstream easily and can maintain a steeper concentration gradient.

Figure 2 . 17 :

 217 Figure 2.17: The local mass flux at the surface of the reservoir ϕ(θ, T) at different Reynolds numbers, while holding the same Schmidt number Sc = 5 and the same dimensionless shell solute permeability P = 3.70. ϕ(θ, T) is amplified by increasing Re. It is almost uniform, with a localized drop at the separation points when Re is large. ϕ(θ, T) is expressed in lattice units.

Figure 2 . 18 :

 218 Figure 2.18: Sh(θ, T) at the surface of the cylindrical core-shell reservoir in cross-flow at Re =1, 80 and 180 for Sc = 5 and P = 3.70.

Figure 2 . 19 :

 219 Figure 2.19: Sherwood number reported by Karniadakis [77] for a cylinder without a shell whose surface is sustained at constant flux. The Schmidt number is Sc = 0.7. Karniadakis considered heat transfer from the cylinder, and, here, Sh and Sc are respectively substituted for the Nusselt number Nu and the Prandtl number Pr.

Figure 2 . 20 :

 220 Figure 2.20: Evolution in time of the average Sherwood number Sh(T) at Re = 1, 80 and 180, corresponding to the three flow regimes. Sh(T) reaches a plateau at larger Reynolds numbers, whose value increases with the flow strength. The Schmidt number is Sc = 5 and the shell solute permeability is P = 3.70.

Figure 2 .

 2 Figure2.21 depicts the evolution in time of the released mass of solute R(T), Eq.(1.22), for the three flow regimes. The Schmidt number is Sc = 5 and the shell permeability is P = 3.70. The release curves exhibit the same characteristic features as observed in the first chapter for a spherical capsule in absence of flow. The solute first diffuses from the core to the shell. During this pure diffusion process, the solute is still protected by the shell and the flow has then no impact on the release. At T ≈ 0.015, the three curves detach from each other because the solute reaches the surface of the reservoir and starts being advected by the flow. Afterwards, R(T) evolves approximately as √ T for a short time (see Fig.2.21b). This follows the Higuchi model[START_REF] Higuchi | Rate of release of medicaments from ointment bases containing drugs in suspension[END_REF], which has originally been derived to describe drug release from particles in absence of external flow. However, this model approximates well the computed data at short-term only. At long-term, the release rate gradually slows down as the reservoir exhausts its content and R(T) deviates from the square root model. Higuchi's model cannot then be used to describe the whole release curves, as also discussed in the first chapter for a spherical capsule in absence of flow. The release is faster at larger Re due to the mass transfer enhancement by forced convection.

Figure 2 . 21 :

 221 Figure 2.21: Release rate of the solute from the cylindrical core-shell reservoir as a function of T (a), and √ T (b) computed at Re = 1, 80 and 180. The Schmidt number is Sc = 5 and the shell permeability to solute is P = 3.70.

  1.7), with discontinuities of the concentration radial derivative due to the boundary conditions Eqs. (2.6) and (2.7) set at the shell interfaces.

Figure 2 . 22 :

 222 Figure 2.22: Solute concentration profiles computed for P = 2.22, Re = 80 and Sc = 1. The location of the shell is represented by the gray rectangle. The solute diffusion coefficient is the same in the core and the fluid (D c = D f ), and it is smaller in the shell (D s < D f ).This leads to discontinuities of the concentration radial derivative at the shell interfaces, and to the slowdown of the release process.

Figure 2 . 23 :

 223 Figure 2.23: Steady average Sherwood number Sh as a function of P and the corresponding ratio D s /D f in linear (a) and log-log (b) scales. The gray area limit is not studied since it corresponds to the nonphysical case of D s > D f .

Figure 2 .

 2 Figure 2.24 reports the steady values of the average Sherwood number Sh versus the Reynolds number Re in a log-log scale for a core-shell cylinder with dimensionless shell solute permeability P = 3.70, and at Sc = 1. Data are obtained within the range 0.01 ≤ Re ≤ 180 covering the three aforementioned flow regimes.

Figure 2 . 24 :

 224 Figure 2.24: The steady average Sherwood number Sh versus the Reynolds number Re at different flow regimes, with Sc = 1 and P = 3.70. Sh is almost constant for Re ≤ 1, while it strongly increases with the flow strength for Re > 1.

Figure 2 . 25 :

 225 Figure 2.25: Sherwood number Sh vs. the Reynolds number Re (a), and the Schmidt number Sc (b) for different dimensionless shell solute permeabilities P. The computed data are fitted to extract the correlation power exponents m = 0.43 and n = 0.40, which visibly do not depend on P.For comparison purposes, data from simulations with constant surface concentration on a cylinder without a shell are also shown (down-pointing triangles), along with the correlations of Churchill-Bernstein[START_REF] Churchill | A correlating equation for forced convection from gases and liquids to a circular cylinder in crossflow[END_REF] and Hilpert[START_REF] Hilpert | Wärmeabgabe von geheizten Drähten und Rohren im Luftstrom[END_REF][START_REF] Bergman | Incropera's Principles of Heat and Mass Transfer[END_REF]. The presence of the shell shifts Sh upwards to larger values.

Figure 2 .

 2 25 suggests a correlation of the form: Sh = cRe 0.43 Sc 0.40 , (2.35)

Figure 2 . 26 :

 226 Figure 2.26: The correlation prefactor c = Sh/ (Re 0.43 Sc 0.40 ) vs. the dimensionless shell solute permeability P and the corresponding ratio D s /D f , computed for various combinations of the Reynolds and Schmidt numbers. Most data points collapse into a line of slope of l = -0.40.

Figure 2 .

 2 Figure 2.26 suggests the following correlation to estimate mass transfer from a single cylindrical core-shell reservoir placed in a cross-flow, Sh = 1.99P -0.40 Re 0.43 Sc 0.40 .(2.36)

  2.27.

Figure 2 . 27 :

 227 Figure 2.27: Sh as a function of the blockage ratio B for a cylindrical core-shell reservoir with dimensionless shell solute permeability P = 3.70 and at Sc = 5. The gray solid lines are fits of the numerical data with an equation of the form α exp (βB γ ), where the parameters depend on Re. Extrapolation of the data gives the expected values of Sh for an unconfined cylindrical reservoir (B = 0).

Figure 3 . 1 :

 31 Figure 3.1: Numerical setup of the studied problem. Subscripts "in" and "out" respectively refer to the properties of the inner encapsulated fluid, and the outer surrounding fluid. Continuity of the concentration and mass flux is considered on the capsule's membrane.

Figure 3 . 2 :

 32 Figure 3.2: Scheme illustrating the basis of the numerical method and the coupling between the fluid, mass and structure solvers.

Figure 3 . 3 :

 33 Figure 3.3: (a) Meshes of the particle obtained after 1, 2, 3 and 4 subdivisions of a regular icosahedron, corresponding respectively to icospheres with 80, 320, 1280 and 5120 triangular faces. (b) Spring network connecting the Lagrangian nodes of 6 adjacent triangular faces.

Fig 3 .

 3 5c. This figure shows periodic oscillations in time, meaning the capsule rotates at constant angular velocity ω = θ.

Figure 3 . 5 :

 35 Figure 3.5: (a) Snapshots of the capsule showing the position of the point P 0 that is initially located at r = (-R, 0, 0) during the first period (the indicated times correspond to the number of iterations divided by 100). (b) Normalized position of the point P 0 over time. (c) Instantaneous angular position θ of P 0 . The wall velocity is set to U = 0.05 (in lattice units).

Figure 3 . 6 :

 36 Figure 3.6: (a) The angular velocity ω of the capsule plotted versus the shear rate γ. The green solid line represents the analytical solution ω = 1 2 γ. Excellent agreement is achieved between the numerical simulations and the analytical solution. (b) Angular velocity ω as a function of the spring stiffness e s , with snapshots of the capsule obtained at e s = 0.1 (soft capsule) and e s = 7 (rigid capsule).

Figure 3 . 7 :

 37 Figure 3.7: Concentration profiles obtained at three different dimensionless times with the LB and FD methods, when the surrounding fluid is at rest (Re = 0). Both methods lead to the same solution.

Figure 3 . 8 :

 38 Figure 3.8: Streamlines in the vicinity of the capsule at Re = 1. (a) Three-dimensional representation of the streamlines. (b) Streamlines in the x-z plane passing by the center of the capsule. (c) Streamlines in the y-z plane exhibiting the spiraling nature of the flow induced by the spin of the particle. Same streamlines topology is observed over the range 0.01 ≤ Re ≤ 1.

4 and 1 ,

 1 and at dimensionless time T = 0.52. The Schmidt number is set to Sc = 10 for the three simulations. The relaxation times in the lattice Boltzmann equations (3.15) and (3.16) are τ ν = 1 and τ D = 0.55, respectively. These values are carefully chosen to ensure both the accuracy of the numerical solution and the stability of the lattice Boltzmann scheme. For Re = 0.01, solute transport is dominated by diffusion. The resulting concentration contours are concentric and centered with the capsule. At higher Reynolds numbers, advection becomes the dominant mass transfer mechanism and the solute is efficiently transported by the shear flow. Consequently, the concentration isocontours are stretched in the direction of the flow. The principle direction of diffusion is rotated towards the x-axis as the Reynolds number is increased, because the solute has less time to diffuse in the transverse direction.

Figure 3 . 9 :

 39 Figure 3.9: Snapshots of the concentration field computed at Re = 0.01, 0.4 and 1. The Schmidt number is Sc = 10 and the dimensionless time is T = 0.52.

Figure 3 . 10 :

 310 Figure 3.10: Local mass transfer quantities computed at Re = 0.01 (left), Re = 0.4 (middle) and Re = 1 (right). From top to bottom: normalized surface concentration c s /c 0 , surface mass flux ϕ and local Sherwood number Sh. The Schmidt number is Sc = 10 and the dimensionless time is T = 0.52.

Figure 3 . 11 :

 311 Figure 3.11: Instantaneous average Sherwood number at Re = 0.01, 0.4 and 1. The Sherwood number increases with the Reynolds number due to the mass transfer enhancement induced by forced convection and shear. The Schmidt number is Sc = 10.

5 2.438Pe 0.187 5 ≤ Pe ≤ 100 9 100 . ( 3 . 33 )•

 55100100333 + 0.580Pe 0.5 Pe ≤ -16.128Pe -0.349 Pe ≥ Correlation of Polyanin and Dil'man:

Figure 3 . 12 :

 312 Figure 3.12: Steady Sh vs. Re and Pe computed at Sc = 10 when considering continuous (circles) and Dirichlet (triangles) boundary conditions at the surface of the capsule. The correlations of Wang and Brasseur [26], Polyanin and Dil'man [92], and Longest and Kleinstreuer [94] are also shown for comparison purposes. The Sherwood number is dramatically altered by the type of boundary conditions set at the particle surface. The figure suggests Sh ∝ Re 0.59 for continuous boundary conditions and Sh ∝ Re 0.19 for Dirichlet boundary conditions in the advection-dominated regime (Pe > 1).

Figure 3 . 13 :

 313 Figure 3.13: Concentration field around the capsule computed when considering Dirichlet boundary conditions at its surface at Sc = 10. (a) Re = 0.01; (b) Re = 0.4; (c) Re = 1.The qualitative behavior is similar as the one reported for continuous boundary conditions in Fig.3.9, however, here, the overall concentration is higher since the surface of the capsule is sustained at constant concentration c 0 , leading to larger surface gradients.

Figure 4 . 1 :

 41 Figure 4.1: Scheme representing the microfluidic device designed by McFaul et al.[START_REF] Mcfaul | Cell separation based on size and deformability using microfluidic funnel ratchets[END_REF] to sort large sets of polydisperse particles (or cells) based on size and deformability. The device is constituted of multiple rows of constrictions whose width gradually reduces from bottom to top. Particles are transported by the flow toward the constrictions, and naturally, the bigger and less deformable particles will remain in the lower rows, while the smaller and more deformable particles will easily reach the upper rows of the device. The flow direction can be reversed if particles are blocked into the constrictions. In this way, by repeating flow inversions, particles will finally be sorted by size and deformability.

Figure 4 . 2 :

 42 Figure 4.2: Microfluidic device "Multiplex Fluidic Plunger" developed by Myrand-Lapierre et al.[START_REF] Myrand-Lapierre | Multiplexed fluidic plunger mechanism for the measurement of red blood cell deformability[END_REF] exploiting constrictions to characterize red blood cells (RBCs) deformability. This device can be used, for example, to diagnose diseases that rigidify RBCs, such as malaria or hemoglobinopathies. It consists of 34 tapered microchannels (only few shown here) where RBCs enter one by one and get deformed and stuck at the constriction entrance. Diseased RBCs being less deformable, they require a higher pressure difference P 2 -P 1 to pass to the post-constriction tube. They can then be identified and counted by the device.

Figure 4 . 3 :

 43 Figure 4.3: Experiment conducted at BMBI Laboratory by Le Goff et al.[START_REF] Goff | Squeezing bio-capsules into a constriction: deformation till break-up[END_REF] to characterize the breakup of biological and bioartificial capsules using a microfluidic constriction. The capsule is first progressively squeezed into the constriction due to the imposed pressure gradient between the inlet and the outlet of the channel. The pressure gradient is then increased step by step until the breakup of the capsule's membrane. Immediately after the breakup, the encapsulated solute is released and transported by the flow to the postconstriction chamber.

Figure 4 . 4 :

 44 Figure 4.4: Mechanically-triggered solute release from a soft capsule using a microfluidic constriction. Snapshots show the deformation of a capsule having tension-dependent permeability passing and squeezing through a narrow constriction at various times. The colors in the left panel indicate the local tension in the membrane, which is higher at the front of the capsule. The right panel shows the resulting membrane permeability and the concentration of the released solute (purple). The solute is mainly released at the front of the capsule, where the membrane permeability is larger. Image adapted from Luo et al. [31].

Figure 4 . 5 :

 45 Figure 4.5: Flow of a suspension of rigid particles in a microfluidic constriction at various neck-to-particle size ratios. Particles clog the channel at D/d = 1.02 and 2.07, in contrast to the case of D/d = 4 where no clog is observed, but the flow of particles is continuous despite the volume fraction is large (≈ 60%). Images are adapted from Marin et al. [32].

Figure 4 . 6 :

 46 Figure 4.6: Motion of a single capsule flowing through a narrow constriction. (a) Geometry of the constriction used by Kusters et al. [106]. (b) Phase diagram giving whether or not the particle can pass the constriction or gets stuck at its entrance depending on the confinement degree d/R 0 and the applied body force generating the flow, F m . (c) Phase diagram showing the passage and blockage states of the capsule as a function of the neck length L and the confinement degree d/R 0 . Adapted from Ref. [106].

Figure 4 . 7 :

 47 Figure 4.7: The numerical setup used to study the flow of capsules through a microfluidic constriction. The direction of the flow is from the left-lower pre-constriction chamber to the right-upper post-constriction chamber.

  Similar flow patterns are obtained for the other values of the aperture in the range 1.25 ≤ D/d ≤ 3.25, but with lower velocities in the constricted part.

Figure 4 . 8 :

 48 Figure 4.8: Flow field and streamlines in the constriction computed for Re = 0.1 and at D/d = 1.25, which is the smallest aperture considered in this study.

Figure 4 . 9 :

 49 Figure 4.9: Motion of multiple capsules (N = 38) in a microfluidic constriction for a given flow strength. The left panel shows the effect of the constriction aperture D/d, and the right panel the effect of particle deformability Ca. Easy passage of the particles is achieved at large D/d and Ca. Arch structures form at the entrance of the constriction for small apertures and rigid particles. At large Ca and narrow constrictions, the particles cross the neck one by one forming a regular train, and they exit the constriction as a jet.

Figure 4 .

 4 Figure 4.10: The "clog" and "no clog" states of multiple capsules (N = 38) obtained when varying their deformability Ca and the constriction width D/d, while holding the same applied flow strength. The deformability character of the particles reduces the threshold of the transition to the no clog state below the value D/d = 3 measured experimentally for rigid particles by Marin et al. [32]. Simulations with different initial conditions are represented with different symbols.

Figure 4 . 11 :

 411 Figure 4.11: Dynamics of the number of escapees N (t) for rigid particles with Ca = 0.001 (a) and for soft particles with Ca = 0.1 (b) at various constriction widths D/d. The indicated values of t correspond to 10,000 time iterations.

Figure 4 . 12 :

 412 Figure 4.12: Evacuation time T evac needed for 38 capsules to pass the constriction as a function of the capillary number Ca (a) and the aperture D/d (b). T evac is a nonlinear decreasing function of both Ca and D/d. The indicated values of T evac correspond to 10,000 time iterations.

Articles 1 .

 1 C. Bielinski, N. Le, B. Kaoui, Unsteady mass transfer from a core-shell cylinder in crossflow, Physical Review Fluids 6, 023501 (2021) 2. C. Bielinski, B. Kaoui, Characterisation of capsule membrane permeability, (Under review) 3. C. Bielinski, O. Aouane, J. Harting, B. Kaoui, Squeezing multiple soft particles into a constriction: transition to clogging, (Submitted) 4. C. Bielinski, L. Xia, B. Kaoui, Unsteady mass transfer from a spherical capsule under shear flow, (In preparation) Proceeding 1. C. Bielinski, O. Aouane, J. Harting, B. Kaoui, Squeezing multiple soft particles into a constriction: transition to clogging, Comptes-rendus de la 24e Rencontre du Non Linéaire Paris 2021 2. C. Bielinski, B. Kaoui, Transfert de masse à partir d'une particule noyau-coque sous écoulements stationnaires et non stationnaires, Comptes-rendus de la 23e Rencontre du Non Linéaire Paris 2020

18 )

 18 with M b and M the mass of solute in the surrounding bulk and the capsule, respectively. Because the external solution is well-stirred, the concentration of solute in the bulk c b is assumed to be uniform, and, thus,M b (t) = V b c b (t), where V b = 4 3 π (R 3 b -R 3) is the volume of the bulk, excluding the capsule. The mass of solute inside the capsule at time t is M (t) = 4π R 0 c(r, t)r 2 dr, which gives, Assuming the concentration is continuous at the shell/bulk interface, i.e., c b (R, t) = c(R, t),

Table 1 .

 1 and 1.12b, respectively. 1: Parameters and goodness of fit for the first-order, the Weibull and the loglogistic models used to fit the release curve computed by the FD solver in Fig.1.11. The log-logistic equation has the highest R 2 adj and the smallest RMSE. It produces the best fit of the numerical data.

	Model	Parameters	R 2 adj	RMSE
	First-order	a F = 1.1421	0.9851	0.0185
	Weibull	a W = 0.8733	0.9924	0.0132
		b W = 1.2395
	Log-logistic	a LL = 0.6376	0.9993	0.0040
		b LL = 2.0639

They are perfectly fitted (R 2 adj > 0.999) by the surfaces of equations,

Table 2 .

 2 1: Dependency of the computed steady average Sherwood number Sh on the grid resolution at the largest considered Reynolds number Re = 180, and at Sc = 5 and P = 3.70.

		Sh	Deviation (%)
	2401 × 481	20.69	-
	2001 × 401	20.54	0.72
	1601 × 321	20.30	1.88
	1001 × 201	19.87	3.96

Table 3 .

 3 1: Angular velocity of the capsule at various wall velocities U ranging between 0.01 and 0.05. The ratios ω/γ are very close to the expected theoretical value 1/2, which validates the implementation of the fluid-structure coupling. All quantities are expressed in lattice units.

	-4 ) ω/γ

(a) (b)

Acknowledgements

My PhD thesis has been a very rewarding experience that would not have been possible to achieve without the help and guidance I received from many people.

This appendix provides the derivation of Eqs. (1.35) and (1.40), which are classically used to estimate the effective solute diffusion coefficient of spherical capsules by fitting experimental solute release and absorption curves, respectively. Both equations are derived assuming the capsules as homogeneous spheres that are immersed into a well-stirred solution of limited volume, whose concentration is uniform due to stirring, but evolves in time. The derivation of Eq. (1.10) modeling the well-stirred condition is also given in this appendix.

Solute absorption

For the case of solute absorption, the concentration C t at the surface of the homogeneous spherical capsule at time t is,

where C 0 is the initial concentration in the bulk, R is the radius of the capsule, D is the solute diffusion coefficient, and α = V b /V is the ratio of the volumes of the bulk and the capsule. The q n 's are the nonzero positive roots of tan q n = 3q n /(3 + αq 2 n ). The starting point of the derivation of Eq. (A.1) is the analytical concentration profile within the capsule given in the book by Crank [START_REF] Crank | The Mathematics of Diffusion[END_REF] The Mathematics of Diffusion by Eq. (6.36),

with C eq the expected concentration at equilibrium, i.e., when t → ∞. The surface concentration C t is obtained by evaluating Eq. (A.2) at r = R,

Mass conservation gives the following relationship between C eq and C 0 ,

Inserting this expression of C eq in Eq. (A.3) and dividing by the initial concentration C 0 yields the final result,

Solute release

For the case of solute release from a homogeneous spherical capsule to the bulk solution, the concentration at the surface of the sphere is given by,

with C eq the expected concentration as t → ∞. The derivation of this equation is very similar to the one of the absorption case. The concentration profile within the sphere is obtained by making the following substitutions in Eq. (A.2), see Eq. (6.38) in the book by Crank [START_REF] Crank | The Mathematics of Diffusion[END_REF],

where, here, C 0 is rather the initial concentration inside the sphere. These substitutions lead to,

Evaluating this equation on the surface of the sphere (r = R) yields,

which leads to,

Inserting the above expression of C 0 in Eq. (A.10) gives,

Simplifying and rearranging terms lead to the desired result,

Appendix B: Sherwood number in the limit of Re = 0

This appendix gives a proof of the classical result Sh 0 = 2 for the case of mass transfer from a sphere of radius R sustained at constant concentration c s , and that is immersed into a stagnant fluid (Re = 0).

Mass transfer from a spherical particle to a fluid at rest is governed by the diffusion equation, This proves that the Sherwood number for a sphere maintained at constant concentration and immersed into a stagnant fluid is Sh 0 = 2. It should be noted that some authors use the radius of the particle as the characteristic length in the definition of Sh, instead of its diameter. With this alternative definition, the Sherwood number in the limit of Re = 0 is rather Sh 0 = 1.