Paul Breugem

Modelling and Simulations of Contacts in Particle-Laden Flows

Keywords: Local Lubrication Models, Fluid-Structure Interactions, Coupled VP-DEM, Ellipsoidal Particles

Particle-laden flows can be found in many industrial applications such as slurry transport or the chemical industry in general. In mixtures made of solid particles emerged in a viscous fluid, particle interactions play an essential role in the overall mixture viscosity.

The suspension phenomenon is caused by short-range hydrodynamic interactions, known as lubrication. Lubrication forces are usually underestimated due to their singularities and the spatial discretization of the numerical schemes. In this thesis, we propose a lubrication model for a coupled volume penalization method and discrete element method solver that estimates the unresolved hydrodynamic forces and torques in incompressible Navier-Stokes flows. Corrections are made locally on the surfaces of the interacting particles without any assumption on the global particle shapes. The final version of the local lubrication model can be used for suspension of convex particles without any tabulations.

The numerical method has been validated against experimental data with spherical and ellipsoidal particles. With spherical particles, the lubrication model performs as well as existing numerical models that are limited to this specific particle shape. The model compatibility with convex particles has been validated by comparing simulations using ellipsoids to experimental measurements we made.

Modélisation et simulations numériques des contacts dans des écoulements chargés en particules

Résumé :

Les écoulements chargés en particules sont présents dans de nombreuses applications industrielles telles que le transport de boues ou l'industrie chimique en général. Dans des mélanges constitués de particules solides immergées dans un fluide visqueux, les interactions entre particules jouent un rôle essentiel dans la viscosité globale du mélange.

Le phénomène de suspension est causé par des interactions hydrodynamiques à courte distance, connues sous le nom de lubrification. Les forces de lubrification sont généralement sous-estimées en raison de leur nature et de la discrétisation spatiale du problème.

Dans cette thèse, nous proposons un modèle de lubrification qui estime les forces et couples hydrodynamiques non résolues par un solveur couplant la résolution des équations de Navier-Stokes incompressible par une méthode de volumes penalisés, à la résolution de la dynamique des particules par une méthode aux éléments discrets. Les corrections des contraintes hydrodynamiques sont faites localement sur la surface des particules en interaction sans aucune hypothèse sur la forme générale des particules. La version finale du modèle de lubrification proposée peut être utilisée pour des suspensions de particules convexes sans aucune tabulation. La méthode numérique a été validée avec des particules sphériques et des ellipsoïdes, en comparant des simulations à des données expérimentales.

Dans le cas de particules sphériques, le modèle de lubrification est aussi précis que les modèles de lubrification existants qui sont limités à ce type de géométrie. La compatibilité du modèle avec des particules convexes a été validée en comparant des simulations, utilisant des ellipsoïdes, à des mesures expérimentales que nous avons réalisées.

Mots clés :

Modèle local de lubrification, Interactions fluide-structure, Couplage VP-DEM, Ellipsoïdes.

Résumé

L'étude d'écoulements de particules a de nombreux intérêts pour des applications industrielles et de recherche, telles que la fabrication de béton, le transport de boues ou de sédiments, le traitement des eaux usées, ou la fabrication de certains plastiques.

La simulation numérique d'écoulement de particules solides immergées dans un fluide visqueux apporte des informations importantes pour la compréhension des phénomènes physiques en jeu, ou peut servir de support pour l'optimisation de procédés industriels.

Depuis le début des années soixante, de nombreuses méthodes numériques ont été développées afin de simuler précisément ce type d'écoulement. Les méthodes les plus avancées simulent le mouvement des particules ainsi que l'écoulement du solvant. L'approche la plus courante pour résoudre l'écoulement est de coupler un solveur fluide, résolvant les équations de Navier-Stokes, avec un solveur des équations de Newton-Euler résolvant les mouvements des particules. Bien que ces méthodes permettent de résoudre finement l'écoulement, elles ne sont pas capables de résoudre intégralement les interactions hydrodynamiques entre des particules pratiquement en contact, et cela indépendamment de la nature du maillage choisi (maillage adapté ou non aux interfaces fluides-particules).

Quand une particule de rayon a est sur le point d'entrer en collision avec un obstacle, la présence de l'obstacle est perçue par la particule par l'intermédiaire du solvant par la modification des forces hydrodynamiques. L'apparition d'interactions hydrodynamiques de courtes distances, communément appelées forces de lubrification, est la conséquence du drainage du fluide dans la zone interstitielle entre la particule et l'obstacle. La composante normale et tangentielle des forces de lubrification évoluent comme -1 et log( ) respectivement, avec a la distance de séparation entre la particule et l'obstacle. La résolution numérique de ces singularités n'est possible que pour certains écoulements de Stokes, et a un coût de calcul important. Dans le cas général, les forces de lubrification sont résolues partiellement. C'est pour cela que des modèles de lubrification, utilisant des solutions analytiques provenant de la théorie de la lubrification, sont ajoutés aux méthodes numériques.

Dans cette thèse, nous avons développé un modèle de lubrification pour des écoule-RÉSUMÉ ments de particules rigides immergées dans un fluide newtonien visqueux. Alors que la majorité des modèles de lubrification existants sont limités à des particules sphériques ou à des écoulements en régime de Stokes, notre modèle peut s'adapter à des écoulements de particules convexes en régime inertiel incompressible. Afin de valider le modèle, ce dernier a été intégré à un solveur partitionné couplanet la résolution des équations de Navier-Stokes incompressible par une méthode de pénalisation (VP), à la résolution du mouvement des particules par la méthode aux éléments discrets (DEM) (voir Chap. 2).

La méthode numérique.

L'écoulement est modélisé par les équations de Navier Stokes incompressible pour un fluide newtonien de viscosité µ et de densité ρ :

           ∂u i ∂t + ∂ (u i u j ) ∂x j = - 1 ρ ∂p ∂x i + µ ρ ∂ ∂x j ∂u i ∂x j +χλ (u τ,i -u i ) , ∂u i ∂x i = 0, (1)
avec (i; j) = {1, 2, 3} 2 , u i étant les composantes de la vitesse, et p la pression. Ce système d'équations est résolu sur un maillage cartésien uniforme couvrant la totalité du domaine d'intérêt. Les cellules solides de maillage (à l'intérieur des particules) sont repérées à l'aide d'une indicatrice χ qui renvoie le signe de la level-set globale. La level-set est une fonction signée qui donne en chaque élément du maillage la distance entre cet élément et la surface de la particule la plus proche, et cette distance est comptée négativement quand l'élément est à l'intérieur de la particule.

A l'interface fluide-particule, on considère une condition de non-glissement qui est implicitement imposée par le therme de pénalisation χλ (u τ,i -u i ), avec le coefficient de pénalisation λ choisi arbitrairement grand [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF] . Les vitesses pénalisées aux cellules solides à l'interface fluide-particule sont corrigées par la méthode IPC (Image Point Correction) [START_REF] Hovnanian | An accurate cartesian method for incompressible flows with moving boundaries[END_REF] afin d'obtenir une pénalisation globale consistante et de second ordre en espace. La méthode IPC consiste à corriger les vitesses pénalisées de manière à ce que les vitesses ainsi que leurs dérivées satisfassent la condition de non-glissement à l'interface.

Les équations de Navier-Stokes sont discrétisées en espace en utilisant un arrangement colocalisé au centre des mailles des variables primitives (p, u i ). L'intégration en temps se fait à partir du schéma de projection classique de Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] et Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF] (voir Sec. 2.2).

La dynamique des particules est calculée à partir des équations de Newton-Euler (voir RÉSUMÉ Sec. 2.3) appliquées au centre de masse des particules :

     m i dU i dt = F i , dJ i Ω i dt = T i , (2)
avec m i la masse, J i la matrice d'inertie, U i et Ω i les vitesses de translation et de rotation de la particule P i . La force F i et le couple T i sont la force et le couple résultant des collisions, des effets hydrodynamiques et de la gravité appliqués à la particule P i . Les forces et couples générés par la collision de la particule avec un obstacle (mur ou autre particule) sont modélisés par le modèle de sphères molles [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] (voir Sec. 2.3.2). Les forces et couples hydrodynamiques sont décomposés en deux contributions : Le modèle de correction local des forces de lubrification.

     F hyd i = F solv i + F deg i , T hyd i = T solv i + T deg i . (3)
Le modèle de correction des forces et couples de lubrification que nous avons développé est fondé sur les mêmes hypothèses de base que la théorie de lubrification pour des particules sphériques. On considère deux particules isolées dans un espace infini. Le fluide est considéré au repos ainsi qu'une des particules, alors que la seconde particule est animée d'un mouvement constant de translation en direction de la particule au repos. Ce mouvement entraine la compression du fluide entre les deux particules ce qui est à l'origine de la composante principale de la force de lubrification. On considère aussi, dans un second temps, la translation d'une particule à proximité d'une seconde particule initialement au repos. Ce mouvement entraine le cisaillement du fluide interstitiel ce qui génère une force et un couple hydrodynamiques sur les deux particules.

RÉSUMÉ

Si on suppose que la distance entre les deux particules est très petite, le mouvement du fluide dans la zone interstitielle peut être approché par les équations de Stokes. Ainsi il est possible, avec un peu de calcul, d'exprimer la force et le couple hydrodynamiques agissant sur un voisinage du point de contact à la particule initialement en mouvement.

Comme les particules sont isolées, toutes les forces hydrodynamiques du système sont des forces de lubrification (voir Sec. 3.1, pour la théorie de lubrification décrite en détails). On obtient donc ainsi les expressions des forces et couples de lubrification pour un écoulement en régime de Stokes.

Le modèle utilise l'ordre dominant des forces et couples théoriques, comme correction des forces et couples hydrodynamiques, de la manière suivante :

       F hyd i = F solv i + j∈ 1,N \{i}
F lub i,j + F lub i,wall ,

T hyd i = T solv i + j∈ 1,N \{i} T lub i,j + T lub i,wall , (4)
avec F lub i,j et T lub i,j l'ordre dominant des force et couple de lubrification entre les particules P i et P j définis par la théorie de lubrification (voir Sec. 3.3 pour les expressions exactes).

Ces corrections peuvent être appliquées directement pour des particules sphériques. Dans le cas de particules convexes, les corrections sont calculées à partir de sphères virtuelles.

Pour chaque couple de particules en interaction, deux sphères sont construites de manière à ce que leurs surfaces approchent au mieux les surfaces des particules au voisinage du point de contact.

Validations du modèle de lubrification.

Le modèle de lubrification a été validé, dans un premier temps, avec des particules sphériques (voir Chap. 4) pour lesquelles il existe une large bibliographie de mesures expérimentales. Dans le cas le plus simple d'une unique particule impactant un mur, le modèle local de lubrification est aussi précis que les modèles de lubrification existants.

Cependant, notre modèle est en moyenne 10% plus couteux que les modèles de référence spécifiquement conçus pour des suspensions de sphères. En revanche, le modèle local n'est pas limité à des particules sphériques et ne nécessite pas la tabulation de paramètres.

La validation pour des particules convexes est intrinsèquement plus difficile que pour des sphères (voir Chap. 5). Dans beaucoup de cas, les rotations et couples peuvent être négligés pour des particules sphériques ce qui n'est pas possible pour des ellipsoïdes. L'étude expérimentale des interactions entre un ellipsoïde avec un obstacle est à ce jour très peu documentée. Avec le soutien de TU-Delft, nous avons pu effectuer nos propres RÉSUMÉ mesures d'un ellipsoïde impactant un mur. La comparaison de trajectoires simulées à ces premières mesures expérimentales a permis de donner une première estimation de la qualité des corrections des forces de lubrification pour des particules non sphériques.

Introduction

Mixtures of solid particles in a liquid, i.e. suspensions, are a type of two-phase flow that appear in various industrial applications as well as in natural flows. For instance, sedimentation in rivers or near-coast influences the flow. By understanding the dynamics of the bed-load transport, the evolution of the river-stream can be forecast to prevent abnormal erosion or obstruction of waterways [START_REF] Sun | CFD-DEM simulations of current-induced dune formation and morphological evolution[END_REF][START_REF] Frohlich | Numerical simulation of sediment transport in open channel flow[END_REF] . Macroscopic suspensions are also relevant in several industrial applications, such as nuclear waste processing, water treatment [START_REF] Santos | A stochastic model for particulate suspension flow in porous media[END_REF][START_REF] Galaguz | Modeling of particle filtration in a porous medium with changing flow direction[END_REF] , slurry transportation [START_REF] Kuchiic | Cfd modeling for pipeline flow of fine particles at high concentration[END_REF] , reinforced plastics manufacturing, or the animation industry [START_REF] Wan | Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method[END_REF][START_REF] Delmotte | Hydrodynamics of suspensions of passive and active rigid particles: A rigid multiblob approach[END_REF] . For dense mixtures such as uncured solid rocket fuel or concretes [START_REF]Simulation of Fresh Concrete Flow[END_REF] , a high concentration of solid particles is desired without compromising the rheological properties and the flowing behavior of the mixture. Accurate numerical methods are then valuable supports to optimize manufacturing processes.

Over the years, several numerical methods have been developed to simulate particleladen flows at different scales. Particle-laden flows can be simulated at a macroscale with methods seeing the mixture as a non-Newtonian fluid. If the characterization of the mixture rheology is possible, these methods provide a general picture of the mixture flows at a reduced cost. However, the particle-fluid interactions need to be fully resolved to clearly describe the flow. Therefore, microscale methods are preferred as they fully resolve the dynamics of the particles as well as the flow around each of them.

Microscale methods have historically started with adaptations of molecular models, such as Stokesian Dynamics (SD) [START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF][START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] or Force-Coupling Method (FCM) [START_REF] Climent | Numerical simulations of random suspensions at finite Reynolds numbers[END_REF][START_REF] Maxey | Localized force representations for particles sedimenting in Stokes flow[END_REF] . Both methods depend on a truncated multipole expansion of the Stokes equations. Thereby, they are usually suited only for specific particle shapes and are inherently restricted to Stokes flows. Thanks to modern improvements of computer performances, Direct numerical simulation (DNS) has emerged as an appealing alternative. By solving the governing equations (Navier-Stokes equations) directly without any further assumptions, DNS enables microscale simulations of arbitrary particle shape at the numerical method accuracy level. The DNS methods can be divided into two classes. The first one considers boundary-fitted approach [START_REF] Joseph | Direct simulation of fluid particle motions[END_REF][START_REF] Johnson | 3D simulation of fluid-particle interactions with the number of particles reaching 100[END_REF] , where the fluid fills the meshed domain. If the mesh ele-INTRODUCTION ments are small enough, this approach enables an accurate computation of the boundary layers and interactions between particles. However, as the particles move, the domain needs to be constantly adapted leading to possibly complex and expansive remeshing issues. The computational cost of the re-meshing limits three-dimensional simulations to about a hundred particles [START_REF] Johnson | 3D simulation of fluid-particle interactions with the number of particles reaching 100[END_REF] . In contrast, the second class of DNS methods considers non-boundary-fitted meshes which are more suited to numerical simulations using many particles. In this approach, the whole domain is discretized on a Eulerian fixed grid where the particles are embedded. Several techniques are available to simulate the particulate flow using fictitious domain methods [START_REF] Martin | A numerical method for fully resolved simulation (FRS) of rigid particle-flow interactions in complex flows[END_REF][START_REF] Pan | Distributed Lagrange multiplier methods for incompressible viscous flow around moving rigid bodies[END_REF][START_REF] Sharma | A fast computation technique for the direct numerical simulation of rigid particulate flows[END_REF][START_REF] Simeonov | Modeling mechanical contact and lubrication in direct numerical simulations of colliding particles[END_REF][START_REF] Wachs | A DEM-DLM/FD method for direct numerical simulation of particulate flows: Sedimentation of polygonal isometric particles in a Newtonian fluid with collisions[END_REF][START_REF] Shao | A fictitious domain method for particulate flows with heat transfer[END_REF][START_REF] Merlet | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] , encompassing lattice Boltzmann methods [START_REF] Ladd | Lattice-Boltzmann simulations of particle-fluid suspensions[END_REF] , immersed boundary methods (IBM) [START_REF] Feng | Proteus: a direct forcing method in the simulations of particulate flow[END_REF][START_REF] Uhlmann | An immersed boundary method with direct forcing for the simulation of particulate flows[END_REF] and volume penalization method (VP) [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF] .

Although microscale methods can capture flow details, they are challenged by shortrange hydrodynamics. When two particles are moving toward each other they start to interact throughout the fluid as the separation distance a (where a denotes the particle radius and ≥ 0) becomes small. This well-known lubrication effect is due to the draining of interstitial fluid in the gap between the two interacting particles. The normal and tangential components of the lubrication force diverge as the particles collide ( tends to 0) as -1 and log( ), respectively. From a physical point of view, the lubrication forces do not reach infinitely large values as the interstitial flow collapses when the gap becomes smaller than roughly the fluid molecule size. The void generated in the gap enables the mechanical collision of the particles. This critical state can be reached when the particle inertia is non-negligible compared to the fluid viscous friction.

At vanishing Reynolds number, DNS methods are able to resolve these lubrication forces if the grid spacing is small enough. Typically, the grid spacing should be at least smaller than 10 -3 a to capture lubrication effects [START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF] . For weakly inertial flows, collisions between particles might occur making impossible to run DNS simulation without introducing a collision model and approximate the lubrication effects. The computational costs make long-term simulations of concentrated suspensions prohibitive. Most simulations are performed on meshes with a grid spacing of about 10 -1 a, meaning that the accuracy of the description of the lubrication effects drops as the particles come in near contact. Therefore, a numerical lubrication model is usually introduced to balance the unresolved lubrication forces.

In the simplest lubrication model, the theoretical lubrication force known for two spherical particles (normal [START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF] and tangent [START_REF] O'niell | On the slow motion of a sphere parallel to a nearby plane wall[END_REF] components) is added to the computed hydrodynamic force [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] . This approach considers that lubrication cannot be captured by the flow solver. However, by adding the theoretical lubrication force, the lubrication par-

INTRODUCTION

tially captured by the solver is counted twice. A more accurate approach is suggested in Stokesian Dynamics, where hydrodynamic forces are split into a long-and a short-range actions. The short-range actions are also known as lubrication. Hence long-range actions are provided by the solver and short-range actions are modeled via the lubrication theory, avoiding an overestimation of the total hydrodynamic forces. However, this decomposition is not easy to adapt to DNS since the long-range hydrodynamics is not explicitly known. A common technique [START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF][START_REF] Nguyen | A fictitious-domain simulation of solid-liquid flow with subgrid lubrication force correction; a sphere falling onto a plane surface[END_REF] is to perform off-line simulations of two isolated particles using a refined mesh without any lubrication model, for different configurations and normalized separation distance , to estimate the unresolved part of the lubrication forces. From these off-line simulations and the lubrication theory, a tabulation of lubrication corrections is created. Hence the lubrication model for on-line simulations estimates the correction to perform from the tabulation. The accuracy of this method depends on the quality of the tabulation and the diversity of the off-line samples used to create it.

Thus, these methods are usually used for mono-disperse suspensions of spherical particles since only the tabulation of a single parameter is required. Theoretically, the tabulation of more general suspension models is also possible [START_REF] Cox | The motion of suspended particles almost in contact[END_REF] . However, the computational cost involved to generate an accurate multi-parameter tabulation makes general suspensions unachievable in reasonable CPU time. This thesis is focused on lubrication methods applied to non-Brownian suspensions. This manuscript is organized as follows:

In the first chapter, several methods handling problems caused by unresolved lubrication are reviewed. The proposed brief review of the state of the Art aims to highlight the challenges and limitations facing the lubrication correction methods. This chapter also draws the context of this thesis research and places the local lubrication correction model, that we created, among the existing methods.

In the second chapter, the Lagrange-Euler method used in this manuscript to solve particle-laden flow is described in depth. Incompressible Navier-Stokes equations for a Newtonian viscous fluid have been solved using a volume penalization (VP) method coupled with a discrete element method (DEM) [START_REF] Radjai | Discrete-element Modeling of Granular Materials[END_REF] , which resolves the rigid particle dynamics.

In the third chapter, the local lubrication correction model is introduced. The completed theoretical background of the model is first detailed before describing the model implementation into the numerical framework introduced in the previous chapter.

In the last two chapters, the local lubrication correction model is tested with several cases. The model is first validated with spherical particles, where extensive benchmarks are available in the literature. The validation is then extended to ellipsoidal particles

INTRODUCTION

where we carried out our own measurements to complete scarce references.

Chapter 1

State of the Art

Lubrication correction techniques are closely linked to the numerical methods used to solve the particle-laden flows. Nevertheless, they are all based on the lubrication theory established in the first half of the XX th century. Stimson et al. [START_REF] Stimson | The motion of two spheres in a viscous fluid[END_REF] first exhibited analytical solutions of the lubrication forces generated by two particles moving toward each other in a viscous fluid at vanishing Reynolds number. The squeezing motion of a particle toward an obstacle (wall or another particle) gives the dominant and singular component of the lubrication forces. Therefore, this motion has been the focus of theoretical and early numerical studies [START_REF] Brenner | The slow motion of a sphere through a viscous fluid towards a plane surface[END_REF][START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF][START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF][START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-Reynolds-number flow[END_REF][START_REF] Jeffrey | Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow[END_REF][START_REF] Jeffrey | The calculation of the low Reynolds number resistance functions for two unequal spheres[END_REF] .

The second singularity of the lubrication forces is induced by the shearing motion of the interstitial fluid between a particle and a nearby obstacle. This motion is generated either by the translation [START_REF] O'niell | On the slow motion of a sphere parallel to a nearby plane wall[END_REF][START_REF] O'neill | Asymmetrical slow viscous fluid motions caused by the translation or rotation of two spheres. part ii: Asymptotic forms of the solutions when the minimum clearance between the spheres approaches zero[END_REF] or the rotation [START_REF] Cox | Slow viscous motion of a sphere parallel to a plane wall-ii couette flow[END_REF][START_REF] Cooley | On the slow rotation of a sphere about a diameter parallel to a nearby plane wall[END_REF] of the particle. In addition to the second component of the lubrication force, the shearing motion creates a lubrication torque.

Analytic solutions of the lubrication effects have widely been established for a system of two rigid spherical particles in a creeping flow. However, extensions of the lubrication theory can be found for spherical particle with rough surfaces [START_REF] Jenkins | Hydrodynamic interaction of rough spheres[END_REF] , deformable drops [START_REF] Serayssol | The elastohydrodynamic collision of two spheres[END_REF][START_REF] Schonberg | The lubrication force between two viscous drops[END_REF] , or non-spherical particles [START_REF] Cox | The motion of suspended particles almost in contact[END_REF] . Cox et al. [START_REF] Cox | The slow motion of a sphere through a viscous fluid towards a plane surface-ii small gap widths, including inertial effects[END_REF] extend the lubrication theory to unsteady Navier-Stokes flows of axisymmetric particles, widening significantly the Reynolds number range where the theoretical lubrication solutions are applicable. This chapter is focusing on the main lubrication models developed during the last five decades and their integration in their historical numerical frameworks. First, numerical methods limited to the particle dynamics are described. Subsequently, some lubrication correction strategies for coupled fluid-particle methods are discussed.

Particle-based methods.

Accurate simulations of large particle-laden flows have been the ultimate objective for most of the numerical methods developed over the years to simulate suspensions. The applications range from improving our understanding of the physics involved in dense particle flows, to provide reliable solutions for industrial purposes. Due to limitations of computational resources (mainly fast memory access), a class of numerical methods has been developed to simulate particle flows without explicitly solving the fluid phase flow.

Granular media simulations.

The most common approach for granular media simulations is the Discrete Element Method (DEM), and was first introduced by Cundall [START_REF] Cundall | A discrete numerical model for granular assemblies[END_REF] as an efficient approach to simulate rock mechanics. The particle locations and rotations are obtained by solving the Newton-Euler equations of conservation. Only the forces and torques created by collisions and gravity are considered. Hydrodynamic effects are neglected as the particles are assumed to be large and dense enough to not be influenced by the fluid flow (usually gas). This simplification enables to limit the simulation to the particle dynamics. Hence, the problem has 6N degrees of freedom, with N the number of particles, and can be efficiently solved with a good scalability.

From the initial work of Cundall, an extensive family of numerical methods has been developed to simulate granular media [START_REF] Radjai | Discrete-element Modeling of Granular Materials[END_REF][START_REF] Hovi | Physics of dry granular media[END_REF][START_REF] Kishino | Powder and Grains[END_REF][START_REF] Hinrichsen | The Physics of Granular Media[END_REF] . Among them, the soft-sphere collision model (see Sect. 2.3.2) efficiently mimics the particle deformations under stress. Recent works of Haustein et al. [START_REF] Gladkyy | Discrete element modeling of deformable particles in yade[END_REF] include the resolution of the particle deformations into the DEM framework.

When particle chemical properties and interaction laws are added to a DEM, the numerical method is then often referred to as Molecular Dynamics (MD) [START_REF] Luding | Introduction to discrete element methods: Basic of contact force models and how to perform the micro-macro transition to continuum theory[END_REF] . MD was historically introduced to understand the properties of assemblies of molecules in terms of their structure and the microscopic interactions between them. MD simulations act as a bridge between microscopic length and time scales and the macroscopic world of the laboratory that chemists usually observe. Because of the particle length-scale ( 1nm) [START_REF] Rasmussen | Effective particle size from molecular dynamics simulations in fluids[END_REF] usually considered, lubrication forces are negligible compared to molecular interaction forces such as potential forces [START_REF] Hessa | Thermomechanical properties of the wcalennard-jones model system in its fluid and solid states[END_REF] or Van der Waals bonds.

Stokesian Dynamics.

Suspensions are an obvious limitation of DEM, as the hydrodynamic effects play a role essential in the mixture (fluid-particles) dynamics. In the late eighties, Brady et al. [START_REF] Brady | Stokesian dynamics[END_REF] ) introduced a new molecular-dynamics-like approach to solve suspensions without explicitly simulating the fluid flow around the particles. The Stokesian Dynamics (SD) has been designed to solve suspensions in the Stokes regime without neglecting lubrication effects or many-body interactions (influences of interacting particles on particles further away). This method is known for capturing accurately both the near-and far-fields physics.

For N particles suspended in an incompressible Newtonian fluid, the particle dynamics are obtained by the coupled N -body equations of motion

m • dU p dt = F h + F p , (1.1)
where m is the generalized mass/moment-of-inertia matrix of dimension 6N × 6N , U p is the particle linear/rotational velocity vector of dimension 6N , and F h and F p are the hydrodynamic and external force-torque vectors of dimension 6N . As the fluid flows in the Stokes regime, the fluid equations of motion are linear and the hydrodynamic forces and torques acting on the particles are given by

F h = -R FU • (U p -u ∞ ) + R FE : E. (1.2)
The velocity of the bulk linear flow u ∞ acting on the particle is evaluated at the particle center. The externally imposed rate of strain tensor is denoted E, and R FU (x) and R FE (x)

are the resistance matrices which depend on the configuration x -position and orientation -of the particles. In other words, R FU • (U p -u ∞ ) and R FE : E are the hydrodynamic forces/torques on the particles owing to their motions relative to the fluid and owing to the imposed shear flow, respectively. The resistance matrices R SU (x) and R SE (x) are also introduced similarly to R FU (x) and R FE (x), and relate the particle stresslet S h (the symmetric first moment of the force density on a particle) to the velocity and the rate of strain. The combination of the resistance matrices, called the grand resistance matrix, is denoted

R = R FU R FE R SU R SE , (1.3)
and its inverse M = R -1 is called mobility matrix.

Conventional SD exploits the fact that hydrodynamic interactions among particles can be decomposed into long-range mobility interactions and short-range lubrication interactions [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] . The long-range interactions are computed by expanding the force density on the surface of each particle in a series of moments. To minimize the computational 1.2. MULTIPHASE SIMULATIONS.

CHAPTER 1. STATE OF THE ART cost only the two first moments are considered. This level of truncation has been shown to give very accurate results for many hydrodynamic problems [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF][START_REF] Cardinaels | Lubrication analysis of interacting rigid cylindrical particles in confined shear flow[END_REF][START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF][START_REF] Maury | A many-body lubrication model[END_REF] .

The combination of this truncated multipole expansion and the Faxen's laws [START_REF] Faxén | Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist[END_REF][START_REF] Faxén | Der Widerstand gegen die Bewegung einer starren Kugel in einer zähen Flüssigkeit, die zwischen zwei parallelen ebenen Wänden eingeschlossen ist[END_REF] is used to form the far-field grand mobility matrix M ∞ . By construction the far-field grand resistance matrix R ∞ = (M ∞ ) -1 includes an approximation of the dominant many-body interactions but lacks lubrication. Therefore, near-field interactions are introduced into the resistance tensor R ∞ in a pairwise additive fashion. The exact two-body resistance interactions R 2B , are known from the lubrication theory [START_REF] Jeffrey | Calculation of the resistance and mobility functions for two unequal rigid spheres in low-Reynolds-number flow[END_REF] . However, the far-field twoparticle resistance interactions included in R 2B have already been integrated into the grand resistance matrix R ∞ , Therefore, the two-body interactions already included in R ∞ , denoted R ∞ 2B , are subtracted [START_REF] Brady | Dynamic simulation of hydrodynamically interacting particles[END_REF] as well, and the approximation of the grand resistance matrix becomes

R = R ∞ + R 2B -R ∞ 2B . (1.4)
The resulting grand resistance matrix depends only from the position and orientation of the particles. Once R known, the particle velocities can be obtained from the force via Eq. (1.1) and (1.2) or vice versa.

Direct resolution of Eq. (1.1) involves a costly O (N 3 ) calculation cost, due to the inversion of the grand mobility matrix M ∞ , and requires large memory space to store the resistance matrix. However, recent techniques from Sierou et. al. [START_REF] Sierou | Accelerated Stokesian dynamics simulations[END_REF] significantly reduce the cost of the SD to O (N log(N )). The largest drawback of Stokesian dynamics is that the expressions used for the hydrodynamic interactions imply an unbounded medium.

SD can be extended to non-spherical particles by approximating the particles by a combination of several spheres or by computing the resistance matrix of the specific particle geometries [START_REF] Claeys | Suspensions of prolate spheroids in Stokes flow. part 1. dynamics of a finite number of particles in an unbounded fluid[END_REF] .

Multiphase simulations.

In the early days of the XXI st century, improvements in the computational science enable the development of numerical methods solving the particle dynamics as well as the fluid dynamics. This extends the simulation reach to weakly-inertial flows of particles and moves the problem from accurately modeling the fluid flow to modeling the fluid-particle coupling. Numerical methods for modeling two-way solid-fluid coupling can be loosely classified into two categories: partitioned [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] and monolithic [START_REF] Qiu | On thin gaps between rigid bodies two-way coupled to incompressible flow[END_REF] . Partitioned methods typically evolve the fluid and the solids separately, using the results of one as boundary conditions for the other in an alternating one-way coupled fashion. These methods are the most popular as they reuse existing codes targeted to either solids-only or fluids-only problems. However, the partitioned methods face stability issues due to, for instance, the different time scales involved in the flow or the added-mass instability. Monolithic methods aim to more fully two-way couple the fluid and solids alleviating a number of the aforementioned issues [START_REF] Qiu | On thin gaps between rigid bodies two-way coupled to incompressible flow[END_REF] . Most of the lubrication correction techniques described below have been developed for partitioned coupling.

1.2.1 Direct Numerical Simulation.

Direct numerical simulation (DNS) is by construction the most accurate approach to solve particle-laden flows. The equations of motion of the mixture are directly solved numerically without any other approximations. Since no models (lubrication or collision models) are introduced, only errors from numerical methods are included in the solutions.

The accuracy comes with a high computational cost. In practice, direct numerical simulations are limited to some diluted suspensions in a Stokes regime where the additional cost of the numerical resolution of particle collisions can be avoided. Because of this constraint, direct numerical simulations are never strictly performed but are considered with a collision model [START_REF] Wan | Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method[END_REF][START_REF] Delmotte | Hydrodynamics of suspensions of passive and active rigid particles: A rigid multiblob approach[END_REF] . For the sake of clarity, DNS using a collision model will be referred to as quasi-direct numerical simulation (Q-DNS).

Since no lubrication corrections are introduced, the fluid domain needs to be mapped with grid elements small enough to capture interstitial flows. In most cases, unresolved lubrication effects are arbitrarily neglected [START_REF] Wan | Direct numerical simulation of particulate flow via multigrid FEM techniques and the fictitious boundary method[END_REF][START_REF] Delmotte | Hydrodynamics of suspensions of passive and active rigid particles: A rigid multiblob approach[END_REF] . However, lubrication forces cannot be neglected at low particle Stokes number (St d = ρ p U d/9µ 10 3 , where d is the particle diameter, U is the particle velocity, ρ p is the particle density, and µ is the fluid dynamic viscosity.) because of the strong influence of the viscous effects on the particle dynamics. Therefore, several techniques have been suggested to address the resolution of the lubrication effect without introducing a lubrication model.

Mesh refinement.

Face to the massive cost of simulations using uniform size elements small enough to accurately resolve lubrication effect, mesh refinement techniques appear as an obvious approach for cost-efficient Q-DNS. Hence, the mesh is only locally refined in interstitial regions while larger grid elements are considered away from interacting particles. Ideally, space and time meshes need to be both refined to take into account the smaller time scale of the lubrication forces compared to the far-field hydrodynamics [START_REF] Hu | Direct simulation of flows of solid-liquid mixtures[END_REF] .

Mesh refinement is often used with boundary-fitted meshes [START_REF] Joseph | Direct simulation of fluid particle motions[END_REF][START_REF] Johnson | 3D simulation of fluid-particle interactions with the number of particles reaching 100[END_REF][START_REF] Balachandar | Contact problems for particles in a shear flow[END_REF] . This approach enables a more accurate computation of the boundary layers compared to most ficti-Figure 1.1: Example of body-fitted mesh with mesh refinement. Illustration from Johnson et al. [START_REF] Johnson | 3D simulation of fluid-particle interactions with the number of particles reaching 100[END_REF] tious domain approaches [START_REF] Qiu | On thin gaps between rigid bodies two-way coupled to incompressible flow[END_REF] . However, as the particles move, the domain needs to be constantly adapted leading to possibly complex and expansive remeshing issues. The computational cost of the re-meshing limits three-dimensional simulations to about a hundred particles [START_REF] Johnson | 3D simulation of fluid-particle interactions with the number of particles reaching 100[END_REF] .

Nevertheless, Q-DNS with mesh refinement is, as far as we know, the best approach to simulate mixture made of arbitrary particle shapes without neglecting lubrication forces.

Fluid flow decomposition.

Recent works from Lefebvre et al. [START_REF] Merlet | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] enable DNS of Stokes suspensions. The method lies on the idea that Stokes solutions can be decomposed into different contributions, like the SD. Instead of decomposing the total hydrodynamic forces acting on interacting particles into a near-and a far-field (see SD), the linearity of the Stokes equations is used to decompose the velocity field u into a singular flow u sing , which contains lubrication phenomena, and a remainder u reg , which is regular. The regular remaining part of the flow is approximated using a fluid solver. The singular part is decomposed further without approximation over the set of pairs of close particles so that the resolution of the singular field is reduced to the resolution of the Stokes flow around two isolated particles.

Hence, the velocity field of a suspension of N spherical particles is decomposed as

u = u reg + i,j u sing i,j
where u sing i,j is the velocity field around two isolated particles in the same configuration that the two particles P i and P j of the considered suspension. The velocity field u sing i,j is only configuration-dependent (positions, orientations, and relative dimensions of the particles). Therefore, it can be tabulated "off-line" using a fluid solver and refined meshes.

The accuracy of the computation of the lubrication forces depends then directly on the quality of the tabulation: the accuracy of the DNS used for the tabulation and the range of configurations covered by the tabulation. Even if the tabulation needs to be only done once, its computation cost limits, in practice, the method to mono-disperse suspensions of spheres.

Compared to SD, the velocity field decomposition preserves all many-body interactions and not only the dominant effects [START_REF] Merlet | An accurate method to include lubrication forces in numerical simulations of dense Stokesian suspensions[END_REF] . The linearity of the Stokes equations is essential in this method. Therefore, the approach cannot be generalized to Navier-Stokes flows.

Contrarily to methods described in the next section (Sec. 1.2.2), this field decomposition is not considered as a lubrication model because it does not make any simplifications of the lubrication forces acting on the particles. Only lubrication effects resolved by the DNS used for the tabulation are considered, and so errors on the lubrication forces are numerical scheme errors and not intrinsic to a model.

Lubrication Models.

Suspensions involved in industrial applications are rarely limited to less than a few hundred particles. To scale up simulations at a reasonable cost, lubrication effects still need to be modeled. In the past two decades, many lubrication models have been developed to fit the requirement of specific applications or solver properties.

In cases of Stokes flows, the lubrication theory states that collisions between particles occur in an infinite time. Therefore, Verdon et al. have taken advantage of this property by correcting the particle velocities to avoid overlapping (i.e. particle collisions) [START_REF] Verdon | Contact problems for particles in a shear flow[END_REF][START_REF] Verdon | Modified Lees-Edwards boundary conditions and viscous contact for numerical simulations of particles in a shear flow[END_REF] at each time step. The resulting lubrication model is easily scalable and does not require an explicit characterization of the lubrication forces.

When particle inertia is not neglectable, particle collisions might occur. Hence, the expression of the lubrication forces and torques are approximated from the lubrication theory. The most basic approach consists in adding the dominant orders of the lubrication forces to the particle dynamic equations. The lubrication corrections are then activated on interacting particles, in a pairwise fashion, when they are closer to each other than a critical distance ∆L lub . The critical distance is often set arbitrarily [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF][START_REF] Kempe | Collision modelling for the interface-resolved simulation of spherical particles in viscous fluids[END_REF][START_REF] Lin | Numerical investigation of lubrication force on a spherical particle moving to a plane wall at finite Reynolds numbers[END_REF] such that 0 < ∆L lub < d, with d the particle diameter. As far as we know, there is no characterization of the minimal gap length between two interacting particles below which lubrication effects appear. By setting arbitrarily ∆L lub , lubrication forces are likely to be added to the fully resolved hydrodynamics. To avoid overestimating the lubrication, ∆L lub can be 1.3. SUMMARY. CHAPTER 1. STATE OF THE ART tabulated such that the total hydrodynamic forces tend to the theoretical lubrication forces as the gap vanishes [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF][START_REF] Breugem | A combined soft-sphere collision/immersed boundary method for resolved simulations of particulate flows[END_REF][START_REF] Nguyen | A fictitious-domain simulation of solid-liquid flow with subgrid lubrication force correction; a sphere falling onto a plane surface[END_REF][START_REF] Gazanion | Numerical modelling of finite-size particle collisions in a viscous fluid[END_REF] . Such a lubrication model is described in depth in Sec. 3.2 for mono-disperse suspension of spheres. The number of parameters needed for the tabulation of ∆L lub depends on the particle shapes and orientations. For monodisperse flows of spherical particles, the grid resolution is the only tabulated parameter of ∆L lub . For poly-disperse flows, the relative aspect ratio between the pair of particles must be considered as well, and the particle relative orientations become additional parameters for non-spherical particles.

To limit approximation errors from the lubrication model in the solution, a mesh refinement approach can be combined with a lubrication model [START_REF] Phan-Thien | Lubrication approximation in completed double layer boundary element method[END_REF] . By refining the mesh in interstitial regions, the use of the lubrication model is reduced to smaller ∆L lub than on a uniform mesh. Furthermore, this hybrid approach decreases the need for heavy refinement making many-particle simulations more cost-effective than a Q-DNS with mesh refinement.

By correcting the hydrodynamic forces, a pair of interacting particles at the time, many-body interactions are not corrected. For dense suspensions, the dominant manybody interactions can be corrected using an approach inspired by the SD. Gallier et al.

have proposed a model using tabulated resistance matrices to compensate the unresolved lubrication forces [START_REF] Gallier | A fictitious domain approach for the simulation of dense suspensions[END_REF][START_REF] Gallier | Rheology of sheared suspensions of rough frictional particles[END_REF] . By construction, this model should be strictly limited to noninertial flows. However, Gallier et al. have shown the solutions remain fairly accurate for weakly-inertial flows.

Most of the lubrication correction methods mentioned above have been developed and tested for suspensions of spheres. With the improvements in computational sciences in the last decade, interest for suspensions of non-spherical particles is growing, and a few techniques for spherical particles have been extended to spheroidal particles [START_REF] Breugem | Numerical study of the sedimentation of spheroidal particles[END_REF][START_REF] Harting | Accurate lubrication corrections for spherical and non-spherical particles in discretized fluid simulations[END_REF] .

Summary.

The resolution of the short-range hydrodynamics is critical to accurately simulate suspensions. The direct resolution of the flow equations first appears as the best method to fully capture the flow details. Unfortunately, this method requires large computational resources, even with advanced remeshing techniques.

The alternative approach is to introduce a lubrication correction model based on the lubrication theory solutions. Approximations and assumptions inherent to the lubrication models limit the global numerical methods to specific flows. The main limitations lie in the flow regime (Stokes regime) and/or particle shapes (spherical particles). Extensive researches have been made on the suspension of spherical particles and have provided accurate and robust methods to model lubrication effects. Supporting by recent improvements in computational sciences, the research community is now looking to extend lubrication models to more complex particle flows (turbulent flows, complex particle shapes, ...). This thesis study aims at developing a lubrication model compatible with incompressible Navier-Stokes flows of convex particles.

Chapter 2 Numerical Resolution of Particle-Laden Flows

The motion of particle-laden flows is governed by the coupling of the equations of motion for the continuous and dispersed phase. The solid phase (particles) is described by the Newton-Euler equations, while the fluid phase is described by the Navier-Stokes equations. In this chapter, the partitioned numerical method used to solve particle-laden flows is described in depth.

After introducing the notations and the physical problem, the resolution of the incompressible Navier-Stokes equations by a volume penalization method is detailed. The dynamics of the rigid particles are solved with a discrete element method coupled with the fluid dynamics.

2.1 General description of the problem.

       ∂u i ∂t + ∂ (u i u j ) ∂x j = - 1 ρ ∂p ∂x i + µ ρ ∂ ∂x j ∂u i ∂x j + χλ (u τ,i -u i ) , ∂u i ∂x i = 0, (2.1)
where (i; j) = {1, 2, 3} 2 , u i are the velocity components, p is the pressure. The density and dynamic viscosity of the Newtonian fluid are denoted ρ and µ. The system of Eq.

(2.1) is solved on a uniform Cartesian mesh of the three dimensional domain

D = D f ∪ D s ,
where N spherical particles P i forming the solid domain D s = N i=1 P i with P i ∩ P j = ∅ for i = j and D f = D \ D s is the fluid domain. The interface between the solid and fluid phases is denoted Γ s = N i=1 Γ i . Each particle P i is assumed to be homogeneous with a density ρ p,i .

A no-slip boundary condition is implicitly imposed at the interface Γ s , by the penalty term χλ (u τ,i -u i ). Indeed, as reminded by Angot et al. [START_REF] Angot | A penalization method to take into account obstacles in incompressible flows[END_REF] , solving the penalized Eqs.

(2.1) is equivalent to solve the incompressible Navier-Stokes equations in the fluid domain and to enforce a no-slip boundary condition at the boundary Γ s when λ → +∞. The computation of the penalty term is detailed later (see Sec. 2.2.4).

Representation of the particles.

The dynamics of each rigid particle is described by the Newton-Euler equations:

m i dU i dt = F i , (2.2)
dJ i Ω i dt = T i , (2.3)
for a given particle P i of mass m i , inertia matrix J i , linear velocity U i of the mass center and rotational velocity Ω i . The computation of the total force F i and torque T i acting on the particle P i are detailed in Sec. 2.3. The location of the mass center and the orientation of the particle are given by X i and Θ i .

The surface of each particle is meshed with N p elements (Fig. 2.2). These meshes are used to compute the level-set function and the hydrodynamic forces acting on the particle (Sec. 2.3.1).

Figure 2.2: Representation of the particle surface mesh and the particle dimensions in the particle coordinate system e x , e y , e z .

Particle classification.

Particle surface is defined by all (x, y, z) ∈ R 3 such that:

x 2 a 2 + y 2 b 2 + z 2 c 2 = 1, (2.4)
where a, b, and c are the three dimensions of the ellipsoid (Fig. 2.2). The particle is called a spheroid when b = c. For spheroidal particle, the particle aspect ratio is defined as AR = a/b. The spheroid is an oblate if AR < 1, a prolate if AR > 1 and a sphere otherwise. The equivalent sphere, of diameter D eq = 2 3 √ abc, is the sphere with the same volume as the particle.

Computation of the global level-set.

The global level-set function LS is defined on all cells of the background grid (where the penalized Navier-Stokes equations are solved) as the minimal signed distance of the given cell center to the particle surfaces. The level-set is positive on fluid cells and negative on solid cells. For spherical particles, the level-set at the cell center X is given by

LS(X) = min i∈ 1,N ( X -X i -a i ) . (2.5)
The ellipsoid surfaces are not easily characterized in the global coordinate system.

Therefore, the level-set is computed via the particle surface mesh or via the resolution of a minimization problem.

• Orthogonal projection method:

Via the particle surface mesh, the level-set at a given cell center x is computed as follows:

1. The closest particle surface element is found among all the particles of the system.

2. The position x is orthogonally projected on the closest particle surface element.

This projection is denoted X Γ .

3. The level-set LS(x) is then given by x -X Γ and the sign of LS(x) is positive if X Γ is closer to the particle mass center than x. Otherwise, the sign of LS(x) is negative.

• Minimization approach:

The other approach is based on the minimization of the distance between a given grid cell and the analytic surface of the ellipsoid. The grid cell location needs to be in the particle coordinate system e x , e y , e z . Its position is denoted X p = (x p , y p , z p ). The closest point X = (x, y, z) on the particle surface to X p can be find by minimizing the objective function

F (X) = X -X p 2 + λ x 2 a 2 + y 2 b 2 + z 2 c 2 -1 2 , (2.6)
where λ > 0 is the Lagrangian multiplier chosen arbitrary large. When the contact point X is known, the value of level-set is given by X -X p and its sign is the sign of (X p -X)•X. The global level-set LS(x) is then defined as the minimal level-set computed among the closest particles to the given grid point.

Hence, the particles are located on the grid by the characteristic function χ defined on each cell center location x such that χ(x) = 1 if LS(x) ≤ 0, and χ(x) = 0 otherwise.

Computation of the minimal distance to an obstacle.

When all the particles are spherical, the computation of the minimal distance between a particle and a wall or another particle is elementary. However, the problem becomes more complex if ellipsoidal particles are considered. The localization of the contact points for a pair of particles can be computed by solving a minimization problem or by geometrical constructions.

• Minimization problem:

The distance between two particles is obtained by locating the contact points at the surface of both particles. The contact point locations X c,i = (x c,i , y c,i , z c,i ) and X c,j = (x c,j , y c,j , z c,j ) in the particle coordinate systems are computed by minimizing the following objective function:

F (X c,i , X c,j ) = R i X c,i + X i -(R j X c,j + X j ) 2 + λ i x 2 c,i a 2 i + y 2 c,i b 2 i + z 2 c,i c 2 i -1 2 + λ j x 2 c,j a 2 j + y 2 c,j b 2 j + z 2 c,j c 2 j -1 2 , (2.7)
where R i and R j are the rotation matrices from the particle P i and P j coordinate systems to the global coordinate system. The Lagrange parameters λ i and λ j are positive and are arbitrarily chosen large.

The objective function Eq. (2.7) can be minimized using a gradient method. Preliminary tests using the gradient descent method and Newton's method have shown a high sensitivity of the solution to the initial guess and a slow convergence speed. The convergence speed can be improved with adaptive methods. Nevertheless, the geometrical construction method described below has proven to be the more robust and efficient method to compute the minimal distance.

• Geometrical constructions:

The efficient alternative method proposed by Lin et al. [START_REF] Lin | On the distance between two ellipsoids[END_REF] is based on the iterative construction of virtual spheres. The method can be summarized as follows (see Fig. 2.3):

1. The search algorithm starts from two arbitrary points (X c,i , X c,j ) k on the surface of the two particles. These two points are assumed as the nearest points to the other particle surface.

2. Two spheres are constructed completely inside the ellipsoids such that the sphere and ellipsoid surfaces are tangent at the current nearest points (X c,i , X c,j ) k .

3. A new guess of the contact points (X c,i , X c,j ) k+1 is then found by the intersection of the line generated by the centers of the two spheres and the ellipsoid surfaces.

4.

If not converged, go back to step 2. Convergence is obtained when the line generated by the centers of the two spheres is co-linear with the ellipsoid surface gradients at the contact points (X c,i , X c,j ) k+1 .

The procedure converges faster as the radius of the virtual spheres increases, and the initial guesses of (X c,i , X c,j ) k are close to the contact points.

Figure 2.3: Two-dimensional sketch of the iterative method of Lin et al. [START_REF] Lin | On the distance between two ellipsoids[END_REF] used to compute the distance between two ellipsoids.

The computation of the minimal distance for a pair particle-wall is more straight forward than for a pair of particles. Indeed, the outgoing normal unit vector of Γ i at the contact point X c,i is co-linear and has an opposite sign to the wall outgoing unit surface vector. Hence, the contact point

X c,i = (x c,i , y c,i , z c,i ) is given by                x c,i = -x w a 2 i 2 , y c,i = -y w b 2 i 2 , z c,i = -z w c 2 i 2 , (2.8)
with e w = (x w , y w , z w ) the wall outgoing unit surface vector in the particle coordinate system.

Virtual spheres.

For non-spherical particles, the collision model (Sec. 2.3.2) and lubrication correction models (Sec. 3.2 and 3.3) are evaluated indirectly via virtual spheres. These virtual spheres are defined such that the sphere surfaces fit as much as possible the ellipsoid surfaces at the contact points. Hence, the radius of the virtual sphere R vs,i is given by the Gaussian radius of curvature [START_REF] Breugem | Numerical study of the sedimentation of spheroidal particles[END_REF] of the ellipsoid

P i at the contact point X c,i = (x c,i , y c,i , z c,i ): R 2 vs,i = a 4 i b 4 i + b 4 i (c 2 i -a 2 i ) x 2 c,i + a 4 i (c 2 i -b 2 i ) y 2 c,i 2
a 6 i b 6 i c 2 i . (2.9)
Subsequently, the center of the virtual sphere X vs,i , in the particle coordinate system, is defined as

X vs,i = X c,i -R vs,i N c,i , (2.10)
where N c,i is the outgoing normal unit vector of Γ i , in the particle coordinate system, at the contact point X c,i . The coordinates of the virtual sphere center are then projected in the global coordinate system using the particle rotation matrix.

Discretization of the governing equations.

The Navier-Stokes equations (2.1) are discretized on a uniform Cartesian mesh using a cell-centered collocated arrangement of the primitive variables (p, u i ). Face-centered velocities v i are also introduced in addition to the cell-centered velocities u i , to eliminate odd-even decoupling which can lead to large pressure variations in space [START_REF] Bozkurttas | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] .

The equations are integrated in time using a classical projection scheme introduced by Chorin [START_REF] Chorin | Numerical solution of the Navier-Stokes equations[END_REF] and Temam [START_REF] Temam | Sur l'approximation de la solution des équations de Navier-Stokes par la méthode des pas fractionnaires II[END_REF] . The system (2.1) is then solved in the following four steps.

Prediction.

In this first step, the momentum equation is solved to obtain an intermediate virtual velocity u * i , starting from a guess for the pressure field q. We choose an incremental scheme therefore q = p n . The convective and viscous terms are discretized in time by a second-order Adams-Bashforth scheme and an implicit Crank-Nicolson scheme, respectively. Hence, the following momentum Eq. (2.11) is solved at the cell-nodes:

u * i -u n i ∆t n + 1 2 ∆t n + 2∆t n-1 ∆t n C n i -C n-1 i = - 1 ρ ∂q ∂x i + µ 2ρ (D * i + D n i ) , (2.11)
where The virtual face-centered velocities v * i are calculated as follows [START_REF] Bozkurttas | A versatile sharp interface immersed boundary method for incompressible flows with complex boundaries[END_REF] :

C i = ∂v j u i ∂x j and D i = ∂ ∂x j ∂u i ∂x
                       ûi = u * i + ∆t ρ ∂q ∂x i cc , v1 = γ w ûP,1 + (1 -γ w ) ûW,1 , v2 = γ s ûP,2 + (1 -γ s ) ûS,2 , v3 = γ b ûP,3 + (1 -γ b ) ûB,3 , v * i = vi - ∆t ρ ∂q ∂x i fc , (2.12)
where ûP , ûW , ûS , and ûB are the velocities computed at the nodes P, W, S, and B, respectively (see Fig.

Projection.

At the end of the prediction step, the velocities u * i are a priori not divergence free. The projection step aims at finding the pressure field p n+1 such that ∇ • u ∼ 0. In particular, the equation

u i -u * i ∆t = 1 ρ ∂p ∂x i (2.13)
is solved with the constraint that ∇ • u = 0. Hence, the Poisson equation for the pressure

correction p is 1 ρ ∂ ∂x i ∂p ∂x i = 1 ∆t ∂v * i ∂x i , (2.14)
with Neumann boundary conditions on the pressure on all external boundaries (of D).

Spatial derivatives ∂ ∂x are approximated by a second-order central finite-difference scheme. The Poisson problem Eq. (2.14) is solved using the Generalized Minimal Residual method (GMRES) from Saad et al. [START_REF] Saad | GMRES: A generalized minimal residual algorithm for solving nonsymmetric linear systems[END_REF] . Preconditioners1 are used to reduce the number of iterations needed to ensure a divergence free velocity field, especially for low Reynolds flows.

Correction.

From the solution p of the Poisson Eq. (2.14), the pressure and velocity fields (p n+1 , u) are updated as

p n+1 = q + p , u i = u * i -∆t 1 ρ ∂p ∂x i cc , v i = v * i -∆t 1 ρ ∂p ∂x i fc , (2.15)
so that u is divergence free. Gradients are discretized by second-order central finite difference scheme. Cell-centered and face-centered velocities are updated separately to improve the accuracy and stability of the prediction step.

The initial pressure is set arbitrarily to p 0 = 0.

Penalization.

This last step is performed after the resolution of the particle dynamics described in Sec. 2.3 (see Fig. 2.8). The velocity field u is penalized using the particle velocities to compute u n+1 :

u n+1 i = u i + χλ∆tu τ,i 1 + χλ∆t , (2.16)
where λ is the penalty factor (arbitrary chosen here such that λ = 10 8 ). The characteristic function χ is defined for each cell center location x such that χ(x) = 1, if x ∈ D s and χ(x) = 0 anywhere else.

On each solid element (χ(x) = 1), the velocity is penalized using the solid velocity

u τ = U i + (X i -x) × Ω i , with x ∈ P i .
This method is first-order accurate in space.

To improve the consistency of the solution, corrections are made on solid elements at the interface Γ s . The interface characteristic function χ Γ is defined for each cell center location x such that χ Γ (x) = 1, if x ∈ D s and the given cell has at least one neighbor in D f . Everywhere else, χ Γ (x) = 0. The velocity corrections are made using the second-order ghost-cell approach known as Image Point Correction method (IPC) [START_REF] Hovnanian | An accurate cartesian method for incompressible flows with moving boundaries[END_REF] . The main idea of the method is to find the right correction u τ,i for all solid cells at the fluid-solid interface (χ Γ = 1) to impose the desired velocity at the interface Γ s .

Figure 2.5: Two-dimensional sketch of the IPC method. Velocity at the solid ghost point G is corrected using its symmetric F. The velocity at F is obtained by interpolation of the velocity of its four neighboring fluid points (eight neighbors in three-dimensional cases). The orthogonal projection of G at the particle surface is denoted B.

The IPC method, illustrated in Fig. 2.5, corrects the velocity at the solid ghost point G by

u τ = u B + φ ∂u ∂n φ=0 , (2.17)
where u B is the solid velocity at the surface point B, n is the outward normal unit vector of the interface. The level-set function φ is the signed distance from the boundary of the bodies with negative sign within the particles and positive one elsewhere. In particular, φ = 0 at the interface Γ s , and |φ| is the minimal distance to Γ s . Where there is no particle overlapping, φ is equivalent to the global level-set function LS (see Sec. 2.1.2).

The velocity u B can be easily calculated from the velocity U and rotational velocity Ω of the particle mass center. The gradient (∂u/∂n)| φ=0 is determined from u B and the velocity u F at F the symmetric point of G. The distance between B and F, as well as the outward normal vector n can be easily computed with the level-set function φ.

The velocity at the point F is computed by interpolating the velocities of the eight closest fluid neighbors of F. If one of the eight neighbors is G, the surface point B is considered instead of G. Interpolations are made using a second-order Lagrange interpolation scheme.

Numerical stability: time step adaptation.

To ensure the stability of the whole numerical method, the time step ∆t is adapted such that the Courant -Friedrichs -Lewy condition is satisfied: ∆t = β τ ∆x/V max where β τ ≤ 1 is chosen arbitrarily, ∆x is the characteristic length of the grid cells and V max is the maximum of the velocity absolute value computed on the grid cells.

When

Dynamics of the particles.

The particle dynamics is solved by a discrete element method (DEM) which is primarily devoted to multi-contact interactions for a large suspension [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] . The dynamics of each rigid particle is obtained by the Newton-Euler equations of conservation:

m i dU i dt = F hyd i + F coll i + F ext i , (2.18
)

dJ i Ω i dt = T hyd i + T coll i , (2.19

Particle hydrodynamics.

To numerically compute the hydrodynamic force and torque acting on the particle P i , the surface Γ i is meshed using N p,i elements. The k-th element of the mesh of Γ i is denoted P k i and s k is its surface. The set L L L(P k i ) is defined as all the particles P j , j = i such that the distance between the surface of P j and the center of P k i is lower than a lub . The distance a lub is defined as the narrowest gap width between the center of the P k i and a potential nearby obstacle for the solver to fully resolve hydrodynamic interactions (see The total hydrodynamic force and torque on the particle P i are given by

     F hyd i = F solv i + F deg i , T hyd i = T solv i + T deg i .
(2.20)

The force F deg i and torque T deg i are defined as:

                 F deg i = p∈Γ i , L L L(p) =∅ (σ • n) dS, T deg i = p∈Γ i , L L L(p) =∅ a i n × (σ • n) dS.
(2.21)

These two components of the hydrodynamics are underestimated by the numerical simulation due to the insufficient number of grid elements in the gap between the particle

                 F solv i = p∈Γ i , L L L(p)=∅ (σ • n) dS, T solv i = p∈Γ i , L L L(p)=∅ a i n × (σ • n) dS.
(2.22)

The resolved hydrodynamics are computed by numerical integration of the fluid stress σ acting on all elements P k i far enough to nearby obstacles:

             F solv i ≈ k∈ 1,N p,i , L L L(P k i )=∅ (σ • n) s k , T solv i ≈ k∈ 1,N p,i , L L L(P k i )=∅ a i n × (σ • n) s k , (2.23)
The fluid stress σ is interpolated of the pressure and velocities (p, u i ), using a secondorder Lagrange scheme at the center of P k i .

2.3.2

The soft-sphere collision model.

In a Stokes flow, the contact between obstacles is theoretically impossible due to the lubrication singularity. However, a collision model needs to be considered for physical and numerical purposes.

From a physical perspective, the lubrication effect alone cannot explain some phenomena such as the rebound of particle onto a wall, which occurs at the Stokes number above roughly [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] St d = ρpdUc 9µ ≈ 10, with ρ p , d, and U c the particle density, diameter, and impact velocity, respectively. The collision model is seen here as a "low cost" model to mimic the effect of the particle deformation on its trajectory. Furthermore, real surfaces are never perfectly smooth, meaning that contact can occur at the particle roughness scale or at the lubrication film breakdown scale (when the separation distance is close to the fluid particle characteristic size).

From a numerical point of view, a cutoff distance for the lubrication forces is usually introduced to ensure the model stability. This cutoff distance prevents the lubrication force Eq. (3.38) to diverge as the particle collides with an obstacle. A control on the time step to avoid contact due to time discretization errors is not an option. Indeed, an adaptive time step would certainly improve the capture of the lubricating effects, but in case of colliding particles at the Stokes number regime St d > 10, the time step would keep decreasing asymptotically to zero. Therefore, a lubrication cutoff distance and a collision model need to be introduced for stable and finite time simulations.

The collision model chosen here is based on the soft-sphere approach used by Costa et al. [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] . The deformation of particles during contact is modeled by the overlap between a particle and an obstacle (particle or wall). From the overlap measurement, normal and tangential contact forces are computed using a mass-spring-dashpot system and a Coulomb-type threshold for the tangential component. For a given particle P i , collision

       F coll i = j =i F i,j + F i,wall , T coll i = j =i T i,j + T i,wall , (2.24)
where F i,j is the collision force of the interacting particles P i and P j , F i,wall is the collision force of P i with a wall. T i,j and T i,wall are the corresponding collision torques. The force and the torque on P i resulting from the particle-particle interactions between P i and P j are defined using a local system of coordinates (e n , e t ) (Fig. 2.7):

F i,j = F n + F t , T i,j = ae n × (F t ) , (2.25)
with

F n = -δ n δ n δ n k n -γ n (U i,j • e n ) e n , F t = min ( -δ t δ t δ t k t -γ t (U i,j • e t ) e t , µ c F n ) e t , (2.26)
where a is the radius of

P i , δ n δ n δ n (respectively, δ t δ t δ t ) is the normal (respectively, tangential)
overlap, k n (respectively, k t ) is the normal (respectively, tangential) stiffness, µ c is the coefficient of sliding friction, and γ n (respectively, γ t ) is the normal (respectively, tangential) damping coefficient of the spring-dashpot model. The relative velocity of the two particles U i,j at the contact point is given by U

i,j = U i + aω i × e n -(U j -aβω j × e n ).
The normal overlap distance δ n δ n δ n is given by

δ n δ n δ n = max (0, a (1 + β) + col (a + βa) -X i -X j ) e n , (2.27)
where e n = X i -X j X i -X j as shown Fig.

e t = -δ t δ t δ t k t -γ t (U i,j -(U i,j • e n ) e n ) -δ t δ t δ t k t -γ t (U i,j -(U i,j • e n ) e n ) .
The parameters of the spring-dashpot model γ n , γ t , k n , and k t are calculated from the coefficient of normal (respectively, tangent) restitution ξ max,n (respectively, ξ max,t ) of dry collision and the contact time τ c , as follows:

                       k n = m * π 2 + ln 2 (ξ max,n ) τ 2 c , γ n = - 2m * ln (ξ max,n ) τ c , k t = m * t π 2 + ln 2 (ξ max,t ) τ 2 c , γ t = - 2m * t ln (ξ max,t ) τ c , (2.30)
with the effective mass

m * = m i m j m i + m j , and
m * t = m * K 2 K 2 + 1 ,
with K 2 = 2/5 the normalized radius of gyration for spherical particles.

The characteristics of the elastic properties of the particles are ξ max,n , ξ max,t , and τ c . As noticed by Izard et al. [START_REF] Bonometti | Modelling the dynamics of a sphere approaching and bouncing on a wall in a viscous fluid[END_REF] , the relation between τ c and k n is unusual, but several studies show that the normal stiffness can be underestimated without modification of the dynamics of a dry system. Such assumption allows to reduce the simulation time since the collision characteristic time will be larger than the particle characteristic deformation time.

To ensure the stability of the model and the conservation of the momentum, Costa et al. [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] advised that the time step of the overall numerical algorithm ∆t has to be chosen as a multiple of the contact time τ c (at least during the collision). This condition guarantees a zero overlap at the end of the collision and allows the fluid to adapt itself to sudden changes in velocity of the colliding particles.

The force F i,wall and the torque T i,wall are assumed to be equivalent to the asymptotic case β → +∞ and m P j → +∞.

For non-spherical particles, the ellipsoids are approximated as spherical particles with the same mass as the whole particles and with a radius corresponding to the local Gaussian curvature at the contact point [START_REF] Breugem | Numerical study of the sedimentation of spheroidal particles[END_REF] . During collisions, the radii of the approximating spheres remain constant simplifying the problem to the collision between two unequal spheres. The centers of the approximating colliding spheres are stored at the time step before the gap width becomes negative and updated during the collision using the particle velocity and the rotation matrix introduced above. These constraints prevent nonconverged computations of the virtual spheres (Sec. 2.1.2) to destabilize the collision model.

Numerical resolution of the particle dynamics.

Among the forces acting on the particle, the short-range hydrodynamics (lubrication)

and collision forces have time scales smaller than the time scale associated with fluid flow. The particle dynamics is therefore computed at a smaller time step δt = ∆t nt , with n t arbitrary chosen large, to accurately integrate (in time) short-range interaction forces. Numerical simulations have shown that changing n t between 10 2 and 10 3 does not affect significantly the results.

Since the motion of the particles occurs at a time step smaller than ∆t, it is not necessary to re-compute the resolved hydrodynamic forces F solv i at each sub-time step δt. Hence, particle dynamics is solved at each sub-time step δt with updated short-range interactions and "frozen" resolved hydrodynamics. of the whole system from its state n. The algorithm starts at the estimation of the time step ∆t which is the time elapse between the system state n and n+1. The particle dynamics is solved in n t sub-time steps δt = ∆t nt . m denotes the current sub-state of the system between the state n and n + 1 while the particle dynamics is computed (NB: after the correction step, the state m is equivalent to the state n). Dashed boxes contain the reference to the main equations computed at the given step.

Semi-Implicit Euler scheme.

The discretized expressions of Eqs. (2.18) and (2.19) are then written as:

U m+1 i = U m i + δt m i F solv i n+1 + δt m i F coll i + F ext i + F lub i m .
(2.31)

Ω m+1 i = J m+1 i -1 (J m i Ω m i +δt T solv i n+1 + T coll i + T lub i m , (2.32)
where m denotes the "sub-state" of the system at the time t = t n + mδt with t n the time at the state n of the system. Since the particle dynamics are solved after the correction step but before the penalization (see Fig. Position and orientation are then integrated implicitly as follows:

X m+1 i = X m i + δt U m+1 i , Θ m+1 i = Θ m i + δt Ω m+1 i .
(2.33)

2.3.3.2 Verlet integration method.
Verlet integration is a numerical method used to integrate Newton's equations of motion (Eqs. (2.18) and (2.19)). The Verlet integrator is a second order accuracy scheme with a cost comparable to the semi-implicit Euler method (which is first order accuracy).

The particle positions and velocities are computed at t m+1 from the state of the system at t m via the following steps:

1. Velocities are computed at half time step δt/2 as follows:

                     U m+1/2 i = U m i + δt 2m i F solv i n+1 + δt 2m i F coll i + F ext i + F lub i m , Ω m+1/2 i = J m+1 i -1 (J m i Ω m i + δt 2 T solv i n+1 + T coll i + T lub i m .
(2.34)

2. Particle location and orientation are updated using the half step velocities U m+1/2 i and Ω m+1/2 i :

X m+1 i = X m i + δt U m+1/2 i , Θ m+1 i = Θ m i + δt Ω m+1/2 i . ( 2
                     U m+1 i = U m+1/2 i + δt 2m i F solv i n+1 + δt 2m i F coll i + F ext i + F lub i m+1/2 , Ω m+1 i = J m+1 i -1 J m i Ω m+1/2 i + δt 2 T solv i n+1 + T coll i + T lub i m+1/2
.

(2.36)

Numerical implementations.

The significant advantage of partitioned methods is the reused of existing codes. In this thesis, incompressible Navier-Stokes are solved using the existing in-house solver NaSCar to solve the fluid flow as described in Sec. 2.2. NaSCar2 is devoted to solve 3D-flows in around moving and deformable bodies. This code is written in C language and use PETSc library3 for the resolution of large linear systems in parallel.

The solver of the particle dynamics (Sec. 2.3) has been developed from existing features of NaSCar. The dynamics of solids coded in NaSCar has been completed by adding lubrication and collision models and has been optimized to simulate suspensions of rigid particles.

Chapter 3

Lubrication Correction Models

The accuracy of the computed hydrodynamic forces acting on a particle directly depends on the accuracy of the computation of the fluid flow surrounding the given particle.

In dense suspensions, interstitial flows are often poorly resolved by the direct numerical resolution of the Navier-Stokes equations as the gap between interaction particles can be smaller than the grid resolution. To balance the unresolved hydrodynamics, lubrication correction models are introduced. This chapter starts with the detailed calculation of the theoretical lubrication forces acting two isolated interacting particles. Two lubrication correction models, based on the lubrication theory, are then described.

Theoretical lubrication forces and torques.

Two smooth particles of radii a and βa (with β > 0) immersed in an infinite fluid domain are considered. One of the particle (P i ) is moving at the constant velocity U i near the second stationary particle (P j ) (Fig. 3.1). The surfaces of P i and P j are denoted respectively Γ i and Γ j . The gap between the two particles is denoted a with 0 < 1.

Figure 3.1: Sketch of two particles in interaction.

THEORY CHAPTER 3. LUBRICATION CORRECTION MODELS

In any plane containing the axis (Oz), the particle surfaces Γ i and Γ j can be characterized, near the contact point, as follows:

         z i (r) = a (1 + ) - √ a 2 -r 2 ∼ a + r 2 2a + r 4 8a 3 + O r 6 , z j (r) = -βa + β 2 a 2 -r 2 ∼ - r 2 2βa - r 4 8a 3 β 3 + O r 6 . (3.1)
To simplify expressions, the stretched coordinates are introduced:

             R = r a √ , Z = z a , H = Z i -Z j = 1 + 1 + β 2β R 2 + 1 + β 3 8β 3 R 4 + O 2 , (3.2)
where Z i and Z j are respectively z i (r) and z j (r) from Eq. (3.1) in the stretched coordinates.

If the Reynolds number in the gap on the neighborhood of the contact point A is small (Re = ρU i a /µ 1), then the flow can be modeled by the Stokes equations

                 µ∆u = ∇p, ∇ • u = 0, u Γ i = U Γ i , u Γ j = U Γ j , lim r→+∞ p = 0, (3.3)
where u is the velocity field, p is the pressure field, ρ is the fluid density, µ the fluid dynamic viscosity. The next two sections address the particular cases of the squeezing motion (U i = -U sq i e z ) and the shearing motion (U i = U sh i e x ) of P i near the stationary particle P j (U Γ j = 0).

Squeezing motion of rigid spheres.

In this section the particle surface velocities are U Γ i = -U sq i e z and U Γ j = 0. To simplify the resolution of the Stokes problem, the first line of the system (3.3) is projected in the cylindrical coordinate system (r, θ, z):

   ∂ r p = µ r ∂ z Φ 2 (ψ) , ∂ z p = - µ r ∂ r Φ 2 (ψ) , (3.4)
with

Φ 2 (•) = r∂ r 1 r ∂ r (•) + ∂ 2 z (•)
and the Stokes stream-function ψ is defined [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF][START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF][START_REF] Stimson | The motion of two spheres in a viscous fluid[END_REF] such that (u r , u θ , u To simplify the notations, the system (3.5) is solved using the stretched coordinates (3.2). The governing equations (3.5) then become

z ) = 1 r ∂ z ψ, 0, -1 r ∂ r ψ . By noticing that ∂ r (∂ z p) -∂ z (∂ r p) = 0, the system (3.3) becomes        Φ 4 (ψ) = 0, ψ Γ i = 1 2 U sq i r 2 , ψ Γ j = 0,
       [∂ 2 Z + Υ] 2 ψ = 0, ψ Γ i = 1 2 U sq i R 2 a 2 , ψ Γ j = 0, (3.6)
where

Υ(•) = ∂ 2 R (•) - 1 R ∂ R (•).
As suggested by Jeffrey et al. [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF] , the solution of Eq. (3.6) can be expressed as follows:

ψ(R, Z) = a 2 ψ 0 (R, Z) + ψ 1 (R, Z) + 2 ψ 2 (R, Z) + O 3 , (3.7)
where ψ 0 , ψ 1 , and ψ 2 are the three first harmonics of ψ. By substituting the expression of ψ in Eq. (3.7) into the first line of Eq. (3.6), [∂ 2 Z + Υ] 2 ψ can be seen as a polynomial function in equals to zero. The identification of the polynomial's coefficients gives the following hierarchy of equations:

       ∂ 4 Z ψ 0 = 0, ∂ 4 Z ψ 1 + 2Υ (∂ 2 Z ψ 0 ) = 0, ∂ 4 Z ψ 2 + 2Υ (∂ 2 Z ψ 1 ) + Υ 2 (ψ 0 ) = 0. (3.8)
To close the system (3.8), the boundary conditions on the harmonics of ψ need to be retrieved from the boundary conditions on ψ Eq. (3.6). As suggested by Jeffrey et al. [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF] ,

the boundary conditions on the harmonics of ψ can be obtained via the Taylor expansion of ψ at the contact points on Γ i and Γ j .

On the surface Γ i , the Taylor expansion of ψ at Z = H i is given by

ψ (R, Z) = ψ (R, H i ) + ∂ Z ψ (R, H i ) (Z -H i ) + ∂ 2 Z ψ (R, H i ) (Z -H i ) 2 2 + O Z 3 , = a 2 (ψ 0 + ψ 1 + 2 ψ 2 + O ( 3 )) + a 2 (∂ Z ψ 0 + ∂ Z ψ 1 + 2 ∂ Z ψ 2 + O ( 3 )) 1 8 R 4 + 1 16 2 R 6 + O 3 + a 2 2 ∂ 2 Z ψ 0 + O ( ) 1 8 R 4 + O 2 2 + O 5 , (3.9)
with

H i = 1 + 1 2 R 2 the dominant order of Z i .
From the boundary conditions (3.6),

ψ (R, Z i ) ∼ 1 2 U sq i R 2 a 2 .
Hence, by identifying the coefficients of the Taylor expansion of ψ at Z = H i (seen as a polynomial function in ) evaluated at Z = Z i , the following hierarchy of equations is found:

           ψ 0 (R, Z i ) = 1 2 R 2 U sq i , ψ 1 (R, Z i ) = - 1 8 R 4 ∂ Z ψ 0 (R, Z i ) , ψ 2 (R, Z i ) = - 1 8 R 4 ∂ Z ψ 1 (R, Z i ) - 1 16 R 6 ∂ Z ψ 0 (R, Z i ) - 1 128 R 8 ∂ 2 Z ψ 0 (R, Z i ) . (3.10)
The boundary conditions on the first derivative of the harmonics of ψ can be found using the same approach. On the surface Γ i , the Taylor expansion of

∂ Z ψ at Z = H i is
given by

∂ Z ψ (R, Z) = ∂ Z ψ (R, H i ) + ∂ 2 Z ψ (R, H i ) (Z -H i ) + ∂ 3 Z ψ (R, H i ) (Z -H i ) 2 2 + O Z 3 , = a 2 (∂ Z ψ 0 + ∂ Z ψ 1 + 2 ∂ Z ψ 2 + O ( 3 )) + a 2 (∂ 2 Z ψ 0 + ∂ 2 Z ψ 1 + 2 ∂ 2 Z ψ 2 + O ( 3 )) 1 8 R 4 + 1 16 2 R 6 + O 3 + a 2 2 ∂ 3 Z ψ 0 + O ( ) 1 8 R 4 + O 2 2 + O 5 . (3.11)
Since ∂ Z ψ (R, Z i ) = 0, the identification of the coefficients gives the following hierarchy of equations:

           ∂ Z ψ 0 (R, Z i ) = 0, ∂ Z ψ 1 (R, Z i ) = - 1 8 R 4 ∂ 2 Z ψ 0 (R, Z i ) , ∂ Z ψ 2 (R, Z i ) = - 1 8 R 4 ∂ 2 Z ψ 1 (R, Z i ) - 1 16 R 6 ∂ 2 Z ψ 0 (R, Z i ) - 1 128 R 8 ∂ 3 Z ψ 0 (R, Z i ) .
(3.12)

In the same way, the boundary conditions on the harmonics of ψ (and on their derivatives) can be found at the surface Γ j as follows:

           ψ 0 (R, Z j ) = 0, ψ 1 (R, Z j ) = 1 8β 3 R 4 ∂ Z ψ 0 (R, Z j ) , ψ 2 (R, Z j ) = - 1 128β 6 R 8 ∂ 2 Z ψ 0 (R, Z j ) , (3.13)
and,

           ∂ Z ψ 0 (R, Z j ) = 0, ∂ Z ψ 1 (R, Z j ) = 1 8β 3 R 4 ∂ 2 Z ψ 0 (R, Z j ) , ∂ Z ψ 2 (R, Z j ) = - 1 8β 3 R 4 ∂ 2 Z ψ 1 (R, Z j ) - 1 16β 5 R 6 ∂ 2 Z ψ 0 (R, Z j ) - 1 128β 6 R 8 ∂ 3 Z ψ 0 (R, Z j ) . (3.14)
Hence, each harmonic of ψ can be solved via Eq. (3.8) and the boundary conditions Eq. (3.10), Eq. (3.13), Eq. (3.12), and Eq. (3.14). In particular, the dominant order ψ 0 of ψ is solution of the following system:

                   ∂ 4 Z ψ 0 = 0, ψ 0 (R, Z i ) = 1 2 R 2 U sq i , ∂ Z ψ 0 (R, Z i ) = 0, ψ 0 (R, Z j ) = 0, ∂ Z ψ 0 (R, Z j ) = 0. (3.15)
As suggested by Cooley et al. [START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF] , ψ 0 can be written as follows:

ψ 0 = A 0 Z 3 + B 0 Z 2 + C 0 Z + D 0 , (3.16)
where A 0 , B 0 , C 0 , and D 0 coefficients obtained via the boundary conditions on ψ 0 . Therefore, ψ 0 is obtained by solving the following linear system:

      H 3 i H 2 i H i 1 3H 2 i 2H i 1 0 H 3 j H 2 j H j 1 3H 2 j 2H j 1 0             A 0 B 0 C 0 D 0       =       1 2 R 2 U sq i 0 0 0       , (3.17)
leading to:

                     A 0 = - R 2 H 3 U sq i B 0 = 3 2 R 2 U sq i H i + H j H 3 C 0 = -3R 2 U sq i H i H j H 3 D 0 = 1 2 R 2 U sq i (3H i -H j ) H j 2 H 3 (3.18)
where

H = H i -H j with H i = 1 + 1 2 R 2 and H j = - 1 2β
R 2 the dominant order of respectively Z i and Z j .

The other harmonics ψ 1 and ψ 2 can be solved in the same way as ψ 0 , by assuming their following expressions:

ψ 1 = - 1 10 Z 5 Υ(A 0 ) - 1 6 Z 4 Υ(B 0 ) + A 1 Z 3 + B 1 Z 2 + C 1 Z + D 1 , (3.19
)

ψ 2 = 1 280 Z 7 Υ 2 (A 0 ) + 1 120 Z 6 Υ 2 (B 0 ) - 1 10 Z 5 Υ(A 1 ) + 1 12 Υ 2 (C 0 ) - 1 6 Z 4 Υ(B 1 ) + 1 4 Υ 2 (D 0 ) + A 2 Z 3 + B 2 Z 2 + C 2 Z + D 2 , (3.20)
where According to Jeffrey et al. [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF] , only the coefficient A 2 is useful for computing the force on the spheres. Therefore, only the first three harmonics of ψ are exhibited. Furthermore, only ψ 0 is needed to compute the dominant order of the lubrication force, as shown in the next section below.

A 1 , B 1 , C 1 , D 1 , A 2 , B 2 ,

Expression of the lubrication force acting on the inner region.

The lubrication force is due to the hydrodynamic effect of the interaction of the given particle with a nearby obstacle. Thereby, the force is given by the action of the fluid stress on a particle surface when this particle and the interacting obstacle are isolated in an unbounded domain (flow field at rest far from the particles). The lubrication force is directly obtained by integrating the fluid stress σ acting on the particle surface:

σ = -pI + µ ∇u + (∇u) , (3.21)
with I the identity matrix.

Due to the symmetry of the flow, only the stress creating a force acting parallel to the axis of revolution (Oz) is relevant. Thus, the lubrication force on P i is

F lub i = Γ i (σ • n) • e z dS e z , (3.22)
with n the outgoing unit vector of the surface element dS of the particle P i as shown in et al. [START_REF] Happel | Low Reynolds Number Hydrodynamics[END_REF] . First, the fluid stress is projected in the intrinsic coordinate system (n, t, n × t) (Fig. 3.2), the normal fluid stress

(σ • n) becomes σ • n = (-p + 2µ∂ n u n ) n + µ (∂ t u n + ∂ n u t ) t, = -pn + 2µ [∂ n u n n + ∂ t u n t] + µ (∂ n u t -∂ t u n ) t, (3.23
)

with u = u n n + u t t = - 1 r ∂ t ψ n + 1 r
∂ n ψ t is the fluid velocity. Via straightforward calculations, Eq. (3.23) becomes:

σ • n = -pn + 2µ∇ 1 r ∂ t ψ + µ r Φ 2 (ψ) t, (3.24)
where

Φ 2 (•) = r∂ r 1 r ∂ r (•) +∂ 2 z (•) = ∇ 2 - 2 r ∂ r .
Thus, the lubrication force can be written as follows:

F lub i • e z = Γ i -p (n • e z ) + 2µ∇ 1 r ∂ t ψ • e z + µ r Φ 2 (ψ) (t • e z ) dS, = Γ i -p (∂ t r) + 2µ∂ z 1 r ∂ t ψ - µ r Φ 2 (ψ) (∂ n r) dS, = Γ i - 1 2r ∂ t r 2 p -r 2 ∂ t (p) + 2µ∂ z 1 r ∂ t ψ - µ r Φ 2 (ψ) (∂ n r) dS, = -π C ∂ t r 2 p ds + π C r 2 (∂ t p) ds -4µπ C ∂ z (∂ t ψ) ds -2µπ C (∂ n r) Φ 2 (ψ) ds, (3.25)
where dS = rdθds, Γ i = [0; 2π] × C with C = [0, R 0 ] a neighbourhood of the contact point A (see Fig. 3.1). To compute the lubrication force, the pressure term needs to be removed from Eq. (3.25). First, on the edges of C, r = 0 leading to

     C ∂ t r 2 p ds = 0, C ∂ z (∂ t ψ) ds = C ∂ t (ru r ) ds = 0. (3.26)
Furthermore, ∂ t p can be written as a function of ψ, by projection of the Stokes equation in the intrinsic coordinate system

∂ t p = ∇p • t, = µ∇ 2 u • t, = -µ (∇ × ζ) • t, = -µ∇ × Φ 2 (ψ) r (n × t) • t, = -µ∂ n (Φ 2 (ψ)) , (3.27)
where

ζ = ∇×u = (∂ n u t -∂ t u n ) (n × t)
F lub i = πµ C r 3 ∂ n Φ 2 (ψ) r 2 ds e z . (3.28)
The inner region, where the solution ψ of Eq. (3.6) is known, is defined by the neighborhood C = [0, R 0 ] of the contact point A (see Fig. 3.1). Before substituting the solution ψ in Eq. (3.28), ∂ n and ds need to be characterized in the inner region using the stretched coordinates.

The element of arc ds is obtained on the surface Γ i , by differentiation and a change of variable as follows:

ds Γ i = d arcsin r a , = d (arcsin (R √ )) , = √ dR √ 1 -R 2 .
(3.29)

Since only the first orders of ψ are known, the element of arc can be approximated by the first orders of the regular expansion of ds at R ∼ 0 (i.e. at the contact point A). Thus, the element of arc ds on the surface Γ i is

ds Γ i = √ 1 + 1 2 R 2 + 3 8 2 R 4 + O 3 dR. (3.30)
The normal derivative ∂ n on the surface Γ i is obtained by differentiation of the outgoing normal unit vector n. Via elementary geometrical relations (Fig. 3.2), the outgoing normal vector on Γ i is given by

n Γ i = R √ e r - √ 1 -R 2 e z . (3.31)
Thus, the normal derivative ∂ n is given by

∂ n Γ i = R∂ R - 1 √ 1 -R 2 ∂ Z . (3.32)
As noticed by Jeffrey et al. [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF] , only the first orders of ∂ n are relevant. Thereby, the regular expansion of ∂ n at R ∼ 0 gives the following approximation of the normal derivative

∂ n Γ i = R∂ R - 1 1 - 1 2 R 2 - 1 8 R 4 2 + O 3 ∂ Z . (3.33)
Via straightforward calculations and simplifications, using Eqs. (3.30) and (3.33), the first orders of the lubrication force in the inner region is

F lub,in i • e z πµaU sq i = 1 aU sq i C r 3 ∂ n Φ 2 (ψ) r 2 ds, = 1 a 2 U sq i C R 3 ∂ n 1 R 2 ∂ 2 R - 1 R ∂ R + 1 ∂ 2 Z (ψ) ds, = 1 f -1 + f 0 + f 1 + O 2 ,
(3.34)

3.1. THEORY CHAPTER 3. LUBRICATION CORRECTION MODELS with                          f -1 = - 1 U sq i C R∂ 3 Z ψ 0 Z=Z i dR, f 0 = - 1 U sq i C R∂ 3 Z ψ 1 -∂ R ∂ Z ψ 0 + 2R∂ 2 Z ψ 0 + R∂ Z ∂ 2 R ψ 0 -R 2 ∂ R ∂ 2 Z ψ 0 Z=Z i dR, f 1 = - 1 U sq i C R∂ 3 Z ψ 2 -∂ R ∂ Z ψ 1 + 2R∂ 2 Z ψ 1 + R∂ Z ∂ 2 R ψ 1 -R 2 ∂ R ∂ 2 Z ψ 1 -3∂ R ψ 0 + 3R∂ 3 R ψ 0 -R 2 ∂ 3 R ψ 0 - R 4 2 ∂ R ∂ 2 Z ψ 0 Z=Z i dR. (3.35)
Therefore, the dominant order (in ) of the lubrication force in the inner region is given by

F lub,in i πµaU sq i = - 6R 4 0 4H 2 0 1 e z + O (1) , (3.36)
with

H 0 = 1 + 1 + β 2β R 2 0 .
When two particles are close to each other ( 1), the lubrication force becomes the dominant force. Furthermore, the magnitude of the lubrication stress σ lub,in i is more and more concentrated at the contact point [START_REF] Serayssol | The elastohydrodynamic collision of two spheres[END_REF] as the gap becomes narrower (see Fig. 3.3 and Eq. 3.37). Thereby, lubrication corrections are critical and complex to perform in numerical simulations.

σ lub,in i ∼ 3U sq i µ β a(1 + β)

Computation of the total lubrication force.

To obtain the lubrication force on the whole surface, the lubrication force on the outer region needs to be added to the contribution from the inner region Eq. (3.36). This can be done by solving Eq. (3.28) with ψ solution of a Stokes problem in the outer region and matching this solution to the solution in the inner region [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF][START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF] . Another approach is to match the inner solution Eq. (3.36) with a solution on the outer region obtained via a numerical simulation [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF] . At the end, both methods give the dominant order of the total lubrication force as follows: [START_REF] Hovnanian | An accurate cartesian method for incompressible flows with moving boundaries[END_REF] 1 e z + O (1) .

F lub i 6πµaU sq i = - β 2 (1 + β)
(3.38)

Shearing motion of rigid spheres.

As for the squeezing motion, two smooth particles of radii a and βa, with β > 0, are immersed in an infinite domain of fluid. The particle P i is now moving near the stationary particle P j at the velocity U i = U sh i e x (see Fig. 3.1). The surfaces of P i and P j are denoted respectively Γ i and Γ j . The gap between the two particles is denoted by a (0 < 1) such that the gap Reynolds number is small (Re = ρU sh i a /µ 1).

According to Kim et al. [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF] , symmetries in the flow field in the gap suggest that (u, p) can be written in the cylindrical coordinate system (e r , e θ , e z ) as follows:

             u r = U sh i U (r, z) cos (θ) , u θ = U sh i V (r, z) sin (θ) , u z = U sh i W (r, z) cos (θ) , p = µU sh i a P (r, z) cos (θ) , (3.39)
where the angle θ is defined such that e r = cos(θ)e x + sin(θ)e y , and U , V , W , and P are unknown functions such that the velocity and pressure fields are solutions of the Stokes problem Eq. (3.3). With a velocity and pressure fields defined as Eq. (3.39), the Stokes problem Eq. (3.3) becomes:

                     0 = ∂ r U (r, z) + U (r, z) + V (r, z) r + ∂ z W (r, z) , 1 a ∂ r P (r, z) = L 2 0 (U (r, z)) - 2 r 2 (U (r, z) + V (r, z)) , - 1 a P (r, z) r = L 2 0 (V (r, z)) - 2 r 2 (U (r, z) + V (r, z)) , 1 a ∂ z P (r, z) = L 2 1 (W (r, z)) , (3.40)
3.1. THEORY CHAPTER 3. LUBRICATION CORRECTION MODELS with L 2 m = ∂ r 2 + 1 r ∂ r -m 2 r 2 + ∂ z 2 .
The non-slip boundary conditions U Γ i = U sh i e x and U Γ j = 0 are also projected in the cylindrical coordinate system and become The resolution of the pressure and velocity fields is similar to the resolution of the Stokes stream-function for the squeezing motion. The unknowns P , U , V , and W are assumed [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF] to be regular expansions in as follows:

                       u r Γ i = 1, u θ Γ i = -1, u z Γ i = 0, u r Γ j = 0, u θ Γ j = 0, u z Γ j = 0.
           P (R, Z) = -3/2 P 0 (R, Z) + -1/2 P 1 (R, Z) + O 1/2 , U (R, Z) = U 0 (R, Z) + U 1 (R, Z) + O ( 2 ) , V (R, Z) = V 0 (R, Z) + V 1 (R, Z) + O ( 2 ) , W (R, Z) = 1/2 W 0 (R, Z) + 3/2 W 1 (R, Z) + O 5/2 . (3.42)
These assumptions on the velocities and pressure are motivated by the consistency of boundary conditions. According to Kim et al. [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF] , boundary conditions Eq. (3.41) imply that U and V are O ( 0 ), so from the equation of continuity in Eq. (3.40), W is scaled by O 1/2 . Then, the scaling for the pressure follows directly the governing equations as

P = O -3/2 .

Substitution of the asymptotic expansion Eq. (3.42) into the governing equations

Eq. (3.40), yields a hierarchy of equations in order of . To make the identification of coefficients for each order of easier, Eq. (3.40) is re-written using stretched coordinates:

                     0 = 1 √ ∂ R U + U + V R √ + 1 ∂ Z W, 1 √ ∂ R P = 1 √ ∂ R 2 U + 1 R ∂ R U + 1 2 ∂ Z 2 U - 2 R 2 (U + V ) , - P R √ = 1 √ ∂ R 2 V + 1 R ∂ R V + 1 2 ∂ Z 2 V - 2 R 2 (U + V ) , 1 ∂ Z P = 1 √ ∂ R 2 W + 1 R ∂ R W - 1 R 2 W + 1 2 ∂ Z 2 W. (3.43)
Hence, the lowest order in for each line of Eq. (3.43) is

               0 = ∂ R U 0 + U 0 + V 0 R + ∂ Z W 0 , ∂ R P 0 = ∂ Z 2 U 0 , - P 0 R = ∂ Z 2 V 0 , ∂ Z P 0 = 0, (3.44)
and, the next order for each line is

                   0 = ∂ R U 1 + U 1 + V 1 R + ∂ Z W 1 , ∂ R P 1 = ∂ R 2 U 0 + 1 R ∂ R U 0 + ∂ Z 2 U 1 - 2 R 2 (U 0 + V 0 ) , - P 1 R = ∂ R 2 V 0 + 1 R ∂ R V 0 + ∂ Z 2 V 1 - 2 R 2 (U 0 + V 0 ) , ∂ Z P 1 = ∂ Z 2 W 0 . (3.45)
To close the systems of equations (3.44) and (3.45), the boundary conditions on the harmonics of P , U , V , and W need to be calculated. In the same way, as for the squeezing motion, the Taylor expansions of P , U , V , and W at the particle surfaces are taken and used with the boundary conditions Eq. (3.41) to obtain the boundary conditions on the harmonics.

For instance, on the surface Γ i , the Taylor expansion of U , Eq. (3.42) at Z = H i is given by

U (R, Z) = U (R, H i ) + ∂ Z U (R, H i ) (Z -H i ) + O Z 2 , = U 0 + U 1 + O 2 + ∂ Z U 0 + ∂ Z U 1 + O 2 1 8 R 4 + 1 16 2 R 6 + O 3 + O 3 , = U 0 + U 1 + 1 8 ∂ Z U 0 + O 2 .
From the boundary conditions Eq. (3.41) on U , U (R, Z A ) ∼ 1. Therefore, by identifying coefficients of the Taylor expansion of U at Z = H i evaluated at Z = Z i , the following hierarchy of equations is obtained

U 0 (R, Z A ) = 1, U 1 (R, Z A ) = -1 8 R 4 ∂ Z U 0 . (3.46)
In the same way, the following boundary conditions are obtained for U , V , W , and P on

                       U 0 (R, Z i ) = 1, V 0 (R, Z i ) = -1, W 0 (R, Z i ) = 0, U 0 (R, Z j ) = 0, V 0 (R, Z j ) = 0, W 0 (R, Z j ) = 0, (3.47)
and The fourth line of Eq. (3.44) implies that P 0 is only a function of R. Therefore, the second and the third lines of Eq. (3.44) can be integrated in the gap, by separating variables leading to the following parabolic profile for U 0 and V 0

                                     U 1 (R, Z i ) = - 1 8 R 4 ∂ Z U 0 (R, Z i ) , V 1 (R, Z i ) = - 1 8 R 4 ∂ Z V 0 (R, Z i ) , W 1 (R, Z i ) = - 1 8 R 4 ∂ Z W 0 (R, Z i ) , U 1 (R, Z j ) = 1 8β 3 R 4 ∂ Z U 0 (R, Z j ) , V 1 (R, Z j ) = 1 8β 3 R 4 ∂ Z V 0 (R, Z j ) , W 1 (R, Z j ) = 1 8β 3 R 4 ∂ Z W 0 (R, Z j ) .
     U 0 (R, Z) = - 1 2 P 0 (Z -Z j ) (Z i -Z) + Z -Z j H , V 0 (R, Z) = 1 2 P 0 R (Z -Z j ) (Z i -Z) - Z -Z j H , (3.49)
with

H (R) = Z i -Z j .
Then, the integration between Z i and Z j of the first line of Eq.

(3.44) gives

W 0 (R, Z i ) -W 0 (R, Z j ) = - H 0 ∂ R U 0 + U 0 + V 0 R dζ, = - H 0 1 R ∂ R (RU 0 ) + V 0 R dζ, = - 1 R ∂ R H 0 (RU 0 ) dζ -R∂ R Z i - 1 R H 0 V 0 dζ, (3.50)
with ζ = Z + Z j . The change of the order of integration and differentiation in the last line is obtained by applying the Leibniz rule, since Z i depends on R. Furthermore, W 0 (R, Z i ) -W 0 (R, Z j ) = 0 as u z = 0 at the particle surfaces. By substituting the expression of U 0 and V 0 from Eq. (3.49) in Eq. (3.50), and via the approximations

∂ R Z i ≈ R and ∂ R H ≈ 1 + 1 β R
, the Reynolds equation is obtained as

R 2 P 0 + R + 3 1 + 1 β R 3 H P 0 -P 0 = -6 1 - 1 β R 3 H 3 . (3.51)
The unique solution [START_REF] Kim | Microhydrodynamics: Principles and Selected Applications[END_REF][START_REF] O'niell | On the slow motion of a sphere parallel to a nearby plane wall[END_REF] of Eq. (3.51) with lim R→+∞ P 0 = 0 is

P 0 = 6R 5H 2 1 - 1 β . (3.52)
Hence, U 0 and V 0 in Eq. (3.49) are known anywhere within the gap:

                 U 0 (R, Z) = - 1 2 dP 0 dR (Z -Z j ) (Z i -Z) + Z -Z j H , P 0 (R) = 6R 5H 2 1 - 1 β , V 0 (R, Z) = - 1 2 P 0 R (Z -Z j ) (Z i -Z) - Z -Z j H . (3.53)
According to O'Niell et al. [START_REF] O'niell | On the slow motion of a sphere parallel to a nearby plane wall[END_REF] , P 1 , U 1 , V 1 , and W 1 can be found by means of Eq. (3.45) after substitution of U 0 , V 0 , W 0 , and P 0 in Eqs. (3.45) and (3.48). However, only U 0 , V 0 , and P 0 are needed to get an expression of the dominant order of the lubrication force and torque.

Expression of the lubrication force and torque acting on the inner region.

As for the squeezing motion, the lubrication force and torque are directly obtained by integration of the fluid stress σ (see Eq. (3.21)) over the particle surface Γ i . Due to the symmetry of the flow, only the component of the force according to the (Ox) is non-zero.

Hence, the lubrication shear force and torque acting on an element of surface dS of Γ i , are given by

       F lub i = Γ i (σ • n) • e x dS e x , T lub i = Γ i (σ • n) × e y dS.
(3.54)

By substitution of (u, p), expressed as Eq. (3.39), and straightforward simplifications, the lubrication force and torque become

                         F lub i = Γ i U sh i µ cos(χ) -∂ z U cos 2 (θ) -∂ r W cos 2 (θ) + ∂ z V sin 2 (θ) + sin(χ) -P a cos 2 (θ) + 2∂ r U cos 2 (θ) -∂ r V sin 2 (θ) dS e x , T lub i = Γ i U sh i µa cos 2 (χ) ∂ z U cos 2 (θ) + ∂ r W cos 2 (θ) -∂ z V sin 2 (θ) -sin 2 (χ) (∂ z U cos 2 (θ) + ∂ r W cos 2 (θ)) -cos(χ) sin(χ) 4∂ r U cos 2 (θ) -∂ r V sin 2 (θ) dS e y , (3.55)
in the cylindrical coordinate system (e r , e θ , e z ). The angle χ is defined such that n • e r = sin(χ) and n • e z = -cos(χ) (see Fig. 3.2). Before expending further Eq. (3.55) with the pressure and velocity expressions of Eq. (3.42), the force and torque are integrated according to θ and expressed using the stretched coordinates Eq. (3.2). Hence, the force and torque become

                               F lub i U sh i µaπ = χ 0 0 cos(χ) - 1 √ ∂ R W + 1 (∂ Z V -∂ Z U ) + sin(χ) -P + 1 √ (2∂ R U -∂ R V ) sin(χ)dχ e x , T lub i U sh i µa 2 π = χ 0 0 cos 2 (χ) 1 (∂ Z U -∂ Z V ) + 1 √ ∂ R W -sin 2 (χ) 1 ∂ Z U + 1 √ ∂ R W -cos(χ) sin(χ) 1 √ (4∂ R U -∂ R V ) sin(χ)dχ e y .
(3.56)

The parameter χ 0 ∈ [0; π] defines the inner region, where (u, p) are solutions of Eq. (3.3).

To use the same definition of the inner region as for the squeezing motion, the variable χ is changed as follows:

     χ = arcsin(R √ ), dχ = √ dR √ 1 -R 2 = 1 + 1 2 R 2 + O 2 dR. (3.57)
Hence, the inner region is defined by C = [0; R 0 ] and the Eq. (3.56) become

                                 F lub i U sh i µaπ = - R 0 0 √ ∂ R W + ∂ Z U -∂ Z V RdR e x + R 0 0 -P √ + 2∂ R U -∂ R V R 2 dR √ 1 -R 2 e x , T lub i U sh i µa 2 π = R 0 0 (∂ Z U -∂ Z V ) + √ ∂ R W R √ 1 -R 2 dR e y - R 0 0 ∂ Z U + √ ∂ R W R 3 dR √ 1 -R 2 e y - R 0 0 [4∂ R U -∂ R V ] R 2 dR e y .
(3.58)

Since the gap is small ( 1) and we want only the first orders of the lubrication force and torque, the squared roots in Eq. (3.58) are approximated as follows:

                                 F lub i U sh i µaπ = - R 0 0 √ ∂ R W + ∂ Z U -∂ Z V RdR e x + R 0 0 -P √ + 2∂ R U -∂ R V R 2 1 + 1 2 R 2 + O 2 dR e x , T lub i U sh i µa 2 π = R 0 0 (∂ Z U -∂ Z V ) + √ ∂ R W R 1 - 1 2 R 2 + O 2 dR e y - R 0 0 ∂ Z U + √ ∂ R W R 3 1 + 1 2 R 2 + O 2 dR e y - R 0 0 [4∂ R U -∂ R V ] R 2 dR e y . ( 3
                                     F lub i U sh i µaπ = R 0 0 [-P 0 R + ∂ Z V 0 -∂ Z U 0 ] RdR e x + R 0 0 -P 0 2 R 3 -P 1 R + 2∂ R U 0 -∂ R V 0 -∂ R W 0 + ∂ Z V 1 -∂ Z U 1 RdR e x + O ( 2 ) , T lub i U sh i µa 2 π = R 0 0 [∂ Z U 0 -∂ Z V 0 ] RdR e y + R 0 0 - 3 2 R 2 ∂ Z U 0 + 1 2 R 2 ∂ Z V 0 -R (4∂ R U 0 -∂ R V 0 ) + ∂ Z U 1 -∂ Z V 1 (1 -R 2 ) ∂ R W 0 ] RdR e y + O ( 2 ) . (3.60)
In particular, the dominant orders (in ) of the lubrication force and torque in the inner

F lub,in i aµπU sh i = R 0 0 [-P 0 R + ∂ Z V 0 -∂ Z U 0 ] RdR e x + O ( ) , (3.61) T lub,in i a 2 µπU sh i = R 0 0 [∂ Z U 0 -∂ Z V 0 ] RdR e y + O ( ) , (3.62)
with U 0 , V 0 , and P 0 given by Eq. (3.53).

As for the squeezing motion, the total lubrication force and torque can be found by matching [START_REF] O'neill | Asymmetrical slow viscous fluid motions caused by the translation or rotation of two spheres. part ii: Asymptotic forms of the solutions when the minimum clearance between the spheres approaches zero[END_REF][START_REF] Jeffrey | The forces and couples acting on two nearly touching spheres in low-Reynolds-number flow[END_REF] the inner and outer solution of the Stokes problem (3.3):

F lub i aµπU sh i = - 24β (2 + β + 2β 2 ) 15 (1 + β) 3 ln ( ) e x + O ( ) , (3.63)
T lub i a 2 µπU sh i = - 8β (4 + β) 10 (1 + β) 2 ln ( ) e y + O ( ) , (3.64)
For a given velocity U i = U sq i = U sh i , the lubrication force created by the squeezing motion is dominant:

F lub,sh i F lub,sq i ∼ 1 ln( ) .
Therefore, lubrication effects created by the shearing motion are often neglected in lubrication correction models (see Sec. 3.2).

Costa Lubrication Correction Model (CLM).

The CLM [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] is a two-parameter model which corrects the normal component of the lubrication force on a spherical particle. The correction is made by adding ∆F lub i = ∆F lub i e n to the computed hydrodynamic force:

∆F lub i 6πµaU sq i =        λ( ∆x ) -λ( ), col ≤ < ∆x , λ( ∆x ) -λ( col ), 0 ≤ < col , 0, otherwise, (3.65)
where e n is defined as represented in Fig. 3.6. The Stokes amplification factor [START_REF] Jeffrey | Low Reynolds number flow between converging spheres[END_REF] λ is defined for the lubrication interaction between a sphere and a wall as λ pw , and for interaction between two spheres as λ pp : (3.66)

     λ pw ( ) = 1 - 1 5 ln( ) - 1
Hence, the total hydrodynamic force is given by

F hyd i = Γ i (σ • n) dS + ∆F lub i . (3.67)
The value of the parameter ∆x is determined by simulating the slow approach of a sphere toward a wall or a second particle, for a given grid cell resolution h. Figures 3.4

and 3.5 represent the total hydrodynamic force acting on a particle approaching a wall or a particle, respectively. Simulations without lubrication correction are compared to the analytical solutions given by Brenner [START_REF] Brenner | The slow motion of a sphere through a viscous fluid towards a plane surface[END_REF] and Cooley et al. [START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF] . The parameter ∆x is defined as the largest value of such that for ≤ ∆x the hydrodynamic force from the simulation without lubrication correction no longer matches the analytical solution.

The CLM parameter for Fig. Figure 3.4: Total hydrodynamic force as a function of in the case of a single particle approaching a solid wall. Simulations results are compared against the analytical solution of Brenner [START_REF] Brenner | The slow motion of a sphere through a viscous fluid towards a plane surface[END_REF] (dashed line). All simulations are run with h = 1/40. The LLCM is activated on all smaller than lub = 2 ∆y a .

Figure 3.5: Total hydrodynamic force as a function of in the case of a single particle approaching another particle. Simulations results are compared against the analytical solution of Cooley [START_REF] Cooley | On the slow motion generated in a viscous fluid by the approach of a sphere to a plane wall or stationary sphere[END_REF] (dashed line). All simulations are run with h = 1/40. The LLCM is activated on all smaller than lub = 2 ∆y a .

The Local Lubrication Correction Model (LLCM).

The lubrication force F lub i and torque T lub i acting on P i are given by

F lub i = j∈ 1,N \{i} F lub i,j + F lub i,wall ,
T lub i = j∈ 1,N \{i} T lub i,j + T lub i,wall , (3.68)
where F lub i,j and T lub i,j are the lubrication force and torque acting on P i , created by the interaction of P i and P j . The lubrication force and torque created by the interaction of the particle P i with a wall are denoted F lub i,wall and T lub i,wall , and are equivalent to the asymptotic case β → +∞.

From the lubrication theory detailed in Sec. 3.1, F lub i,j and T lub i,j are given by

                                                   F lub i,j = F lub,sq i,j + F lub,sh i,j 1 [ col , lub ]( ) , T lub i,j = T lub,sh i,j 1 [ col , lub ]( ) , F lub,sq i,j πµaU sq i,j = - 6R 4 0 4H 2 0 1 e n , F lub,sh i,j πµaU sh i,j = R 0 0 [-P 0 R + ∂ Z V 0 -∂ Z U 0 ] RdR e t ,
T lub,sh i,j

πµa 2 U sh i,j = R 0 0 [∂ Z U 0 -∂ Z V 0 ] RdR e n × e t , U sq i,j = U sq i -U sq j , U sh i,j = U sh i -U sh j + (ω i -ω j β) a, (3.69)
with

             H 0 = 1 + 1 2 aαR 2 0 , α = β + 1 aβ , R 0 = 1 √ 1 -(1 + -lub ) 2 .
(3.70)

The projected relative velocity of the two particles on the directions e n and e t are U sq i,j and U sh i,j (see Fig. 3.6).

The indicator function 1 [ lub , col ] ( ) is non-zero and equals to 1 only if col ≤ ≤ lub .

Hence the lubrication model is switched off when the hydrodynamics is fully resolved ( ≥ lub ), and when the gap disappears ( ≤ col ).

Figure 3.6: Sketch of the two interacting particles with the notations used to evaluate F lub i,j and T lub i,j .

THE LLCM CHAPTER 3. LUBRICATION CORRECTION MODELS

The analytic expressions of F lub,sh i,j and T lub,sh i,j

are not as simple as the expression of F lub,sq i,j

. Therefore, we choose to compute the shearing components by a numerical integration. Furthermore, numerical integration will be necessary if we want to add lower orders of the lubrication forces and torques. The force F lub,sh i,j

and the torque T lub,sh i,j are integrated on n lub sub-sets of the lubrication region [0, R 0 ]:

F lub,sq i,j πµaU sq i,j = n lub k=1 R R R k [-P 0 R + ∂ Z V 0 -∂ Z U 0 ] RdR e t ,
T lub,sh i,j

πµa 2 U sh i,j = n lub k=1 R R R k [∂ Z U 0 -∂ Z V 0 ] RdR e n × e t , (3.71)
where

R R R k = k -1 n lub R 0 , k n lub R 0 .
The local lubrication correction model (LLCM) is then built such that the total hydrodynamic force and torque acting on P i are approximated by

     F hyd i ≈ F solv i + F lub i , T hyd i ≈ T solv i + T lub i . (3.72)
As the lubrication corrections F lub i,j and T lub i,j are the dominant orders of the hydrodynamics acting on the inner region of P i , when the flow in the gap is in the Stokes regime (see Sec.

Re gap i,j = ρa i U i -U j µ
remains small for all P j in interaction with P i (a i is the radius of P i ), during the simulations. In particular, the Reynolds number Re lub needs to be moderate (typically Re lub < 10 3 ):

Re lub = ρ lub Q max µ , (3.73)
with

Q max = max t max (i,j)∈ 1,N 2 (a i U i (t) -U j (t) ) .
This constraint also underlines the limitation of the LLCM to moderate Reynolds number flows. Indeed, inertia effects of the fluid in the gap are not corrected by the LLCM.

Another limitation of the LLCM concerns the many-body interactions, which refers to the 3.3. THE LLCM hydrodynamics action on a particle generated by nearby particles in interactions. As the lubrication corrections are made in an additive pairwise-fashion (see Eq. (3.68)), the only many-body interactions, present in the simulated flow, are the ones resolved by the numerical method (included in F solv i ).

For ellipsoidal particles, the lubrication corrections are evaluated on the virtual spheres introduced in Sec. 2.1.2.4. The approximation of the particle surfaces by the spheres is an additional source of uncertainty from the LLCM. Therefore, the size of the grid elements might need to be reduced for ellipsoidal particles with extreme aspect ratios to limit the lubrication area to contact point neighborhoods where the virtual spheres fit the most the particle surfaces.

Chapter 4

Numerical simulations of spherical particles Interactions between a spherical particle and an obstacle (wall or another particle) have been extensively investigated both numerically and experimentally. The central symmetry makes spherical particles the easiest particle shape to study experimentally. Among the numerous benchmarks available in the literature, three cases have been chosen to validate the LLCM. First, the lubrication correction for a particle interacting with a wall is simulated when the particle is moving normally toward the wall at several impact Stokes number. By this simple case, the dominant component of the lubrication force is validated with or without rebound of the particle. Subsequently, the wall is inclined to validate the shearing component of the lubrication forces. Finally, particle-particle interactions are simulated for two particles in a sheared flow.

Falling particle onto a wall.

A single particle is immersed in a domain [4d] 3 , with d = 2a the particle diameter, uniformly meshed with cubic elements of size ∆x = ∆y = ∆z = hd. Periodic boundary conditions are considered on the lateral faces of the domain. A no-slip boundary condition is imposed at the top and bottom walls (y constant). Simulations have been performed on a [8d] 3 domain (with h = 1/40) and have given identical solutions than on the smaller domain. The fluid is initially at rest and the particle is dropped without initial velocity such that the gap size from the bottom wall is given by a init as shown Fig. 4.1. The gravity field g acts on the y-direction.

To have an experimental reference for comparison, we chose the same configuration as one of the two cases from Harada et al. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] . Particle and fluid properties are presented in Tab.

       Re d = ρU T d µ ≈ 25.7, St d = ρ p U T d 9µ ≈ 3.27. (4.1)
A posteriori computations show that the maximum gap Reynolds number is

Re gap i,wall = ρa i max ( U i ) µ ≈ 0.92,
for ≈ 0.31. Therefore, the lubrication forces form the LLCM should reasonably approx-

The simulation using CLM (Fig. 4.2) has been run with ∆x = 0.06 and h = 1/40.

Aside the computation of the hydrodynamics Eq. (3.67) for the CLM, the resolution of the fluid and particle dynamics is identical to the approach used with the LLCM. For the Harada case, CLM results are in good agreement with experimental measurement. As shown in Fig. 4.2, the Pearson's correlation coefficient is equal to 0.9661,

E 2 (U CLM ) = 2.287 • 10 -3 m • s -1 , E ∞ (U CLM ) = 1.229 • 10 -2 m • s -1 .
The accuracy of the LLCM is comparable to the CLM. The simulated velocities remain close to experimental data (Pearson's correlation coefficient equals to 0.9633,

E 2 (U LLCM ) = 2.385 • 10 -3 m • s -1 , and E ∞ (U LLCM ) = 1.232 • 10 -2 m • s -1
). The computational time is about 10% larger with LLCM than CLM, without considering the computational cost of the tabulation. However, the computational time of the tabulation needed for the CLM can be significant. Furthermore, the LLCM can be extended to non-spherical particles while the CLM will require the tabulation of several parameters at a considerable CPU cost. The computational cost of the tabulation is also a CLM's limitation for polydisperse flows of spherical particles. Therefore, the LLCM is a reliable and efficient alternative to tabulated methods for complex particle flows.

Most of this additional cost comes from the identification of the set L L L(P ) (see definition in Sec.2.3.1) for each surface element of the particle. A very basic approach has been used here which consists in computing the distance of each surface element to the wall. More cost-effective methods (optimized storage of the surface elements, for instance) can be found and would probably reduce the cost difference between the LLCM and the CLM.

Due to the lack of experimental measurements of the total hydrodynamic forces for a falling particle, total hydrodynamic forces obtained by numerical simulations are compared to an analytic model (model H) introduced by Harada et al. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] . The model H is based on the Stokesian dynamic, where hydrodynamic force

F MH i = F hyd i + F lub i on the
particle is modeled by the dominant order of the lubrication force f l , the added-mass force f a , and the Basset history force f Ba .

The dominant order of the lubrication force is obtained by integration of the interstitial pressure between the particle and the wall [START_REF] Adams | The cohesive forces between particles with interstitial liquid[END_REF] :

f l = 6πµa U sq i . (4.2)
Since the motion of the particle is unsteady, the added mass force changes near the wall as follows [85][86] :

       f a = m dU i dt + 1 2 dm dt U i , m = 2 3 πρa 3 1 + ∞ i=0 3a (i+1) f 0 f 1 ...f i , (4.3)
where f i is recursively defined such that f 0 = 2a( + 1) and f i = f 0 -a 2 /f i-1 . The Basset history force for a spherical particle is given by

f Ba = 6a 2 √ πρµ t -∞ dU i dt dt √ t -t . (4.4)
Hence, the total hydrodynamic force of the model H is

F MH i = f Ba + f a + f l .
Numerical simulations have shown that velocities computed with the LLCM tends to experimental measurements as h decreases. Table 4.2 shows that the order of the numerical method used is close to a second order in L 2 and L ∞ on velocity when no lubrication corrections are applied. The particle dynamics for h = 1/20 is only fully resolved for large , where the lubrication forces are not dominant. Therefore, errors on the particle position and velocity are unexpectedly small for h = 1/20 compared to smaller h. Comparisons with Harada's measurements also show that the LLCM gives better results than model H for h ≤ 1/40, even before any lubrication corrections (for the model H: Pearson's correlation coefficient equals to 0.9418,

E 2 (U MH ) = 3.2 • 10 -3 m • s -1 , E ∞ (U MH ) = 1.3 • 10 -2 m • s -1 ).

Grid sensitivity analysis.

Using the same configurations as in Table 4.1, simulations have been performed with four different background grid resolutions h = 1/20, h = 1/40, h = 1/80, and h = 1/100.

The total hydrodynamic forces obtained are compared to the model H (Fig. 4.3).

Far from the wall ( > 0.1), higher grid resolutions (h ≤ 1/80) improve the resolution of the total hydrodynamic forces. Indeed, as shown in the Fig. 4.3, the hydrodynamic forces tend to the hydrodynamic forces of the model H which fits experimental velocities (Fig. 4.2). As the particle goes closer to the wall ( < 0.1), simulated hydrodynamic forces are lower than the solution of the model H, even before any lubrication correction is applied (see h = 1/80 and h = 1/100). Starting from the same initial conditions, the lower hydrodynamic forces induce higher particle velocities (than model H velocities) close to the wall. Therefore, the response of the LLCM is stronger close to the wall with steeper and is ignored here for the simulations using finer grid mesh resolution than h = 1/40.

Velocity h U 1/100 -U h 2 U 1/100 -U h ∞ order L 2 order L ∞ 1/20
Therefore, the global order of convergence on velocity and position of the LLCM could be closer to a second-order, since the computation of the hydrodynamics rely more and more (as h decrease) on a numerical method without lubrication correction, which is a second order of convergence on position and velocity (see Table 4.2). In addition, the LLCM is more accurate if lub is adapted to the grid resolution h than if lub is fixed. For instance, with h = 1/100, the relative errors E 2 (U LLCM ) and E ∞ (U LLCM ) are increased by 19% and 12% by using lub = 0.1 instead of lub = 2∆y a (= 0.04). The value of lub depends on the spatial discretization of the Navier-Stokes equations and needs to be adapted to the numerical method. This configuration can be easily done by running a sensitivity analysis from lub , as shown in this section.

Analysis of the LLCM reliability during a collision.

To complete our model reliability assessment, several simulations had been performed at higher Stokes numbers and compared to experimental data from Joseph et al. [START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] .

The conservation of the energy by the LLCM and the collision model have been checked by simulating the normal collision of a steel particle with a wall (see Fig. [START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] . By definition, ξ n = -U R /U T where U R is the particle rebound velocity and U T is the terminal velocity of the particle before the of a single particle impacting a wall in respect of the particle Stokes number at the impact. Filled markers represent results obtained using the LLCM using the same experimental setup than Joseph et al. [START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] (data represented by hallow markers). The black curve represents a correlation made on experimental data proposed by Legendre et al. [START_REF] Daniel | Experimental study of a drop bouncing on a wall in a liquid[END_REF] When the fluid does not influence the dynamics of the particle while colliding with an obstacle (non-viscous fluids for instance), the effective coefficient ξ n is maximum and is In this section, the bouncing motion of a single particle colliding onto a wall is simulated. This experiment is a common benchmark for confirming that lubrication and collision model return realistic bouncing dynamics by comparing simulation results against Gondret et al. [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] measurements. Furthermore, small differences in rebound velocity are amplified after its temporal integration, and therefore more noticeable in the particle trajectory [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] .

The particle of diameter d is immersed in a domain [8d, 30d, 8d] and falls under its own weight (g = 9.81m • s -1 ) onto a planar surface (Fig. 4.1). The domain is mapped using a uniform Cartesian mesh with grid size ∆x = ∆y = ∆z = hd. Dry and wet collisions have been simulated using the particle and fluid properties contained in Tab. 4.3 and 4.4. The particle is dropped with an initial velocity close to the particle sedimentation velocity to ease the convergence of the velocity to its terminal value before the particle impact the wall. and position to Gondret et al. [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] measurements for dry and wet collisions, respectively.

4.2 Oblique impact of a particle on a wall.

To validate the tangential component of the lubrication force, the oblique impact of a particle on a wall is the last case to validate our model. The same experimental set-up as Joseph and Hunt [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] , for oblique particle-wall collisions in the air and aqueous solutions, has been used.

A single particle is immersed in a domain [19d, 5d, 3d] using a similar set-up as described in Sec. 4.1. The fluid is initially at rest and the particle is dropped without initial velocity as shown in Fig. 4.11. To reproduce the experiment of Joseph and Hunt [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] the particle is linked to O by a virtual string of length L = 10 cm. The tension of the string is a virtual force F T i collinear to

e 0 = X O -X i X O -X i
is added to the particle dynamic Eq. (2.18). The signed value of F T i • e 0 is found such that: [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] experiments. and rebound Ψ out = tan (ψ out ) angles obtained from oblique collisions between steel and glass spheres in the air or water. In practice, Ψ in and Ψ out have been assessed using the following equations:

U i • e 0 = 0, ∀ X O -X i ≥ L + a.

Glass

Ψ in = - V C,in • e x V C,in • e y , Ψ out = V C,out • e x V C,out • e y , (4.6)
with V C,in and V C,out the particle velocities at the contact point,

V C = U i -a Ω i × e y ,
just before and just after collision, respectively.

As highlighted by Joseph and Hunt [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] , piezoviscous effects in wet collision modify the coefficient of friction. Therefore, the modified (wet) friction coefficient (see Tab. For wet collisions (Fig. 4.12 bottom), hydrodynamic effects are no longer negligible.

Therefore, the accuracy of the lubrication correction can be assessed. The overall numerical simulations are in good agreement with the experimental data for the entire range of incidence angles. Hunt [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] experimental measurements (by circles).

below, the wall velocity U wall is set such that the particle Reynolds number is Re d = 0.75.

As exhibited by Yan et al. [START_REF] Morris | Hydrodynamic interaction of two particles in confined linear shear flow at finite reynolds number[END_REF] , the initial relative positions of both particles have an influence on the particle trajectories. In this section, Yan's test cases are simulated in three-dimensions using spherical particles. distinguishes repelling and bypass trajectories is found to be independent of the Reynolds number [START_REF] Morris | Hydrodynamic interaction of two particles in confined linear shear flow at finite reynolds number[END_REF] . This δ y is lower than the critical value (δ y = 2.9d) found by Yan et al.. This difference can be explained by the grid resolution used for this simulation (h = 1/30), and due to resource limitations we were enabled to run finer resolutions. The optimal option would be to measure experimentally this critical value and use it as reference.

The a priori under estimation of the critical value aside, the particle motions are negligible Reynolds number, shear-thickening can be explained by inertial effects [START_REF] Picano | Shear thickening in non-brownian suspensions: an excluded volume effect[END_REF][START_REF] Ness | Shear thickening regimes of dense non-brownian suspensions[END_REF] within the flow. However, standard models fail to simulate DST in the Stokes regime, suggesting another physical effect (different from inertia) must be involved in the DST. Illustration from Mari et al. [START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF] .

Among the recent researches on the sources of DST, Mari et al. [START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF] have identified the minimal set of physical ingredients needed to simulate DST. In particular, they have shown the essential role of the mechanical friction in the DST phenomena and introduced a time scale critical for a shear-rate dependent rheology.

Contrary to Brownian suspension, there is no others force scales, aside from the hydrodynamic forces, in non-Brownian suspensions in the Stokes regime. Therefore, an additional force scale must be introduced to define the second time scale of the flow independently from the inverse shear rate. The minimal model proposed by Mari et al. [START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF] , modifies the Coulomb's law as follows:

F t ≤ µ c ( F n -F CL ) for F n ≥ F CL , 0 otherwise, (4.7)
with the same notations as in Sec. 2.3.2. The shear rate dependence is given by the ratio γ/ γ0 , with γ0 ≡ F CL /6πµa 2 , a is the radius of the particle, and F CL is the critical load force (the additional force scale of the flow). Simulations using a Stokesian-Dynamics like approach have shown the critical-load model above mimics the DST phenomena (Fig. After checking that the numerical method can accurately simulate the motion of an isolated ellipsoid, lubrication corrections are tested with an experimental benchmark. The impact of a prolate particle onto a wall aims to validate the combined use of the LLCM and collision model.

Isolated particle in a sheared flow.

Before looking at lubrication correction, the numerical method described in Chap. 2 has been tested with an isolated spheroidal particle in a simple shear flow (Fig. 5.1).

The analytical solution for the angular velocity Ω z is given [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] in the inertialess regime (Re eq = 0) by

Ω z = - γ a 2 + b 2 b 2 cos 2 (Θ z ) + a 2 sin 2 (Θ z ) , (5.1)
where γ is the imposed shear rate, the semi axes a and b are the polar and the equatorial radius of the spheroid (Fig. where D eq is the particle equivalent diameter. The computational domain of size [10D eq ] 3

is mapped with an uniform Cartesian mesh. No-slip boundary conditions are assumed on the particle and wall surfaces. Figure 5.2 shows that simulation results, obtained with a grid spacing of ∆x = ∆y = ∆z = D eq /30, are in good agreement with the analytical solutions Eq. (5.1). The particle tumbling period around the spanwise axis is given by

T = 2π γ AR + 1 AR .
Figure 5.2: Evolution of the normalized rotation velocity Ω z of a prolate (red) and oblate (blue) in a shear flow. Simulation results (dashed lines), using h = 1/30 grid resolution, are compared to Jeffery et al. solutions [START_REF] Jeffery | The motion of ellipsoidal particles immersed in a viscous fluid[END_REF] . from the source into the side of the tank. The second mirror reflects the light passing through the side of the tank back toward the camera. Hence, two orthogonal views can be recorded using the single-camera available.

To better take into account the spatial distortions on the recorded frames the calibration has been made using a checkerboard pattern placed on each face of the tank and in the middle of the tank as shown Fig. 5.6. From the recording of these six positions of the Particle locations and velocities are obtained from recorded measurements in two steps.

First, each recorded frame is processed to locate the particle. Secondly, raw data from the image processing step are processed to reduce noise and to compute the particle velocities.

The position of the ellipsoid mass center placed thereon is found, on every recorded frame, throughout the following procedure (see Fig. 5.7):

1. Convert RGB image to gray-scale and mask the surrounding.

2. Thresholding: mask pixels which have a value above the threshold value.

3. Find particle outline, defined by its center position and two radii (polar and equatorial radius), using Elliptical Hough Transform [START_REF] Xie | A new efficient ellipse detection method[END_REF] . The particle location and radii found based on the particle outline are given with sub-pixel accuracy, and the radii were found to be in agreement with the measured size of the particle. The procedure described above is applied to every recorded frame of each measurement.

After the image processing of every recorded frame, the particle mass center and the particle orientation are measured in a two-dimension plane of motion (Oxy). The noise of the raw data has been reduced using the moving average method (Fig. 5.8 -a).

Transnational and rotation velocities in the plane (Oxy) are computed from the smoothed trajectory using the central difference scheme. To reduce the noise, the moving average method was also applied to the velocities (Fig. 5.8 -b).

Figure 5.8: Evolution of the vertical position (a) and velocity (b) of a prolate particle colliding with a wall. Raw data from the image processing are smoothed using moving average method.

Measurements of the particle properties.

Collision properties of the particle are directly linked to the particle and wall materials.

For the sack of simplicity collision parameters have been measured using a spherical particle. Experiments have been carried on both in air (dry collision) and water (wetcollision). During each recording, the temperature of the fluid has been monitored. The air temperature was 23 • C for all dry collisions, and the water temperature was 16 • C for wet collisions.

The normal coefficient of restitution ξ max,n of the resin particle has been measured with a 3D-printed sphere of 1cm diameter. The particle has been dropped such that the particle incoming velocity was smaller than 1m • s -1 to avoid velocity dependency of ξ max,n due to plastic deformation [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] (see Fig. 5.9). The average coefficient of restitution for five collisions is ξ max,n = 0.899 ± 0.013. The tangential coefficient of restitution ξ max,t can be assessed the same way than for the normal coefficient of restitution. In this case, the particle is impacting on an oblique wall.

The measurements need to be done at several wall inclinations to identify the inclination belonging to the stick regime (where the coefficient can be properly evaluated). These measurements have not been made due to limited time and access to the TU-Delft lab.

Therefore, only the normal component of the collision force will be considered in all simulations using ellipsoid particles and tangential deformations are assumed elastic (i.e.

ξ max,t = 1). As the mixing of glycerol and water could not be done perfectly, glycerin density and dynamic viscosity have been measured in the lab. An Ostwald viscometer (also known as U-tube viscometer) has been used to measure the kinematic viscosity of the glycerin. The dynamic viscosity is then directly calculated from the measured fluid density. Viscosity values contained in Tab. 5.1 are the average viscosities measured on three samples of both glycerin mixtures at the same temperature (which also was the glycerin temperature during the recording of the particle trajectories). Based on Tab. 5.1 and tabulated measurements [START_REF]Association. Physical Properties of Glycerine and Its Solutions[END_REF] , true glycerol mass concentrations are estimated a posteriori at ∼ 26% and ∼ 47% instead of the targeted concentrations 25% and 50%.

For each fluid, the particle trajectory has been recorded for five drops. A higher instability of the particle trajectory has been observed at large Galileo number (water and 25% glycerin). For instance, in water, two trajectories out of five are mainly contained in the camera recording plane. In the three other cases, the particle slightly rotates along the vertical axis (Ω y = 0) as shown in Fig. 5.10. In addition, a translation of the particle along the Oz direction can be observed in some cases. This could be due to the Magnus effect generated by a rotation of the particle along its main axis (Ω x = 0) initiated by the release mechanism. However, these rotations cannot be measured without markers on the particle surface. Placing markers requires to accurately measure the particle roughness to assess the paint influence on the surface, which we could not carry out. Similar behaviors are observed with 25% glycerin. However, all five recording using 50% glycerin show a stable trajectory of the particle as shown in Fig. 5.11. Therefore, these latter measurements would be mainly used to validate the LLCM as the experience has been successfully replicated several times.

Replicable experiments at high Galileo number (Ga 500) might be possible with more sophisticated release mechanisms. However, instabilities are expected above a critical Galileo number as for it is the case with spherical particle [START_REF] Breugem | Numerical study of the sedimentation of spheroidal particles[END_REF] .

Comparison with numerical simulations.

Numerical simulations have been performed on a [10D eq , 23D eq , 10D eq ] domain (Fig. The particle type (prolate) and initial position-orientation have been chosen such that Figure 5.12: Sketch of the domain used for the numerical simulations.

the particle trajectory remains in the vertical plane, orthogonal to e z (i.e. the camera recording plane). The particle is dropped close to the wall and Θ z ≈ 0. To initialize the simulations, the particle linear and rotational velocities in the plane orthogonal to z-direction (U x , U y , and Ω z ) are linearly interpolated in time from the recorded velocities.

The other velocities are not constrained but computed by the simulations. As in Sec. 4.1, The relative errors L 2 and L ∞ of the particle position X i obtained by a simulation compared to the experimental particle position X exp are denoted respectively E 2 (X i ) = X exp -X i 2 and E ∞ (X i ) = X exp -X i ∞ . Due to unsatisfactory performance of the collision model, only the dates before t ≈ 3.05s (when the particle impact the wall with LLCM) are considered to compute the following errors. With the LLCM, the errors on the positions are E 2 (X i ) = 2.7mm and E ∞ (X i ) = 4.4mm, compared to E 2 (X i ) = 3.6mm and E ∞ (X i ) = 5.7mm with the CLM. Hence, the LLCM reproduces more accurately the particle trajectory than the CLM. Furthermore, CLM overestimates 125 is larger than with 50% glycerin and the particle rebounds are clearly visible on recorded trajectories. As discuss in Sec. 5.2.1.4, the experiment has not been successfully reproduced experimentally. Therefore, we will not quantitatively compare the simulated and recorded trajectories. However, no collisions are observed in the simulated trajectory using CLM contrary to simulation using LLCM. This observation emphasizes that CLM overestimates lubrication force (as with 50% glycerin). While incoming particle velocities and orientations are identical in both simulations, the recorded orientation before impact is different to the simulated orientation (the recorded angle Θ exp z ∼ 2.7 • is wider than the simulated angle Θ simu z ∼ -0.8 • , for instance). This can explain the smaller rebound obtained with LLCM compared to the recorded ones.

Conclusion

This thesis study was focused on modeling short-range hydrodynamic interactions, also called lubrication, between rigid particles suspended in a viscous Newtonian fluid.

As the main result of this work, a local lubrication model based on the lubrication theory has been developed to correct convex particle hydrodynamics. In order to validate the local lubrication correction model (LLCM), a partitioned fluid-particle solver has been developed (Chap. 2). The fluid flow, modeled by the incompressible Navier-Stokes equations, is solved using the volume penalization method (VP), while the particle dynamics is solved by a discrete element method (DEM). The two-way coupling between the DEM and VP is made in two steps during an iteration. First, the hydrodynamic forces are computed via the interpolation of the fluid stress on the particle surfaces. Subsequently, the particle dynamics is fed back to the fluid solver by the penalization of the particle solid velocities on domain grid used to solve the Navier-Stokes equations. The LLCM is computed during the evaluation of the particle hydrodynamics.

On particles in close interaction, unresolved hydrodynamic forces are approximated by solutions found in the lubrication theory. Corrections are made in a pair-wise fashion.

Therefore, many-body interactions are not corrected but partially resolved by the solver on particle surface elements outside lubrication areas.

The LLCM has been first validated using spherical particles (Chap. 4). From the experimental measurements of Harada et al. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] , the suspension of a single particle above a wall has been successfully simulated. Comparison with an existing lubrication model [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] has shown the LLCM is as accurate as lubrication models specifically designed for spherical particles. A large range of incoming particle Stokes number has also been simulated to validate the combined corrections of the LLCM and the soft-sphere collision model [START_REF] Westerweel | Collision model for fully-resolved simulations of flows laden with finite-size particles[END_REF] made on the particle dynamics. Particle rebounds are accurately reproduced when the particle impact wall with a normal [START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] or inclined [START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] angles. However, the LLCM comes with an average 10% additional cost, reducing the model relevance for mono-disperse suspension of spherical particles compared to existing models.

Even if the LLCM is slow with spherical particles, it was initially designed for suspensions of more complex particles. Experimental study of interacting particles is more challenging for non-spherical particles than spherical ones. Consequently, experimental studies of these suspensions are at this day still scarce. To assess the accuracy of the LLCM, we carried out preliminary experimental measurements. The trajectory of a single ellipsoid impacting a wall has been measured in water and glycerin mixtures. Preliminary results have highlighted the numerical and experimental challenges to accurately record or simulate the particle path. Comparisons between simulated and recorded trajectories have shown the LLCM more accurately reproduces the particle dynamics, while CLM overestimates the lubrication forces.

The benchmark proposed for this thesis are a preliminary work to validate the LLCM.

A proper validation would require to precisely track the particle positions-orientations in the three dimensions. Measurements on this upgraded experimental set-up must also cover a large range of particle aspect ratio, impact Stokes numbers, and particle orientation at the impact. The main challenge of the suspension of non-spherical particles remains in establishing a minimal set of accurate experimental benchmarks which would be used to validate the numerical methods already available. Until that day, lubrication models can also be validated on some specific suspension where proper DNS simulations can be performed and used as references. The LLCM will be validate numerically another study following this thesis.

This work opens several perspectives.

The first obvious benefit of a lubrication model adapted to non-spherical particles is to simulate natural particle-laden flows, such as sediment transport, with a more realistic representation of the particle shapes. For instance, studies of ellipsoids in a sheared flow as shown that particles tend to orientate themselves along preferential directions [START_REF] Gillissen | Dynamics of prolate ellipsoidal particles in a turbulent channel flow[END_REF] .

Therefore, ellipsoidal suspensions could give rise to different viscosity or flow behaviors than what can be observed with spherical particles.

As long as the particle surface stays convex around the contact points, lubrication corrections from the LLCM remained accurate. Therefore, the LLCM can be a reliable model for deformable spherical particles such as droplet or bubbles, despite its higher computational cost compared to existing models.

Many similar lubrication models have been developed in the past two decades improving constantly the simulation accuracy. In hindsight, the question of how far should the compromise between accuracy and cost go becomes relevant. Highly accurate simulations could make sense in research applications where simulations are used to understand the underlying physics. However, their high costs make these simulations unrealistic for industrial applications, even with modern improvements in computational sciences. It is also not clear how costly, relatively to each other, the different methods (Chap. 1) are to simulate a given benchmark at a given accuracy. After decades of researches on numerical methods for suspensions, a critical review of the models which have already been developed would be as relevant as investing in yet another lubrication model.

  résolues par le solveur, et F deg i et T deg i les contributions complémentaires, non capturées par le solveur. Les contributions résolues sont données par l'intégration de la contrainte du fluide sur la surface de la particule. La seconde contribution est la conséquence des erreurs commises par la discrétisation du problème en espace et en temps ainsi que par la méthode numérique. Elle est donc par nature inaccessible et inévitable. L'objectif d'un modèle de lubrification est d'estimer F deg i et T deg i afin de minimiser l'influence des approximations numériques sur les forces et couples hydrodynamiques totaux (voir Sec. 2.3.1).

CHAPTER 1 .

 1 STATE OF THE ART 1.1. PARTICLE-BASED METHODS.

CHAPTER 1 .

 1 STATE OF THE ART 1.3. SUMMARY.

2. 1 . 1

 11 The governing equations.

Figure 2 . 1 :

 21 Figure 2.1: Sketch of the decomposition for the domain D into a fluid domain D f and the solid particles D s .

  j are the convective and diffusive terms, respectively. The current time step and the time step at the previous iteration are denoted ∆t n (= ∆t) and ∆t n-1 , respectively. Spatial derivatives ∂ ∂x are approximated by a second-order central finite difference scheme. The gradients computed at cell-centers and face-centers are denoted ∂ ∂x i cc and ∂ ∂x i fc .

  2.4). The weights γ w , γ s , and γ b are corresponding to linear interpolation for the west, south and back face velocity components, respectively. The same approach is considered for the opposite faces.

Figure 2 . 4 :

 24 Figure 2.4: Sketch of a grid cell with the notations used for the spatial discretization of the governing equations.

Fig. 2 .

 2 Fig. 2.6).

Figure 2 . 6 :

 26 Figure 2.6: Representation of the two parameters lub and col of the lubrication correction model. Typically, a lub ∼ 2 -3 grid cells.

P

  i and its surrounding obstacles to properly capture the fluid flow. Therefore, a local lubrication correction model is introduced to balance the degraded hydrodynamics (see Sec. 3.3). Lubrication correction on the force and torque are denoted F lub i and T lub i . The remaining of the hydrodynamics, F solv i and T solv i are obtained via the flow solver, as follows:

Figure 2 . 7 :

 27 Figure 2.7: Contact of two particles with notations associated to the soft sphere model.

Figure 2 . 8 :

 28 Figure 2.8: Sketch description of the numerical algorithm used to compute the state n + 1of the whole system from its state n. The algorithm starts at the estimation of the time step ∆t which is the time elapse between the system state n and n+1. The particle dynamics is solved in n t sub-time steps δt = ∆t nt . m denotes the current sub-state of the system between the state n and n + 1 while the particle dynamics is computed (NB: after the correction step, the state m is equivalent to the state n). Dashed boxes contain the reference to the main equations computed at the given step.

2 . 8 )

 28 , the force F solv i n+1 and torques T solv i n+1 are computed from the fields (p n+1 , u), and not (p n+1 , u n+1 ).

.35) 3 .. 4 .

 34 The forces (F coll i + F ext i + F lub i ) and torques (T coll i + T lub i ) are updated assuming the particles are located in X m+1 i , Θ m+1 i and are moving at the velocities U The particle velocities are updated as follows:

( 3 . 5 ) 3 . 1 . 1 . 1

 353111 Resolution of the Stokes stream-function.

C 2 and D 2

 2 coefficients obtained via the boundary conditions on ψ 1 and ψ 2 .

Fig. 3 . 1 .Figure 3 . 2 :

 3132 Fig. 3.1. The following steps aim to express the lubrication force as a function of the

  is the fluid vorticity, and ∇ 2 u = -∇×ζ because the fluid velocity field is divergence free. Hence, by substituting Eqs. (3.26) and (3.27) in Eq. (3.25), the lubrication force on the particle P i is expressed as a function of ψ, solution of Eq. (3.6), as follows:

Figure 3 . 3 :

 33 Figure 3.3: Distribution of the lubrication constraint around the contact point A (with a = 5 mm, U sq i = 1 m • s -1 , and µ = 1).

( 3 .

 3 41) 3.1.2.1 Resolution of the fluid velocity and pressure fields.

3. 1

 1 . THEORY CHAPTER 3. LUBRICATION CORRECTION MODELS the surfaces Z i and Z j :

(3. 48 )

 48 Hence, the first two harmonics of each quantity can now be solved using Eqs. (3.44) and (3.45) with the boundary conditions (3.47) and (3.48).

. 59 )

 59 Finally, the expression Eq. (3.42) of U , V , W and P are substituted in Eq. (3.59) to identify the coefficients of F lub i and T lub i seen as polynomial functions in :

  3.4 is ∆x = 0.06. These configurations have been used in Secs. 4.1 and 5.2.2 to compare the CLM to the local lubrication correction model (LLCM) described Sec. 3.3.

  3.1), F lub i and T lub i are not identical to the degraded hydrodynamics F deg i and T deg i . Nevertheless, simulations presented in Sec. 4.1 show that the approximation F deg i ≈ F lub i can be made as long as the gap Reynolds numbers

Figure 4 . 1 :

 41 Figure 4.1: Sketch of two cross sections of the domain with its configurations and initial location of the particle.

  2). The particle position Y h and velocity U h obtained with a grid mesh resolution h are compared to the particle position Y 1/100 and velocity U 1/100 obtained with a grid mesh resolution h = 1/100. The units of the position and velocity errors are m and m • s -1 .

Figure 4 . 3 :

 43 Figure 4.3: Total hydrodynamic force, according to the y-direction as a function of 1/during the approach phase. For each grid mesh resolution h, lubrication correction is activated for ≤ h .

Figure 4 . 4 :

 44 Figure 4.4: Total hydrodynamic force, according to the y-direction as a function of 1/during the approach phase. For all curves, the lubrication correction is activated for ≤ lub .

Figure 4 .

 4 Figure 4.5 shows that the global order of the numerical method used is first-order in L 2 and L ∞ on position and velocity. However, one of the main features of the LLCM is to rely as much as possible on the solved fluid stress to compute the hydrodynamics,

Figure 4 . 5 :

 45 Figure 4.5: Representation of the global rate of convergence on velocity (a) and position (b) of the particle center of gravity according to the norms L 2 and L ∞ . The particle position Y h and velocity U h obtained with a grid mesh resolution h are compared to the particle position Y 1/100 obtained with a grid mesh resolution h = 1/100.

4. 1 . 3

 13 Sensitivity analysis from the lubrication parameter lub .The LLCM is a model with a single parameter lub which sets the minimum gap length where lubrication corrections are needed. Several simulations have been made for different lub . Figure4.6 represents the total hydrodynamic force normalized by the particle velocity with respect to 1/ . Numerical simulations using a lub ≈ 2∆y are expected to give more realistic solutions since hydrodynamic forces better fit the solution of model H. For smaller lub , hydrodynamic forces are underestimated leading to unrealistic mechanical contact and rebound

Figure 4 . 6 :

 46 Figure 4.6: Total hydrodynamic force, according to the y-direction as a function of 1/during the approach phase for the critical lubrication distance a lub equals to 0.1∆y, 1∆y, 2∆y, 20∆y, 3∆y, or 5∆y (by increasing F/U at the highest 1/ ). Grid mesh resolution is h = 1/40 for all curves.

  4.7). When the particle is elastic (ξ max,n = 1) the energy is conserved during the collision. Otherwise, plastic deformations of the particle are mimicked by the collision model creating the observed energy dissipation.

Figure 4 .

 4 Figure 4.8 represents the distribution of the normalized effective coefficient of normal restitution of a particle versus the impact Stokes number. To be comparable to the measurements, we use the same technique and definition of the coefficient of normal restitution ξ n , detailed by Joseph et al.[START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] . By definition, ξ n = -U R /U T where U R is

Figure 4 . 7 :

 47 Figure 4.7: Evolution of the total energy (energy potential and kinetic) of the system during an elastic (ξ max,n = 1) or plastic (ξ max,n = 0.5) collision of a steel particle with a wall. In both cases the impact Stokes number is St d ≈ 6900.

Figure 4 . 8 :

 48 Figure 4.8: Distribution of normalized effective coefficient of normal restitution ξ n /ξ max,nof a single particle impacting a wall in respect of the particle Stokes number at the impact. Filled markers represent results obtained using the LLCM using the same experimental setup than Joseph et al.[START_REF] Hunt | Particle-wall collisions in a viscous fluid[END_REF] (data represented by hallow markers). The black curve represents a correlation made on experimental data proposed by Legendre et al.[START_REF] Daniel | Experimental study of a drop bouncing on a wall in a liquid[END_REF]

  called dry coefficient of normal restitution ξ max,n . Experimental measurements reported in Fig. 4.8 show that dry collision occurs at Stokes numbers, St d > 1000 where lubrication forces become negligible compared to the particle inertia. At low Stokes numbers St d 200, lubrication forces are dominant. For St d ∼ 10 lubrication effects prevent solid collision with the wall and maintain the particle suspended above the wall. Results with the local lubrication model (Fig. 4.8) are comparable to measurements and have a strong correlation with the experiment at a low Stokes number.4.1.5 Multiple rebounds.

Figure 4 .

 4 Figure 4.9 and 4.10 compare the evolution of the simulated particle vertical velocity

(4. 5 )Figure 4 . 11 :

 5411 Figure 4.11: Sketch of the simulation domain and initial location of the particle. Inside the dashed box, particle characteristic velocities and angles are displayed before and after collision with the wall.

Figure 4 .

 4 Figure 4.12 shows a comparison between the normalized incidence Ψ in = tan (ψ in )

4. 5 )

 5 is considered when the particle is immersed in water. Numerical simulations have been performed with h = 1/40 using N p = 3200 particle mesh elements. Lubrication parameters have been chosen as lub = 2 and n lub = 100.Collisions of the particle immersed in the air are called dry collisions since the interactions between the solid particle and the surrounding fluid are negligible. Simulations of dry collisions have been performed (Fig.4.12, top) to validate the collision model and to assess the collision model accuracy with neglectable lubrication effects.

Figure 4 . 12 :

 412 Figure 4.12: Comparison between the normalized incidence Ψ in and rebound Ψ out angles of steel and glass particles in aqueous solution (a) and in the air (b). Results of simulations are represented by filled dots and are compared to Joseph andHunt[START_REF] Joseph | Oblique particle-wall collisions in a liquid[END_REF] experimental measurements (by circles).

  Both particles are placed initially at the same separation distance δ x = 0.75L, but symmetrically off-axis at different values of δ y . Simulated trajectories of both particles using the LLCM and h = 1/30 grid spacing are plotted on Fig.4.14.

Figure 4 . 14 :

 414 Figure 4.14: Particle trajectories in the (Oxy) plane with initial separations δ y ∈ 0.015L, 0.035L, 0.05L, 0.105L .

Figure 4 .

 4 Figure 4.14 shows two kind of symmetrical trajectories. When δ y 2d, the particles are initially carried by the fluid flow before being repelled from each other by lubrication forces. At larger values of δ y , the particles pass each other to find new equilibria on the centerline in their respective upstream directions. The critical value δ y = 2d which

  accurately simulated. No collisions occur on the bypass and repelling trajectories as expected since St d = Re d /9 < 10 (see Sec. 4.1.4).

4. 4

 4 Shear thickening suspensions.Particles in a Newtonian fluid raise the apparent dynamic viscosity η r of the mixture and also usually give rise to a shear-rate dependant rheology. At high particle volume fractions φ, discontinuous shear-thickening (DST) can be observed (see Fig. 4.15) in inertial non-Brownian suspensions as well as at finishing Reynolds number. The shearthickening is created by the formation of clusters of particles jamming the flow. At non-

Figure 4 . 15 :

 415 Figure 4.15: Shear rate γ dependence of the relative viscosity η r , for a volume fractions 0.45 ≤ φ ≤ 0.56 and N = 1000 particles in an unbounded sheared flow.Illustration from Mari et al.[START_REF] Morris | Shear thickening, frictionless and frictional rheologies in non-Brownian suspensions[END_REF] .

4. 15 )

 15 . Similar results can also be obtained by introducing an electrostatic repulsion model instead of the critical-load model. This approach is better on a physical point of view as the DST is not modelled by the collision model which should in theory never be active in the Stokes regime. Due to limited computation resources, we were not able to simulate, with our framework, a full range of γ/ γ0 to capture the DST for a given volume fraction as Mari et al. (Fig. 4.15). However, we have noticed a dependency of the apparent viscosity to the channel width with a bounded flow.

Figure 4 . 16 :

 416 Figure 4.16: Sketch of the simulation domain.

4. 16 )

 16 . Periodic boundary conditions are considered on the lateral faces of the domain. A no-slip boundary condition is imposed at the top and bottom walls (y constant). The top and bottom wall are moving at the absolute velocity U wall in the opposite direction. The Chapter 5 Numerical simulations of ellipsoidal particles Contrarily to spherical particles, studies of the interactions between a single ellipsoid and an obstacle are scarce. In this chapter, we propose a benchmark to validate lubrication correction models for non-spherical particles.

Figure 5 . 1 :

 51 Figure 5.1: Sketch of the simulation domain and particle initial location.

Figure 5 . 3 :

 53 Figure 5.3: The initial experimental set-up.

Figure 5 . 4 :

 54 Figure 5.4: The release mechanism used during the experiments. The particle is resting on the two wires before being released (in the picture to the right).

Figure 5 . 5 :

 55 Figure 5.5: View on the top of the camera of the improved experimental set-up. Examples of direct and indirect light paths to the camera are represented by the dashed lines. For all measurements, the distance of the camera to the tank was d CT = 59cm.

Figure 5 . 6 :

 56 Figure 5.6: Example of a recorded checkerboard pattern place in the middle of the tank. The three dashed lines are orthogonal to each other.

Figure 5 . 7 :

 57 Figure 5.7: Procedure to find the outline of a spheroidal particle on a raw image (a).

Figure 5 . 9 :

 59 Figure 5.9: Evolution of the vertical velocity of a 1cm particle colliding with a wall (in air).

Figure 5 . 10 :

 510 Figure 5.10: Trajectory of the prolate drops in water. Frames are superposed at a constant 0.2s time interval.

Figure 5 . 11 :

 511 Figure 5.11: Trajectory of the prolate drops in glycerin (50%). Frames are superposed at a constant 0.2s time interval.

5. 12 ) 2 :

 122 , with D eq the equivalent diameter of the particle. Non-slip boundary conditions are assumed on the top and bottom of the domain while periodic conditions are assumed on the other faces. The domain is uniformly mesh with ∆x = ∆y = ∆z = D eq /20 grid elements. The particle surface is mapped using 12800 elements. The simulation configurations are reminded in Tab. 5.2. Particle aspect ratio AR 2 Particle dimension a × b × c 1 × 0.5 × 0Simulation configurations.

Figure 5 .

 5 Figure 5.13 shows the reproduction of the prolate trajectory in 50% glycerin using the LLCM or the CLM and compared to the experimental recording. The first 2.6s of the experimental recording are used for the initialization. Impact Stokes number is estimated from the recorded trajectory at

Figure 5 . 13 :

 513 Figure 5.13: Evolution of the prolate vertical position around the first impact with the wall. The particle is immersed in 50% glycerin mixture. Lubrication corrections are applied only on position below the horizontal green dashed line.

Figure 5 . 14 :

 514 Figure 5.14: Evolution of the prolate vertical velocity around the first impact with the wall. The particle is immersed in 50% glycerin mixture. The simulated velocities in the gray area (t 2.94s corresponding to the initialization phase) are still directly interpolated from the recorded velocities.

Figure 5 .

 5 Figure 5.15 shows the reproduction of the prolate trajectory in 25% glycerin using the LLCM or the CLM and compared to the experimental recording. Particle velocities are interpolated from recorded data until the particle enters the lubrication area (i.e. the gap between the particle and the wall is smaller than 2∆y). Impact Stokes number St Deq ∼

Figure 5 . 15 :

 515 Figure 5.15: Evolution of the prolate vertical position around the first impact with the wall. The particle is immersed in 25% glycerin mixture.

  When a particle gets close to an obstacle, the flows around the particle is modified by the presence of this obstacle. The compression and shearing of the fluid in the gap between the particle and the obstacle raise the hydrodynamic forces acting on both solid surfaces. This augmentation is caused by lubrication forces acting in the opposite direction than the relative motion of bodies. The dominant component of the lubrication forces diverges when the gap width vanishes as -1 along the fluid compression direction. Lubrication effects are essential in the suspension phenomena and challenging to solve numerically due to their singular behavior, and the time and spatial discretization intrinsic to numerical methods. Extensive researches have been made in the past five decades to understand and correct lubrication in particle-laden flows (Chap. 1). The most robust and costeffective approach introduces a model which corrects the lubrication forces in the flow using force expressions from the lubrication theory. However, classical lubrication models are rigorously limited to weakly-inertial suspension of spherical particles. Suspension of perfectly spherical particles is, at best, rare in nature. Therefore, lubrication models need to improve in order to cover a broader diversity of particle shapes. To that end, we went back to the lubrication theory and developed a lubrication model compatible with any convex particle (Chap. 3). Instead of correcting the lubrication forces on the whole particle surfaces, poorly resolved hydrodynamics acting on a neighborhood of the contact points (lubrication areas) are substituted by analytical solutions from the lubrication theory. Hence, local corrections enable the simulation of non-spherical particles and reduce the total hydrodynamics dependency from the lubrication model. Indeed, only the hydrodynamic forces in the lubrication areas are constrained to match analytical solutions found in the Stokes regime, while they are computed from a moderate Reynolds flow everywhere else. For classical lubrication models, the total hydrodynamic forces are modified such that they match the analytical solutions found in the Stokes regime when particles are close to each other even at non-vanishing Reynolds number.

  particles are near contact, the time step ∆t has to satisfy the stability condition of the collision model (see Sec. 2.3.2). Therefore, when lubrication corrections are active (i.e. collision might occur) ∆t is chosen such that τ c = N t ∆t (with τ c the contact time) and ∆t ≤ β τ ∆x Vmax , where N t > 0 an integer (N t = 8, if not explicitly stated otherwise).

  2.7. The tangential overlap distance δ t δ t δ t is obtained by integrating the relative tangential velocity at the point of contact while the Coulomb's law is verified. Therefore, the tangential overlap distance δ t F n e t -γ t (U i,j • e t ) e t ) , (2.28) when the particle is sliding (i.e. F t > µ c F n ), and by • e t ) e t dt, (2.29) when the particle is sticking to the obstacle (i.e. F t ≤ µ c F n ). The rotation tensor R δt moves δ t δ t δ t n to the new local coordinate system at the state n + 1, and

				δ t δ t	n+1 at the time step n + 1 is
	obtained by
	δ t δ t δ t (-µ c δ t n+1 = 1 k t δ t δ t n+1 = R δt • δ t δ t δ t	n +	t n+1	(U i,j
			t n

Table 4 . 1 :

 41 4.1. Simulation configurations.

	Fluid density	ρ	985	kg • m -3
	Fluid dynamic viscosity	µ	0.142	Pa • s
	Particle density	ρ p	1127	kg • m -3
	Particle diameter	d	2.54	cm
	Normal restitution	ξ max,n	0.97
	Contact time	τ c	7.98 • 10 -5	s
	Particle roughness	a col	2 • 10 -4 d	m
	Gravity field	g	9.781	N • kg -1
	Particle terminal velocity	U T	0.146	m • s -1
	Initial position	init	0.4181
	Under these configurations (Tab. 4.1), the fluid characteristic Reynolds and Stokes
	numbers are as follows:

Table 4 . 2

 42 2.13 • 10 -4 4.74 • 10 -4 1/40 3.83 • 10 -4 6.35 • 10 -4

	-	-

: Table of convergence for the particle vertical position and velocity, when the particle dynamics is fully resolved ( ≥ 0.

Table 4 .

 4 5: Particle properties. Numerical simulations using LLCM have shown that the parameter n lub , introduced in the Sec. 3.3, has a limited impact on the results when n lub > 10.

	Steel

Table 5 .

 5 1: Density and dynamic viscosity of the fluids used for the experiments. The last column contains the prolate Galileo number based on the particle equivalent diameter.

	5.2.1.4 Observations.
	Mesurements have been made with a same prolate emerged in 30cm filtered water
	(every particle larger than 25µm are removed from the water) or glycerin. Two mass
	factions of glycerol has been considered: 25% and 50% (target values). Fluid properties
	are contained in Tab. 5.1. Particle Galileo number
		Ga =	ρ p -ρ ρ	gD 3 eq ν 2 ,
	which scales the gravity force to the viscous forces acting on the particle, decreases when
	the glycerol concentration increases. Therefore, the most stable particle trajectory are
	obtained at the highest concentrated glycerin.
		Temperature ( • C) ρ (kg • m -3 ) µ (mPa • s)	Ga
	Water	21.7 ± 0.1		998	0.9601	2172
	Glycerin (25%)	24.3 ± 0.1	1060 ± 2	1.980 ± 0.059 801
	Glycerin (50%)	23.1 ± 0.1	1112 ± 2	5.096 ± 0.065 413

ILU and additive Schwarz method have been used from the PETSc library (see Sec.

2.4 for references).

NaSCar has been developed within the team Memphis at INRIA-BSO.

Portable, Extensible Toolkit for Scientific Computation (PETSc) is an open-source library of data structures and routines developed by Argonne National Laboratory for the scalable (parallel) solution of scientific applications modeled by partial differential equations.

Particle have been printed using Form

printer from Formlabs © using RS-F2-GPBK-04 black resin.

Acknowledgements

This thesis is a result of the collaboration between the university of Bordeaux and the French National Institute for computer science and applied experimental parallel testbed is developed under the Inria PlaFRIM 1 development action with support from LaBRI and

imate the degenerated forces:

Numerical experiments have shown that the number of particle surface mesh elements (starting from about a thousand elements) has a limited impact on the solutions. Therefore, all the following simulations have been performed using N p = 3200 elements for the particle surface mesh. experimental measurements made by Harada et al. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] . Numerical simulations were performed on a uniform Cartesian mesh with a grid spacing h = 1/40 and using N p = 3200 elements for the particle surface mesh. Local lubrication corrections are performed on particle mesh elements closer than 2∆y to the wall. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] experimental measurements and the model H. The LLCM and CLM are activated for smaller than the blue square and red star, respectively.

The relative errors L 2 and L ∞ of the particle velocity U obtained by a simulation compared to the experimental particle velocity U exp provided by Harada et al. [START_REF] Tanaka | Fluid force acting on a particle falling toward a wall[END_REF] are

Simulations have been performed using a h = 1/30 grid resolution and with the LLCM. Experimental measurements are provide by Gondret et al. [START_REF] Gondret | Bouncing motion of spherical particles in fluids[END_REF] Simulations results are comparable to experimental measurements. In both cases, the simulated velocities and positions before the first impact do not perfectly fit the measurements. These small differences are amplified at each rebound. However, the errors remained limited for the first rebounds. The overall accuracy of the simulation could be improved by using finer grids.

Two particles in a sheared flow.

The behavior of two particles suspended in a confined linear shear flow at finite

Reynolds number has been widely studied in the past two decades experimentally [START_REF] Ingber | Interaction of two particles in a shear flow[END_REF][START_REF] Ovarlez | The interaction of two spherical particles in simple-shear flows of yield stress fluids[END_REF] and numerically [START_REF] Morris | Hydrodynamic interaction of two particles in confined linear shear flow at finite reynolds number[END_REF][START_REF] Cardinaels | Lubrication analysis of interacting rigid cylindrical particles in confined shear flow[END_REF] . In this section, the accuracy of the particle-particle interaction is investigated. 4.6: Apparent dynamic viscosity η r for φ = 50%.. With φ ≈ 50% and H = 4d, no DST is observed, and the flow behaves like a Newtonian flow for γ/γ 0 ranging from 10 -2 to 1. This result can be a numerical artifact due to the small length of the domain and few particles (N = 72). Further testing is required. However, if the numerical simulations are accurate, the lack of DST could be explained by the particle cluster structure in a confined channel such that frictions between particles are minimized. Therefore, the critical-load is never reached, and the flow remains Newtonian.

Falling particle onto a wall.

Faced with the lack of simple experiments of an ellipsoid in interaction with an obstacle in the literature, we have made our own measurements with the support of Dr. Wim-Paul Breugem, Dr. Christian Poelma and Wout Cornel from TU-Delft. In this first test case, a single particle (prolate) is impacting a planar surface. First, the experimental procedure used to track the dynamics of a single spheroidal particle is described. Records have then used as references to assess the LLCM accuracy.

Experimental Measurements.

Experimental measurements have been made in two phases. First, the motion of the particle has been recorded in a vertical plane using a single camera. From the observations made during the first measurements, the set-up has been improved to record the particle motion in three dimensions using mirrors.

Set-up and calibrations.

Experiments were performed in a 25×25×50 cm glass tank (see Fig. 5.3). The glass is 5 mm thick. Particles are 3D-printed using stereolithography (SLA) printing technology.

Particles are made of an opaque homogeneous resin of density 1220 ± 10 kg • m -3 , and printed 1 with about 25 µm accuracy. The particle roughness h col is assumed to be in the same order of magnitude than the accuracy of the printer (h col ∼ 25 µm).

A LaVision APX-RX 1Mpx high-speed camera, mounted with a Nikon 35 mm lens, was used to record the motion of the falling particles. The particle motion was recorded at 500 frames per second with a camera exposure of 1/500 and aperture equals to F/11.

Camera calibration has been made using a checkerboard pattern.

The particle is dropped using two fine nylon wires as shown in Fig. 5.4. The wires are initially in tension and the particle is rested on top of them. The tension on the wires is slowly decreased to release the particle with a limited influence (of the release mechanism) on the fluid.

Preliminary measurements have shown that the prolate motion is mainly in the recording plane. However, the particle trajectory is sensitive to the initial conditions and in some cases, the particle slowly rotates along the vertical axis while settling (see Sec. 5.2.1.4).

Therefore, the experimental step-up has been improved to also record the particle motion in a vertical plane orthogonal to the camera recording plane (see Fig.