Prof Ulrich Nieken

Dr Manuel Hopp-Hirschler
email: manuel.hirschler@icvt.uni-stuttgart.de

Prof Dany Vandromme

Dr Thibaut Ménard

Prof. Denis Lebrun

Prof Armelle Cessou

Valérie Florence

Prof Reveillon Julien

Fatemeh Almasi
email: almasif@coria.fr

M.S Shadloo

Abdellah Hadjadj
email: abdellah.hadjadj@insa-rouen.fr

M Ozbulut

N Tofighi

M Yildiz

Ulrich Nieken
email: ulrich.nieken@icvt.uni-stuttgart.de

Mostafa Safdari Shadloo
email: msshadloo@coria.fr

Numerical simulations of multi-phase electro-hydrodynamics flows using a simple incompressible smoothed particle hydrodynamics method

Keywords: Electrohydrodynamics, Thermocapillarity, Multiphase fluid flows, SPH, MPI, HPC List . . . . . . . . . multiphase modelling. . . . . . . . . . . . . Smoothed Particle Hydrodynamics (SPH) Meshless methods Electrohydrodynamics electrohydrodynamics, thermocapillary, multiphase fluid flows

the administrative staff in CORIA Lab. Thank you for your kindness and your help. I specially thank the personnel of the library of INSA of Rouen for their efficiency and goodwill. Thank you Alexis and Cédric from CORIA IT support who responded efficiently to my questions. Many

Résumé

La quantité d'eau contenue dans le pétrole brut est à l'origine de problèmes importants dans l'industrie pétrolière. Le transport efficace du pétrole et du gaz est donc un défi majeur qui doit être relevé par le secteur de l'énergie afin de réduire le coût du transport de l'énergie et de protéger les installations de la corrosion, ce qui conduit à une réduction de l'empreinte carbone totale de l'industrie. Il est donc crucial d'envisager de nouveaux concepts pour étudier la séparation des gouttelettes d'eau en suspension dans un milieu pétrolier. Des forces externes peuvent être appliquées pour améliorer l'efficacité du processus de séparation. Les interactions entre ce système multiphasique et les champs électriques et thermiques sont donc des préoccupations croissantes. Par conséquent, le système susmentionné doit être modélisé comme un système multiphase, multiphysique et multi-échelle. L'objectif de cette thèse est d'examiner la dynamique des fluides en suspension d'un point de vue de physique fondamentale en utilisant la stratégie numérique sans maillage la plus fréquemment utilisée dans la littérature sur la dynamique des fluides numériques (CFD), appelée hydrodynamique particulaire lissée (SPH). L'interaction entre les constituants du système multiphase à l'interface fluide-fluide ainsi que les effets de l'électrohydrodynamique (EHD) et de la thermocapilarité (TC) sur le système doivent être mieux clarifiés. De plus, la nature multi-échelle des échelles spatiales et temporelles impliquées représente un défi de simulation numérique et nécessite un traitement particulier. Cette thèse a pour but d'initier une compréhension numérique des phénomènes couplés EHD-TC en analysant la dynamique des gouttelettes dans des applications multiphysiques et de présenter une validation étape par étape en augmentant les difficultés du problème. Cette thèse démontre que la méthode SPH pourrait être une approche alternative réussie en tant que pointe de la modélisation numérique pour prédire et contrôler la dynamique de problèmes multi-physiques complexes. Pour atteindre cet objectif, un modèle mathématique est proposé et une solution numérique est développée et validée en utilisant des solutions numériques, expérimentales et analytiques précédentes. Un solveur SPH parallélisé par MPI et une approche de calcul haute performance (HPC) sont sélectionnés pour traiter efficacement les aspects de calcul. Plusieurs phénomènes physiques sont simulés, notamment la montée, la déformation et la rupture des bulles. À l'interface fluide-fluide, des interactions hydrodynamiques complexes sont décrites mathématiquement à l'aide de la méthode de la force de surface continue (CSF), dans laquelle le gradient des forces de tension interfaciale et les forces de Marangoni, dues au gradient thermique, sont modélisés. La thermocapillarité génère une contrainte suffisamment importante, même pour des nombres de Marangoni faibles, qui conduit à la migration des gouttelettes. Les phénomènes EHD influencent la morphologie des gouttelettes en fonction de la conductivité électrique et de la permittivité électrique des rapports de la phase dispersée et de la phase continue. Nous avons montré que lorsque ces deux physiques affectent le système simultanément, les principaux paramètres qui caractérisent la déformation, la migration et la rupture de la phase dispersée sont les caractéristiques du champ électrique appliqué et du gradient thermique, le nombre de Marangoni et les paramètres des fluides. Il s'avère que le couplage des phénomènes EHD-TC affecte le système multiphasique et peut être réglé pour améliorer le processus de séparation des phases.

Mots-clés :

Electrohydrodynamique, Thermocapillarité, Ecoulements fluides multiphasiques, SPH, MPI, HPC hot to cold region. The flow propagates inside due to drag and reverse flow forms to hold mass conservation (a ). Droplet on a solid surface with temperature gradient leans toward the cold region due to higher surface tension at the inception of movement and Marangoni circulation forms inside it (b ). In the case of channel flows, the carrier liquid obeys a bulk flow from hot to cold region due to surface tension forces and the confined droplet moves to the opposite direction (toward the hot region) to conserve the mass (c ). Schematic reproduced form [START_REF] Karbalaei | Thermocapillarity in microfluidics-a review[END_REF] . . . . . . . . . . . . . . . The streams are modeled to describe the interaction of these phases between themselves as well as the contribution of the individual phases. Multiphase flows form a major feature of the oil extraction process as freshly extracted crude oil naturally contains water droplets or gas bubbles suspended in oil (see Fig. 1.1).

It is not easy to set a precise historical start point to the studies related to multiphase flows.

It is believed that the buoyancy discovery by Archimedes on 250 BC can be considered a commencement to the multiphase fluid dynamics development. From this moment to development of a two phase pressure drop model by Lockhart and Martinelli [START_REF] Lockhart | Proposed correlation of data for isothermal two-phase, two-component flow in pipes[END_REF] in 1949, multiphase phase flows have been an intertwined part of fluid dynamics. Shortly after and (a) Droplets migration and evaporation in boiling water [30].

(b) Electric fields application to move droplets on a surface enables more efficient experiments. [START_REF] Hardesty | Programmable droplets[END_REF].

(c) Temperature difference propels droplets [24].

(d) Water-in-oil (left) and oil-in-water (right) emulsions [2]. by increasing demand for petroleum in early 70's, many researchers started working on theoretical and experimental studies on two-phase flows. From 1990s the use of numerical models flourished that made it possible to use more complex mathematical models leading a 1.2 Multiscale approach of a multiphase simulation more detailed characterization of multiphase systems. Based on the comprehensive review by Tryggvason et al. [START_REF] Tryggvason | Computations of Multiphase Flows[END_REF][START_REF] Tryggvason | Fluid mechanics with interfaces[END_REF], despite significant progress made in the simulation of flows with variable surface tension, temperature and electric fields coupling, and flows with phase change, more work has to be done before such simulations become ubiquitous.

Coupled systems in addition to requiring the solution of a large number of equations, typically have a significantly broader variety of length and time scales. The new generation of multiphase flow problems are dealing more and more with this multi-physics aspect where incorporating 1) implicit time integration and 2) accurate meshing become crucial.

The more discontinuities a process contains, the more important it is to decide whether to go for explicit time integration with lower computational cost or expensive implicit schemes. The same dichotomy also exists when deciding whether to choose a mesh-based method requiring specific interface characterization strategy or mesh-free method where the interface can be naturally distinguished and thus complex geometries easier to control. In multi-phase multi-physics problems, both of these factors influence the computational cost and in consequence, increase the need for more power full computers and well-parallelized schemes. Among all CFD methods, Shadloo et al. [START_REF] Shadloo | Improved Incompressible Smoothed Particle Hydrodynamics method for simulating flow around bluff bodies[END_REF] has shown that SPH is capable of producing promising results in simulating multiphase-multiphysics fluid dynamics problems.

Multiscale approach of a multiphase simulation

Numerical simulations are controlled experiments that give us insight into the physical phenomena. Multiphase flows have been extensively investigated by means of numerical methods during the past few decades. Computational power has enabled simulations of unprecedented sophistication and detail, and has enabled the resolution of combined phenomena occurring at many different spatial and temporal scales. The main challenge in modeling complex systems is to determine the scale, accuracy and complexity of the model required to achieve acceptable predictive capabilities and reflect these requirements in a stable and efficient computational framework. The advantage of numerical simulations is facilitating the execution of experiments that are hard to setup in laboratory environment. In the case of droplet dynamics, numerical simulation can be easier to establish compared to experimental studies.

The main challenges in simulating droplet dynamics arises due to the resolution of multiple scales involved in the simulation as well as different physics interacting on and across appropriate scales (from the electrical forces acting at the molecular level to the hydrodynamic forces acting at the flow scale). When droplet dynamics are studies experimentally or in the industry level, the scale is the macro-scale, in designing phase separators for instance.

However, when the same problem needs to be numerically simulated, the smaller scales will influence and interfere the behavior of the system at macroscopic scales.

Fluid dynamics systems under a continuum assumption can be modeled using Navier-Stokes equations. However, in multiphase flows, subject to multiple physics, the length scale and the time scale of each physics involved (electric, hydrodynamics, thermal fields) may require additional treatment necessary to represent and capture the physics at all different scales.

Multiphysics couplings of phenomena occur on multiple temporal and spatial scales. When choosing algorithms one needs to take into account that combining existing codes through software often fail to adequately address the coupling physics as one code may violate basic physical principles assumed to hold by the other code, and sometimes algorithms suffer from issues related to disparate temporal and/or spatial scales between coupled physical processes [START_REF]Workshop I: Multiphysics, Multiscale, and Coupled Problems in Subsurface Physics[END_REF].

Objectives

Objectives

The objectives of the thesis are the following:

• The first objective of this thesis is to set up a mathematical model allowing the numerical simulation of a biphasic system using Smoothed Particle Hydrodynamics (SPH) method.

This model must take into account geometrical and physical complexities of the problem. Thereafter, we validate this model using available research data of different physics: first validation carried out in Electrohydrodynamics (EHD) and a second one in thermocapillary (TC) conditions.

• The second objective is to develop a procedure to include the temperature-dependent surface tension coefficient on thermocapillary flow simulations. Within the Siper code, developed by the a team of the ICVT and used for other multiple cases. Indeed, these simulations are relatively expensive in terms of computational resources, in particular those required for the EHD-TC coupled use case. This procedure uses the properties of the heat transfer library at the fluid-fluid interface already functional for heat transfer problems. This temperature-dependency of surface tension must confirm to already established data structure of the code and assures the accuracy of the simulations and conformation with theoretical solutions.

• The third objective is similar to the first one but it concerns a multiphase flow tool to solve the dynamics of the system under EHD-TC conditions. The modelling must allow, first of all, to ensure that the coupling of two physics happens in the conditions that were previously and separately validated. Then, it will be necessary to verify the solution for the set of parameters obtained for specific coupled system of forces.

Layout of thesis

The rest of the text is organised as follows:

Chapter 2 We present in this first section, the fundamental concepts associated with multiphase flows and in particular the problematic of surface tension. We detail in the following parts the equilibrium condition and the description of the thermocapillary droplet migration. Towards the end of this section, we provide the expression of transport equations, their characteristics and finally we breifly develope the context using the state-of-the-art publications.

Chapter 3 To further detail the important role of each force, in particular in the presence of electric field, we precisely detail the principles of the electrohydrodynmics (EHD) theory in the case of leaky dielectric fluids. The choice of this case for the simulations and the complexity of different force interactions are thereafter presented in more details. In particular, the material properties, initial and boundary conditions of a classic system are described. To take a deeper look at the application of EHD in heat transfer enhancement, we develop the last section on the existing literature and the state of the art on our specific topic.

Further on this chapter, the used numerical methods are expounded. We explicitly develop the methods used for discretisation of the equations based on the continuous form of the conservation equations in multiphase flows using the SPH method.

Furthermore, we present the description of the specific used code, the role of MPI parallelisation and HPC framework from the computational point of view. The last part of this section concerns numerical simulation of coupled EHD and thermocapillary phenomena. A mathematical model for such multiphysics systems is proposed, validated according to the previously mentioned principles, on the benchmark problem setups.

Layout of thesis

Chapter 4 The forth chapter presents the numerical results of EHD, TC and coupled EHD-TC cases. Several validation studies are provided to make sure that the numerical results are in agreement with available theoretical, experimental or numerical results in the literature. Each case has been thoroughly examined and then, analysed with supporting physical arguments.

Chapter 5 The fifth chapter includes the conclusions to this work and some perspectives to future works are proposed.

Chapter 2

Literature Review and Theoretical

Framework

This chapter is dedicated to the presentation of the one of the most commonly used grid-less numerical methods in CFD known as Smoothed Particle Hydrodynamics (SPH). We briefly present the basics of this method, its advantages and limitations, and main applications to show how SPH method has successfully been used as an alternative strategy in the numerical simulation of multi-phase multi-physics problems. To this end, we take a close look at the Electro-hydrodynamics (EHD) flows and the thermo-capillary flows.

Governing Equations

The balance equations used in the current thermo-capillary system are derived from the original model proposed by Espanol and Thieulot [39] and previously applied by Shadloo [START_REF] Shadloo | Improved Multiphase Smoothed Particle Hydrodynamics[END_REF] assuming that two immiscible and incompressible fluids are in contact with no phase change. We present the governing equations in Lagrangian framework.

The continuity equation

Dρ Dt = -ρ∇.⃗ u, (2.1)
and the momentum balance for multi-phase system subject to external electric force

ρ D⃗ u Dt = -∇p + η∇ 2 ⃗ u + ⃗ F body + ∇.Π capillary + ⃗ F e , (2.2)
and in the thermo-capillary flow case assuming no viscous dissipation and surface energy contributions, the energy balance Regarding the capillary stress tensor in equation 2.2, which is due to continuous change of density, a mathematical derivation of the sharp interface limit from a diffuse formulation is elaborated in [] from which the capillary stress tensor can be calculated as

ρc p DT Dt = ∇(κ∇T ), (2.3
∇.Π capillary = (γκ⃗ n + ∇ s γ)δ , (2.4)
1The material time derivative is a directional time derivative for a fixed point.

where γ is the surface tension coefficient, κ the curvature and ⃗ n the unit normal vector to the interface. The tangential component of the surface tension is ∇ s γ) and δ the Dirac-Delta function defined as δ = 1 at the interface 0 in the bulk. (2.5) This implies that the capillary stress tensor is constant inside the bulk. The surface tension force can be modeled using the continuum surface force (CSF) model [15].

Continuum Surface Force model (CSF)

Morris et al. [START_REF] Morris | Modeling Low Reynolds Number Incompressible Flows Using SPH[END_REF] introduced the first surface-tension model based on the Continuum Surface Force (CSF) technique proposed by Brackbill et al. [15]. This method essentially computes the local curvature of an interface. The surface tension can be defined as a force per unit volume

⃗ F s = ⃗ f s δ s , (2.6)
with δ s being the surface delta function and

⃗ f s = σ κ ⃗ n + ∇ s σ , (2.7)
where f s is the force per unit area, σ is the surface tension coefficient, κ is the local curvature of the interface, and ⃗ n is the unit normal vector to the interface. The term ∇ s σ links the surface tension to the the Marangoni effect due to change of the surface tension. The unit normal vector can be computed using the color function

⃗ n = ∇c |∇c| . (2.8)
In the presence of an external electric force, the coupled interaction between hydrodynamics and electrostatics of the fluid flow is modeled. In the following section, we discuss the foundations of the ElectroHydroDynamics theory (EHD) starting from the Maxwell equations [] as a set of coupled partial differential equations that, together with the Lorentz force law, form the foundation of classical electromagnetism.

Theoretical elements associated with ElectroHydroDynamics (EHD) theory

EHD studies on Droplet Migration

External electric fields can influence droplet dynamic characteristics such as morphology, velocity, and the manifestation of breakdown and coalescence. Identification and understanding of the underlying principles of Electrodynamics (EHD) can be applied in better controlling and predicting the motion and deformation of droplets.

Initially, fluids were divided into the perfect electric conductor and perfect dielectric. When the electric force imposed on a droplet suspended in a carrier fluid, the droplet starts to deform or displace. Since the electric charge was not introduced into the background fluid, it was assumed that the droplet migration and deformation happens as the result of dielectric and electrostrictive forces. [6,4,[START_REF] Turnbull | Electroconvective Instability with a Stabilizing Temperature Gradient. II. Experimental Results[END_REF][START_REF] Turnbull | Electroconvective Instability with a Stabilizing Temperature Gradient[END_REF] In contrast to the initial theory, Allan and Mason [3] showed experimentally that a droplet subject to the electric field deforms to the oblate shape.

From these observations, Taylor [START_REF] Taylor | Studies in Electrohydrodynamics. I. The Circulation Produced in a Drop by Electrical Field[END_REF] concluded that fluids are partially conductive and partially dielectric with finite electrical permittivity and electrical conductivity values. To explain the oblate shape of the isolate droplet subject to an electric field, Taylor proposed that the electric stress accumulated on the droplet interface has indeed both normal and tangential components. The tangential stress imbalance leads to the oblate droplet. The experiments performed by Vizikia and Torza [START_REF] Vizika | The electrohydrodynamic deformation of drops suspended in liquids in steady and oscillatory electric fields[END_REF][START_REF] Torza | Electrohydrodynamic deformation and bursts of liquid drops[END_REF] confirmed Taylor's model for droplet deformation.

Taylor established a structured basis for the EHD phenomena known as leaky dielectric theory that still continues to develop with analytical, experimental and numerical studies [START_REF] Melcher | Electrohydrodynamics: A Review of the Role of Interfacial Shear Stresses[END_REF]1,[START_REF] Worraker | The effect of temperature-induced variations in charge carrier mobility on a stationary electrohydrodynamic instability[END_REF][START_REF] Worraker | A nonlinear electrohydrodynamic stability analysis of a thermally stabilized plane layer of dielectric liquid[END_REF][START_REF] Saville | ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model[END_REF]10,38,22,[START_REF] Hase | Rhythmic motion of a droplet under a DC electric field[END_REF][START_REF] Jalaal | An experimental study on the motion, deformation and electrical charging of water drops falling in oil in the presence of high voltage d.c. electric field[END_REF][START_REF] Patel | Recent experimental advances in electrohydrodynamic conduction pumping research[END_REF][START_REF] Mandal | The effect of uniform electric field on the cross-stream migration of a drop in plane poiseuille flow[END_REF].

Numerical simulations on the leaky dielectric model commenced by Sherwood [START_REF] Sherwood | Breakup of fluid droplets in electric and magnetic fields[END_REF] who studied the deformation of a spherical droplet subject to electric or magnetic field assuming that the electrical permittivity and electrical conductivity of the disperse and continues phases are finite values. A boundary-integral scheme was used to simulate droplet equilibrium shape, breakup process as well as the steady state leaky dielectric droplet when large electric field and Stokes flow are imposed. Thereafter the Finite Element based numerical simulations of Tsukada [START_REF] Tsukada | Theoretical and experimental studies of circulations inside and outside a deformed drop under a uniform electric field[END_REF] showed good agreement with their experimental results as well as Taylor's analytical solution. To address the quantitative discrepancy between obtained results from Taylor's theory and the experimental observations, Feng [42] Hereafter, we tend to focus our attention on the deformation and migration of droplets in a carrier fluid subject to electric field. The influence of the surfactant on a viscous drop deformation is investigated based on the leaky dielectric model using the level-set method [START_REF] Teigen | Influence of surfactant on drop deformation in an electric field[END_REF]. They found that the deformation is dependent on the permittivity and conductivity ratios of two fluids as well as the direction of the electric field. The direction of the circulation zones can be clockwise or counter clockwise in a prolate deformed shape.

When counterclockwise, surfactant addition causes greater deformation at low surfactant concentration. They reported a strong relationship between the effect of the surfactant of the steady-state deformed droplet and the relative electric permittivity and conductivity of fluids. In a related study [START_REF] Tomar | Two-Phase Electrohydrodynamic Simulations Using a Volume-of-Fluid Approach[END_REF], a coupled level set and volume-of-fluid (CLSVOF) algorithm was used to simulate two-phase electrohydrodynamic flows and it was shown that a wrong choice of interpolation scheme (weighted arithmetic mean) may lead to a transition region thickness dependent electric field in the continuous phase. To control multiphase flow systems using electrohydrodynamics some studies have been performed using mehsless numerical methods. In this context, the behavior of the droplet in a uniform electric field is analysed using Weighted Least Square (WLS) method and the performance of the method is validated with theoretical results [START_REF] Yoshikawa | Meshless method based on weighted least square method for electrohydrodynamic problems[END_REF]. Another meshfree method to simulate multi-physics and multi-phase problems is Smoothed Particle Hydrodynamics (SPH) that is implemented EHD-driven phenomena such as bubble rising and droplet deformation. The results are validated through the comparison between theoretical and experimental results [START_REF] Shadloo | A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid[END_REF][START_REF] Rahmat | Numerical simulation of the electrohydrodynamic effects on bubble rising using the sph method[END_REF].

To comprehend the origins of the electric field forces, one must examine the molecular structure of the system at the microscopic scale. In the absence of electric forces, fluids are electrically neutral, containing zero net charges in the form of ions or free electrons.

Ionization process, on the other hand, occurs at the presence of the electric field within which massive amounts of charged particles are injected to the system. These charged particles are rapidly propelled in the direction dictated by the electric field provided by the electrodes.

The charged particles smash with the neutral molecules of the gas as they move between the electrodes and, on average, "push" the latter in the in the direction of the applied electric field. The Maxwell equations are composed of four equations with each one describing one phenomenon respectively.

Maxwell didn't invent all these equations, but rather he combined the four equations made by Gauss (also Coulomb), Faraday, and Ampère's law [43].

But Maxwell added one piece of information into Ampere's law (the 4th equation),

Displacement Current, which makes the equation complete.

• Gauss's law for static electric fields.

• Gauss's law for static magnetic fields.

• Faraday's law which says a changing magnetic field (changing with time) produces an electric field.

• Ampere-Maxwell's law which says a changing electric field (changing with time)

produces a magnetic field.

We first focus on the Gauss's law for static electric fields which states that the electric field produced by electric charge diverges from positive charge and converges upon negative charge. Moreover, the electric field flux passing through any closed surface is proportional to the total charge contained within that surface. The integral and differential form of this law are shown in 2.1.

Integral Form Differential Form

S ( ⃗ E •⃗ n)ds = q enclosed ε 0 ⃗ ∇ • ⃗ E = ρ ε 0
where the electric charge q produces an electric field ⃗ E. Here the n is a vector with length of one pointing in the direction perpendicular to the surface. The term ⃗ E • n represents the component of the electric field vector that is perpendicular to the surface under consideration and is used to calculate the electric flux. According to the differential form, in an static electric field, the divergence at one point equals to the electric charge volume density ρ at that point divided by the permittivity of free space ε 0 .

In EHD problems, the interactions between the electric field and the electric charges in fluids are considered. The conductivity of a of a fluid depends on the availability of mobile charge that it carries. This determines whether a substance can have a very high conductivity or be a non-conductive insulator. In transport phenomena of EHD, the electric charge distribution is not in equilibrium. The temporal change is taken into account in the Maxwell-Ampèr's law. The integral and differential form of this law are shown in 2.2.The integral form of 2.2 demonstrates that an electric current I or a changing electric flux through a surface produces a circulating magnetic field around any path that bounds that surface.

This implies that the electric currents and changes in electric fields are proportional to the magnetic fields circulating about the areas where they accumulate.

Integral Form Differential Form

C ⃗ B•d ⃗ l = µ M I enclosed + ε 0 d dt S ( ⃗ E •⃗ n)ds ⃗ ∇ × ⃗ B = µ M ⃗ J + ε 0 ∂ ⃗ E ∂t
where ⃗ B and ⃗ E are the vectors of the magnetic and electric fields, µ M is the magnetic permeability that determines a material's response to an applied magnetic field, and is the electric current density. Electric currents and the temporal change of the electric field lead to a magnetic vortex field.

Integral Form Differential Form

C ⃗ E • d ⃗ l = - d dt S ( ⃗ B •⃗ n)ds ⃗ ∇ × ⃗ E = - ∂ ⃗ B ∂t
By applying these laws in the EHD context, the influence of the dynamic charge on the magnetic induction is so small that it can be neglected. Consequently, the electro-magnetic part of the system can be described with a quasi-static electric field model. In addition, no external, time-varying magnetic field is imposed on the system, which means that the coupling between the electric and the magnetic field in Faraday's induction law is omitted [START_REF] Saville | ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model[END_REF].

If the divergence of the maxwell-ampèr law 2.2 is substituted by Gauss's law for electric fields 2.14 and observe that for a magnetic field the divergence of rotation is zero (∇.∇ × ⃗ B = 0), we obtain for the conservation of charge

Dq v Dt + ∇. ⃗ j = 0 (2.9)
Where ⃗ j is the total current density, which is proportional to the electric field vector and defined by the Ohm's law ⃗ j = σ ⃗ E. Here the proportionality factor σ is the electrical conductivity. In isotropic materials σ is a scalar, otherwise it is a tensor.

Volumetric Electric Force

In case of the application of the external electric field, a new term of electric force called as

⃗ F e will be added to the right hand side of the momentum equation 2.2. The calculation of this force reads as

⃗ F e = ∇.T e (2.10)
in which the Maxwell's stress tensor reads

T e = ⃗ D ⃗ E - 1 2 ( ⃗ D. ⃗ E)I (2.11)
where the dielectric displacement vector is defined as

⃗ D = ε ⃗ E (2.12)
with ⃗ E the external electric field and the electrical permitivity ε is the production of vacuum permitivity and relative permitivity :

ε = ε 0 ε r (2.13)
Basically as mentioned before, in this study the electrostatics and hydrodynamics of the system are coupled together. This coupling is through the Maxwell stress tensor. The stress induced in an incompressible liquid medium due to the presence of an electric field can be described as 2.11. That is because in 2.11 the electrodynamic part of the system can be regarded as quasi-static model, dynamic currents values are so low and the induced magnetic field effects are negligible, therefore the contribution from the induced magnetic field is neglected. Here, ⃗ E is the electric intensity, and ⃗ D is the electric displacement vector such that 2.12. It is also noteworthy that in this two phase system both fluids are considered to be leaky dielectric fluids, i.e. electric relaxation time is much shorter than viscous counterpart

(t E << t µ ).
As for the Gauss Law, the divergence of the electric displacement field could be obtained from the free electric charge density q v , that is

∇. ⃗ D = q v (2.14)
Application of 2.11 and 2.14 in 2.10 will result in

⃗ F e = q v ⃗ E - 1 2 ⃗ E. ⃗ E∇ε. (2.15)
The electric force in equation 2.15 consists of two terms. The first term from the left is the electric field force. Due to the interaction of the free charges with the electric field, the electric field force acts along the electric field lines. Thereby the free charges influence the direction of the electric field force. The electric field force thus acts between two charged particles and corresponds to the Coulomb force in the classical sense. For a system with two perfectly dielectric fluids, this term is omitted because in this case there is no free charge [START_REF] Griffiths | Introduction to electrodynamics[END_REF].

The second term is the polarization force. It results from the pairs of charges acting along the normal direction toward the interface and originating from the gradient of permittivity. The polarization force acts between uncharged particles or an uncharged particle and a charged particle. This is possible because charge separation occurs within the uncharged particles, inducing an electric dipole. The equation 2.9 can be further simplified. A homogeneous fluid with constant fluid with constant permittivity and electrical conductivity is assumed.

If the Gaussian law for electric fields 2.14 is substituted into the equation of the current density ⃗ j = σ ⃗ E and the result is then substituted into the equation of charge conservation 2.9, the following yields in both integral and differential form ( see table 2.4).

Here, t e = ε σ is the relaxation time required to return from the deflected state to the equilibrium state. For EHD problems, the time t can be described by the viscous time of

Integral Form

Differential Form

q v = q v 0 exp -t t e qv + q v σ ε = 0 fluid motion t µ = ρL 2 µ .
Here, L is the characteristic length. Depending on the conductivity level of each component, a two-phase system can be classified as

• Dielectric-dielectric (t e >> t µ ), • Dielectric-conductive • Conductive-conductive (t e << t µ ).
The distinction is made by means of the electrical conductivity and permittivity of the fluids involved. Equation (2.18) can be used to describe the retardation of free charge in a fluid can be described.

If a fluid is only weakly conductive, a perfect dielectric material can be assumed

(t e >> t µ ).
If the electrical conductivity is large, then the charge equalization at the surface is almost instantaneous compared to the time scale of the Fluid motion (t e << t µ ).

This means that the conservation of charge from Eq. 2.9 in such a system reaches equilibrium much faster than the fluid can react. Consequently, it is a conducting fluid.

In this work both fluids are considered to be leaky dielectric. In the following section we concisely present this model and its conditions.

Leaky Dielectric Model

In the mid-1960s, Taylor [START_REF] Taylor | Electrically Driven Jets[END_REF] developed a model to describe fluid motions under the influence of an electric field. Leaky dielectric model is based on the following assumptions:

• An important modeling assumption in the multiphase flows is considering that the dispersed phase and the continuous phase weekly conductive fluids forces from the bulk free charge are negligible thereby the coupling between electric field and the hydrodynamic system takes place at the boundaries.

• Fluids must have finite electric conductivity in a quasi-static electric field

• electric relaxation time is much shorter compared to its viscous counterpart (t e << t µ ) Since the fluid motion is much slower than the charge balance, an equilibrium state can be assumed for equation 2.9 ( Dq v Dt = 0). Such a system of two fluids can be classified as conducting-conducting. The conservation of charge 2.9 yields [START_REF] Saville | ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model[END_REF] ∇.(σ ⃗ E = 0).

(2.16)

Having that the rotation of a gradient from a scalar field results in zero 2, and the electric field is rotation-free3, the electric field vector can be expressed by the electric potential φ

⃗ E = -∇φ . (2.17)
Substituting the equation 2.17 in the equation 2.16 results in ∇.(σ ∇φ ) = 0.

(2.18)

2∇ × ∇φ = 0 3∇ × ⃗ E = 0
Using dimensional analysis performed by Saville [START_REF] Saville | ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model[END_REF], we can simplify the equation 2.18 if the t E << t µ as following

σ ⃗ E .n = 0. (2.19)
This equation shows that the electric charge on the interface is continues. This implies that the electric potential is also continues at the interface ( φ = 0). For a fluids with constant electrical conductivity, the equation 2.19 is equivalent of ∇ 2 φ = 0 which returns electric potential. Hence, the electric field can be calculated from equation 2.17.

Looking back to the objective of this section, which is to calculate the electric field force from equation 2.15, one needs to calculate q v using electric field ⃗ E. The volumentric charge density q v , calculated from equations 2.14 and 2.12, reads as

q v = ε∇. ⃗ E + ⃗ E∇ε (2.20)
A more detailed elaboration on the derivation of EHD system of equations can be found in [108,[START_REF] Shadloo | Improved Multiphase Smoothed Particle Hydrodynamics[END_REF]26].

Numerical Simulation of Thermocapillary-driven Flows

Since surface-driven flows have been found to be an important phenomenon in many industrial processes, and in order to expand the SPH application for complex multiphysics phenomenon, Hopp-Hirschler et al. [START_REF] Hopp-Hirschler | A smoothed particle hydrodynamics approach for thermo-capillary flows[END_REF] presented an ISPH model of a thermo-capillary flow driven by a surface tension gradient. For accurate results, this model applies CSF approach including Marangoni forces and is using (i) density-invariant divergence-free (DIDF) [START_REF] Hu | An Incompressible Multi-Phase SPH Method[END_REF] (ii) corrected SPH [14,[START_REF] Shepard | A Two-Dimensional Interpolation Function for Irregularly-Spaced Data[END_REF] and (iii) particle shifting (PS) approaches for multi-phase systems. Their validation started from single-phase flows where they studied the effects of corrected ISPH, DIDF and PS approaches on particle distribution using a Taylor-Green vortex. It was demonstrated that only PS enables a homogeneous particle distribution. In the same study, lid-driven cavity was studied for different Reynolds numbers. The results were compared to reference solutions taken from OpenFOAM. A convergence study showed that by increasing the grid resolution, the SPH model converges to the OpenFOAM results.

Next, the authors validated diffusive mass transport, compared it to analytic solution and found very accurate agreement. Then, to validate the coupling of momentum and energy equation, they investigated non-isotherm single-phase flow.

In a case study of buoyancy-driven cavity, where the fluid is accelerated due to a temperature difference between two walls, it was shown that the fluid density decreases with increasing temperature. Since the gravity is not neglected, a circular flow begins to accelerate the denser fluid downwards and the lighter fluid upwards (see figure 2.1). After an extensive validation and convergence study by comparing SPH and OpenFOAM results, Hopp-Hirschler et al. [START_REF] Hopp-Hirschler | A smoothed particle hydrodynamics approach for thermo-capillary flows[END_REF] investigated multi-phase systems. In thermo-capillary flow, a droplet migration happens as a result of gradients of the surface tension tangential to the interface, caused by a gradient of the temperature, even in the absence of gravity. For sufficiently low Reynolds numbers, the shape of the droplet remains spherical during this motion. The setup of this test case and the convergence study compared to Tong and Browne [START_REF] Tong | An incompressible multi-phase smoothed particle hydrodynamics (sph) method for modelling thermocapillary flow[END_REF] and Ma and Bothe [START_REF] Ma | Direct numerical simulation of thermocapillary flow based on the Volume of Fluid method[END_REF] are shown in figure 2.2.

Fig. 2.2 (left) Schematic setup of thermocapillary migration of a droplet and (right) Droplet migration using different resolutions compared to solution from [START_REF] Tong | An incompressible multi-phase smoothed particle hydrodynamics (sph) method for modelling thermocapillary flow[END_REF] and [START_REF] Ma | Direct numerical simulation of thermocapillary flow based on the Volume of Fluid method[END_REF]. The box length to radius ratio is L/R=4. reproduced from [START_REF] Hopp-Hirschler | A smoothed particle hydrodynamics approach for thermo-capillary flows[END_REF].

In thermo-capillary-driven flows, a non-isotropic pressure contribution, due to the interface, is expected in the momentum balance. This contribution is referred to as capillary stress in general or surface tension in an immiscible system [START_REF] Kirkwood | The Statistical Mechanical Theory of Surface Tension[END_REF]. The capillary stress tensor in the momentum balance covers both normal and a tangential stresses. Therefore, one may break-down the capillary stress into a normal and tangential part. The tangential part, known as the Marangoni force, consists of the gradients of the surface energy along the interface.

First, we studied normal part of the surface tension in detail using static and dynamic test cases. This component of the surface tension is validated using the Young-Laplace law in a two dimensional steady-state problem. The analytic law was captured very well. Deviations of the pressure at the interface are a result of the CSF model and where previously found in other models as well. The error was

The notion of the thermocapillary motion of a flat liquid layer and the motion of a spherical non-deformed droplet was published in the Russian language by Fedosov in 1956 and translated to English in 2013 [41] who solved analytically two thermocapillary cases one of which being of importance for our study. After investigating the motion of a drop suspended in a viscous medium due to temperature gradient, he states that because of the temperature difference at different points of the drop surface, the later cannot remain at rest. Instead, it will move from the warmer regions with the lower surface tension toward the colder regions with the higher surface tension. By doing so, the moving drop will drag the surrounding liquid medium by applying some force to it. Based on Newton's law of action and reaction, there will also be an equal force in opposite direction which the surrounding medium will apply onto the drop. This reactive force will push the drop in the direction of the temperature gradient. Later on, a more recognized description of the thermocapillary migration of a liquid droplet and a gas bobble was studied by Young and colleagues in 1959 [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF] who derived the analytical solution for the terminal velocity of a single suspended bubble subject to temperature gradient as well as the temperature gradient required for the bubble to hold at stationary state. Similar to Fedosov [41], their solution was subject to several assumptions including an always spherical droplet, unbounded domain of carrier fluid, and an imposed Stokes flow conditions. Before discussing the relaxation of these constraints, one must get informed about dimesionless numberrs that could govern the system. The term thermocapillary flow is often used to describe the Marangoni flows induced by surface tension gradient caused from temperature gradient [33]. In case of a suspended droplet in the background fluid, the gradient of the surface tension on the interface, exerts a shear stress on the continues phase and generates a dynamic from the hot to the cold region. In reaction to this movement, the droplet migrates towards the hot region (see (c).Fig 2 .3). This is due to the fact that in most fluids as the temperature increases, the surface tension decreases to obtain a more favorable total energy condition.

The gradients of temperature can be decomposed into two basic modes.

• The first mode refers to the temperature gradients perpendicular to the interface between phases that provoke the production of Bénard-Marangoni circulation zones [START_REF] Pearson | On convection cells induced by surface tension[END_REF] and the emergence of Marangoni instability phenomenon. The flow propagates inside due to drag and reverse flow forms to hold mass conservation (a ). Droplet on a solid surface with temperature gradient leans toward the cold region due to higher surface tension at the inception of movement and Marangoni circulation forms inside it (b ). In the case of channel flows, the carrier liquid obeys a bulk flow from hot to cold region due to surface tension forces and the confined droplet moves to the opposite direction (toward the hot region) to conserve the mass (c ). Schematic reproduced form [START_REF] Karbalaei | Thermocapillarity in microfluidics-a review[END_REF] • The second mode is generated from temperature gradients tangential to the interface through which a capillary flow from low surface tension to high surface tension, with the system aiming to attain minimum total energy, will propagate.

Essentially, surface tension gradients turn the static fluid flow system to a dynamic condition where the hydrodynamic stress at the interface is balanced by the second mode. The aforementioned constraints are gradually relaxed in further studies. For example, the effect of convective heat transport on the thermocapillary droplet migration was studied by Subramian [START_REF] Subramanian | Slow migration of a gas bubble in a thermal gradient[END_REF] for, unlike the initial work of Young [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF], low Reynolds and Marangoni numbers by keeping the inertial terms in the momentum equation. The theoretical and experimental studies got extended to find the migration velocity of a undeformable droplet in low Reynolds numbers [START_REF] Subramanian | Thermocapillary migration of bubbles and droplets[END_REF][START_REF] Meyyappan | The thermocapillary motion of two bubbles oriented arbitrarily relative to a thermal gradient[END_REF][START_REF] Subramanian | The Stokes force on a droplet in an unbounded fluid medium due to capillary effects[END_REF]. Based on the theoretical studies by Balasubramanian and Subramanian [8] at the high Marangoni and high Reynolds number regime, the fluid droplet velocity is proportional to the Marangoni number, whereas for gas bubbles the migration velocity stays independent of the Marangoni number.

Assuming an immiscible two-phase Newtonian, viscous, incompressible, isothermal fluid flow, the corresponding mass conservation and linear momentum conservation in a

Lagrangian formulation for the balk will be written as [START_REF] Shadloo | Improved Multiphase Smoothed Particle Hydrodynamics[END_REF]:

Continuity (mass conservation) : Dρ Dt = -ρ∇.⃗ v (2.21)
Linear momentum conservation :

ρ D⃗ v Dt = ∇.T + ρ ⃗ f b + ⃗ f v (2.22)
where ⃗ f b is the Body force ( gravitational force) and ⃗ f v is the volumetric surface tension force which is the summation of the Surface tension force ⃗ F S and the Electric field force ⃗ F E .

Herein, ρ is the density of each fluid, D Dt is the material time derivative operator4, ∇.⃗ v is the divergence of the velocity vector, T is the total stress tensor which is defined as T = -pI + τ when p is the absolute pressure, I is the identity matrix and τ is the viscous term that is equal to µ(∇⃗ v + (∇⃗ v) T ) for µ being the dynamic viscosity and T the transpose operator. where T and T are the total stress tensor exerted from one fluid to the other at the interface and σ characterizes the surface tension. Similarly, the hydrodynamic tangential stress component must also be balanced at the interface that yields for every t unit tangential vector to the fluid's surface to

n.T.t -n. T .t = ∇σ .t (2.24)
In order to solve this problem, a system of equations should couple the classical Navier-Stoks equations and the heat equation:

∂ T ∂t + u.∇T = κ∇ 2 T (2.25)
Thereafter, the notion of σ (T ) shows the linear dependency of the stress σ to the temperature T that affects the boundary conditions. A concrete example is the effect of a thermal gradient dT dx on the vertical motion of the bubble suspended in a more viscous fluid. The theoretical solution lies in the assumption that the bubble remains spherical and therefore, the normal stress of the interface is completely balanced by the normal component of the surface tension. Considering other boundary conditions and the kinematic condition, Young and Goldstein [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF], proposed the solution for the coupled system of Navier-Stocks and heat transfer. It is found that in the low Re regime, a new force directed from the cold fluid to the warm fluid and the magnitude of 2πR 2 dσ dx will act on the bubble, that is the origin of the thermocapillary instability.

The idea of the capillary motion of fluids due to surface tension gradient was first mentioned in 1855 by Thomson [21] where he develops several experimental scenarios to explain the capillary attraction phenomena at the surface of alcoholic liquids. However, it took more than a decade until a more concrete theoretical explanation for the capillary forces due to surface tension changes presented by Marangoni [START_REF] Marangoni | Ueber die Ausbreitung der Tropfen einer Flüssigkeit auf der Oberfläche einer anderen[END_REF] where he measures the capillary constant of some fluids and develops the surface tension forces in different droplets and continuum phases are examined. In most materials there is a linear relationship between the surface tension and the temperature gradient. The critical temperature where the linear relationship starts to vanish, corresponds to the critical Marangoni number [START_REF] Vyas | An investigation of marangonibenard convection in water based nanofluids[END_REF].

Bénard-Marangoni convection [13], also known as the Thermocapillary (TC) effect, refers to the mass transfer across the fluid-fluid interface when temperature gradient leads to surface tension gradient. The resulting unbalance shear stress generates fluid motion as shown in Fig. 2.4

Previously mentioned investigations assume that the fluid carrying the gas bubble or liquid Literature Review and Theoretical Framework droplet is infinitely extended while in experimental and numerical studies this condition is not applicable and the object is confined in a finite ambient fluid or by adjacent solid. To address this issue, the quasi-steady problem of thermocapillary migration of a bubble normal to free fluid surface with constant temperature or infinite planar solid were investigated [START_REF] Meyyappan | The thermocapillary motion of two bubbles oriented arbitrarily relative to a thermal gradient[END_REF][START_REF] Sadhal | Flow past a liquid drop with a large non-uniform radial velocity[END_REF]. In another study [START_REF] Meyyappan | Thermocapillary migration of a gas bubble in an arbitrary direction with respect to a plane surface[END_REF] the gas bubble was placed parallel to the rigid planar surface.

Contrary to what was previously found [START_REF] Young | The motion of bubbles in a vertical temperature gradient[END_REF], the velocity was smaller than the predicted value. The analytical, experimental and numerical research on this phenomenon has been extensively reviewed [START_REF] Ishii | Thermo-Fluid Dynamics Theory of Two-Phase Flow[END_REF][START_REF] Subramanian | Motion of bubbles and drops in reduced gravity[END_REF]32,[START_REF] Karbalaei | Thermocapillarity in microfluidics-a review[END_REF][START_REF] Schatz | EXPERIMENTS ON THERMOCAPILLARY INSTABILITIES[END_REF].

The droplet can be thermally actuated using two main micro manipulation principles. It might be in direct contact with the heat source [START_REF] Hadland | Thermocapillary migration of bubbles and drops at moderate to large marangoni number and moderate reynolds number in reduced gravity[END_REF]. It can also be actuated using a non-contact mechanism such as electric and magnetic fields [START_REF] Legchenkova | Magnetically inspired deformation of the liquid/vapor interface drives soap bubbles[END_REF], and the optical tweezers that use a laser beam [START_REF] De Saint | Thermocapillary migration in small-scale temperature gradients: Application to optofluidic drop dispensing[END_REF].

In the next section, we elaborate the numerical studies on multiphase systems.

Numerical Studies on Multiphase Flows

Introduction

Despite extensive research and advancements made in the numerical modelling of multiphase flows, such as development of computing power, this topic remains a significant challenge.

In this thesis, the phenomena comprising the mixing of fluids as well as the chemical reactions between the components are excluded from the scope of this thesis. Focusing particularly on two-phase flows, the primary phase which occupies the bulk (higher volume fraction) is called the "continuous phase" and the other phase that exists within the main phase (lower volume fraction) is referred to as the "dispersed phase".

In principle, a multiphase flow system can experience three classes of forces: (i) Line, such as surface tension forces; (ii) Surface, such as viscous forces and pressure forces; and (iii) Volume, such as inertia and buoyancy. The complete scope of total resulting forces on a multiphase system should take into account both external forces (such as the electric or magnetic forces), and the interactions between the forces acting on each phase, as well as the traction force of one phase to the other(s). The complexity of the domain of influence of each force raises challenges in the quantitative analysis of these systems considering the length scales in many orders of magnitudes from sub-micron to centimeters or more in a unique problem. A typical example could be the gas bubble or liquid droplet rising with a radius r in the order of millimeters while the continuous phase could be in the order of centimeters. Considering the effects of externally imposed fields, it is known that only under the static conditions, the electric and magnetic fields influence the flow independently. In such a situation and if the characteristic time scale of the electrostatic process is larger than that of the magnetic phenomena, the latter could be neglected, and thus the electrostatic equations provide an acceptable approximation [START_REF] Saville | ELECTROHYDRODYNAMICS: The Taylor-Melcher Leaky Dielectric Model[END_REF]. In multiphase multiphysics problems are highly non-linear which necessitate dealing with stable/unstable conditions and transient behavior of the bubble/droplet during its deformation and migration. In the same context, another challenging point encountered while modeling turbulent multiphase systems is to define and find apposite eddy length and time scales in different zones of the problem (i.e. core flow, shear layer, and recirculation zones) [START_REF] Kataoka | Basic equations of turbulence in gas-liquid two-phase flow[END_REF].

Lack of solid and comprehensive theories and simulations for predicting and calculating force interactions in multi-physics multiphase flow models need to be properly addressed to obtain a superior interpretation of the experimental observations and pave the way for technological advancements.

A Multiphase Approach

Given that in continuum fluid mechanics, in the length scales of a few micrometers to a few meters, fluids are considered as a single condensed matter, neglecting the intermolecular interactions. An instinctive property of multiphase fluid flows is that despite assuming macro-scale forces, it is equally crucial to confront molecular forces. This dichotomy of approaches comes from the fact that at the fluid-fluid interface, different forces are exerted from each phase to the other, that is to say, the interfacial molecular interactions need to be considered. A classical two-phase system of water and air is shown in Fig. 2.5. In such a situation, a molecule in the fluid balk experiences similar forces from adjacent molecules in every direction since they are all in the same phase (isotropic distribution). Whereas, a water molecule at the interface is only subject to downward traction forces (See Fig. 2

.5).

As a consequence of the action of different forces compared to that of the bulk molecules and hence having access to less potential energy, all molecules at the interface are trapped in a discrimination condition which is defined as the interfacial energy. To attain a more favorable condition, fluid particles tend to minimize surface areas. That is why a water droplet tends to form a spherical shape in its equilibrium condition. The interfacial energy per unit of area is called the surface tension γ. Equivalently, one can imagine that surface tension is a negative surface pressure and therefore can be modeled as the force exerted at the interface per unit length. Surface tension originates at the molecular level, however, its representative parameter in multiphase configurations γ is defined at the macroscopic scale.

CFD methods in multiphase fluids

Conventional numerical methods, also referred to as mesh-based methods, use either front-tracking or front-capturing methods to treat the interface. In the front-capturing procedure, the interface is determined by a three-dimensional distribution of fluid properties. A well-known example is the volume of fluid (VOF), one of the front capturing techniques known for its efficiency in simulating free boundaries embedded in a mesh with Eulerian or Arbitrary Lagrangian-Eulerian cells. Nevertheless, the major drawback of VOF method is that it is computationally costly to capture complex geometries and the results can be more prone to inaccuracy due to mesh regeneration. The front tracking methods on the other hand, Mess-less methods are another class of numerical methods used to simulate multiphase flows. Here focus on the Smoothed Particle Hydrodynamics (SPH) methods.

The standard SPH schemes discretise the mass and momentum conservation which result in diffused density discontinuities at the interface, in particular in high density ratio problems.

To address this issue, considerable advancement of the SPH in multiphase flow simulation was initiated by Colagrossi and Landrini [27] research. They proposed the first multiphase SPH model that does not include density gradient terms , hence the treatment of systems with high density ratios at the interface became possible. Three years later, another approach to facilitate the treatment of density discontinuities was proposed by Hu and Adams [START_REF] Hu | An Incompressible Multi-Phase SPH Method[END_REF]. In this interpolation-based formation, the density of a particle of interest i can be computed only based on its position and the kernel radii. The mass conservation in this scheme is guaranteed characterized in that the temporal integration of the mass conservation equation is not needed. Table 2.5 summarizes the main numerical approaches in multiphase modelling. interfacial flows and bubbly flows.

Models are not independent of numerical parameters and higher computational work compared to Eulerian approach.

Mesh-free methods

In engineering problems there are several meshfree methods used to solve partial differential equations (PEDs). The following list summarizes some of the main meshless methods along with a major work on each method.

• Smoothed Particle Hydrodynamics (SPH) [START_REF] Gingold | Smoothed Particle Hydrodynamics -Theory and Application to Non-Spherical Stars[END_REF].

• Element Free Galerkin Method (EFGM) [12].

• Meshless Local Petrov-Galerkin method (MLPG) [5].

• Reproducing Kernel Particle Method (RKPM) [START_REF] Liu | Reproducing kernel particle methods for structural dynamics[END_REF].

• Point Interpolation Method (PIM) [START_REF] Wang | A point interpolation meshless method based on radial basis functions[END_REF].

• Generalized Finite Element Method (GFEM) [7] SPH is considered to be one of the most suitable meshless methods to simulate multiphase flows, in particular, in the case of two non-miscible fluids interacting with each other, i.e oil and water mixture or air bubbles suspended in water. Being a Lagrangian method, SPH is capable to capture non-diffusive interface. This is a major advantage compared to Eulerian mesh-based methods where one needs to establish specific interface tracking strategies which necessitate re meshing and thus, higher computational cost. In a numerical study to quantify slamming loads experienced by the subsea oil and gas installations, Sasson et al.

[109] compared the SPH and Reynolds Averaged Navier-Stokes (RANS) Volume of Fluid (VOF) method and found that despite the fact that RANS results correlate slightly better with the experimental data, the robustness and quick set up of the SPH simulations makes it as promising as RANS method for this application.

Foundation of SPH

Smoothed Particle Hydrodynamics (SPH) is one of the oldest mesh-less Lagrangian methods which discretizes the domain into a set of nodes (material particles) within which there is not necessarily a pre-defined relationship and hence, particles can freely move on any arbitrarily computational domain. Initially introduced by Gingold and Monaghan and Lucy [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] for astrophysics applications, SPH was soon found to be suitable for fluid dynamics problems where complex geometries, discontinuities or large deformations are involved.

The method was used to study phenomena such as the formation and evolution of stars and galaxies. The collective motion of the bodies is very similar to that of a liquid or a gas and can therefore be modeled with conventional equations of fluid mechanics. By its independence of mesh i.e. removing the necessity to mesh generation and refinement, SPH offers notable efficiency in calculating partial derivatives used in transport equations. With a wide range of applications from movies' special effects to formation of stars [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF], SPH has attracted many researchers significantly in recent years. Zainali et al. [START_REF] Zainali | Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method[END_REF] have investigated the Newtonian and non-Newtonian two-phase flows 5 using ISPH method where they used a cubic spline for the continuum surface force (CSF) model. This approach models the sharp interface between two fluids with a transition region of a finite thickness. Their results

show significantly improvement in the quality of the calculated interface, elimination of the inter-phase particle penetrations, and calculation of more accurate velocity and pressure fields ( see Figure 2.6).

Multi-phase systems could be classified based on the nature of the fluid such as gas, liquid, or bubbly flows [START_REF] Laurila | Numerical study of bubbly flow in a swirl atomizer[END_REF]. In this approach of classification, finite number of phases spread through the volume of continuous phases such as droplets, drops, particles or bubbles.

Two examples of systems with this structure are shown in Figure 2.7.

5The droplet is the dispersed phase and the surrounding media is the continuous phase.

Fig. 2.6

Comparison of ISPH and VOF method [START_REF] Zainali | Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method[END_REF].

Depending on the percentage of volume fraction of the dispersed phase, a hierarchy of coupling forms between the particle(s) and the background media is illustrated in Figure 2.8

Multiphase flows can also be categorized in different flow regimes by taking into account the principal dimensionless number of the system, such as low Reynolds number suspension flows and high Reynolds number turbulent dusty gas flows, among others.

The manifestation scope of the natural phenomena and industrial applications of multiphase flows reflects in a wide spectrum commencing from the scales above the molecular level in biofluids in the human body with applications such as the electrical impedance measurements to full-scale two-phase oil-water separation process in the petroleum industry.

Not only are the multiphase flows naturally present in many environmental phenomena as in stratified and free-surface flows, but they also govern the large majority of processing technologies and contribute to the development of new multi-functional materials. To address different aspects of multiphase flows, several books and review papers have been published. In a comprehensive review, Shadloo et al. [START_REF] Shadloo | Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges[END_REF] developed industrial applications of modeling multiphase flows and associated transport phenomena using SPH method due to its capabilities of handling complex boundary evolution as well as modeling complicated physics. Two-phase flow in porous media [START_REF] Sadeghi | Three-dimensional lattice Boltzmann simulations of high density ratio two-phase flows in porous media[END_REF], Two-Phase Nanofluid flows [37], multiphase turbulent flows [START_REF] Safdari Shadloo | Special topic on turbulent and multiphase flows[END_REF]118] and heat transfer in multiphase flows [START_REF] Keepaiboon | Two-phase flow boiling in a microfluidic channel at high mass flux[END_REF][START_REF] Haghighi | A review on heat transfer characteristics of cryogenic heat pipes[END_REF][START_REF] Shahsavar Goldanlou | Numerical Investigation on Forced Hybrid Nanofluid Flow and Heat Transfer Inside a Three-Dimensional Annulus Equipped with Hot and Cold Rods: Using Symmetry Simulation[END_REF]. .

Smoothed Particle Hydrodynamics Method

Mathematical foundation of SPH

Smoothed Particle Hydrodynamics (SPH) is one of the oldest mesh-less Lagrangian methods which discretizes the domain into a set of nodes (material particles) within which there is not necessarily a pre-defined relationship and hence, particles can freely move on any arbitrarily computational domain. Initially introduced by Gingold and Monaghan and Lucy [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF][START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] for astrophysics applications, SPH was soon found to be suitable for fluid dynamics problems where complex geometries, discontinuities or large deformations are involved. The method was used to study phenomena such as the formation and evolution of stars and galaxies. The collective motion of the bodies is very similar to that of a liquid or a gas and can therefore be modeled with conventional equations of fluid mechanics. By its independence of mesh i.e.

removing the necessity to mesh generation and refinement, SPH offers notable efficiency in calculating partial derivatives used in transport equations. With a wide range of applications from movies' special effects to formation of stars [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF], SPH has attracted many researchers significantly in recent years.

The comment made by Monaghan [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] is a clear and concise starting phrase to introduce two main concepts in SPH that are smoothing function and smoothing radius, as he states:

At the heart of SPH is an interpolation method which allows any function to be expressed in terms of its values at a set of disordered points-the particles.

The term "Any function" could be interpreted as any field variable related to the mechanical properties of the fluid such as mass, momentum, velocity, or pressure from a physical point of view.

Moreover, To determine the appropriate "set of disordered points" a radius of neighbourhood, sometimes called as smoothing radius, needs to be defined inside which the particles have a significant influence on each others hydrodynamic properties.

There are many techniques to find the suitable smoothing radius which is neither too large to be computationally costly, nor too small to exclude important information and result in discontinuities.

Kernel Function

As stated before, family of Lagrangian particle based methods for simulating the nonaxisymetric problems of astrophysics, known as Smoothed Particle Hydrodynamics (SPH), was introduced separately by Lucy [START_REF] Lucy | A numerical approach to the testing of the fission hypothesis[END_REF] and Gingold and Monaghan [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF]. The studies of [START_REF] Shadloo | Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges[END_REF][START_REF] Liu | Smoothed particle hydrodynamics (sph): an overview and recent developments[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] provide a general overview of the SPH approach and its numerous applications.

Here, we give a quick overview of the SPH approaches as well as the pertinent discretizations employed in the current work. This method provides an easier computation process while keeping the reasonable accuracy threshold. Being a mesh-less method, SPH allows the calculation of partial derivatives used in momentum and mass conservation to be done much faster using the concept of interpolating kernel function. The SPH method is based on the fact that any field variable A(⃗ x) (such as fluid density, flow velocity, etc ) which is a function of spatial coordinates, can be represented by its Dirac-Delta distribution. This exact equation could be mathematically formulated as:

A(⃗ x) = Ω A( ⃗ x ′ )δ (|⃗ x -⃗ x ′ |)d ⃗ x ′ , (2.26
A(⃗ x) = Ω A( ⃗ x ′ )W (h, |⃗ x -⃗ x ′ |)d ⃗ x ′ , (2.27)
Having the initial particle spacing L 0 in a Cartesian grid, h values can be chosen using

2L 0 < h < 3L 0 .
To ensure the convergence of the process, the kernel function must meet the following conditions:

• Normalization condition:

Ω W (⃗ x)d⃗ x = 1, (2.28)
• For h going to zero, W tends toward Dirac function

lim h→0 W (h, |⃗ x -⃗ x ′ |) = δ (|⃗ x -⃗ x ′ |), (2.29)
• W is differentiable Moreover, there is also important to mention the compact support condition of the kernel function as below:

W (h, |⃗ x -⃗ x ′ |) = 0 x -x ′ > kh, (2.30)
The first kernel function used by Gingold and Monaghan [START_REF] Gingold | Smoothed particle hydrodynamics: theory and application to non-spherical stars[END_REF] reads as

W (h, |⃗ x -⃗ x ′ | = exp -(⃗ x -⃗ x ′ ) 2 /h 2 /(h √ π). (2.31)
The kernel function can also be considered as a weighting average where the weight function decreases with the increase in distance. When this distance tends to zero, the weight function tends to Dirac delta function. Depending on the distance of the neighboring particles, these exert a stronger or stronger influence on the considered particle.Since from a certain distance this influence goes to zero, only particles within a certain radius are included.These particles carry the index j. With the help of the kernel function, the influence of particle j on particle i is weighted. In this work the following kernel function is mostly used according to Wendland [START_REF] Wendland | Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree[END_REF] as a fourth-order kernel function.

W h, ⃗ x -⃗ x ′ = f w h d      1 -q 2 4 (1 + 2q) for 0 ≤ q ≤ 2 0 for 2 < q (2.32)
Here, d is the dimension and f w is a normalization constant with f w = 3 4 , 7 4π , 21 16π for d = 1, d = 2, and d = 3 dimensions, respectively.

q = |⃗ x-⃗ x ′ |
h is the non-dimensional smoothing length.

The Wendland kernels are computationally more convenient than the higher-order Bsplines, allowing large number of neighbour particles and hence better numerical convergence (note that computational costs rise sub-linear with number of neighbor particles). One could suggest an analogy between the Shape (approximation) function in Finite Element Method, the probability density function in probability theory, and the interpolating kernel function in SPH. They all are used to approximate the value of a function by being multiplied then integrated over the entire space.

This will take us to the first golden rule of SPH by Monaghan that in order to find the best kernel for physical interpretation of an SPH equation, one should assume the kernel function as Gaussian. Initially, the proposed kernel functions where in such a way that each particle had to have interaction with all other particles [START_REF] Gingold | Smoothed Particle Hydrodynamics -Theory and Application to Non-Spherical Stars[END_REF]. By introducing the concept of Neighbouring particles6, kernel functions on a compact support based on a cubic B-spline functions were substituted. Numerically, only a small number of particles in the entire space affect the approximated value of the kernel function. The radius after which the value of W 6particles that are located within the range of the kernel function with respect to the particle of interest. Outside this range, the kernel function has already dropped to zero. falls off rapidly is called the smoothing radius (h). The first-order derivative yields

∂W (h, |⃗ x -⃗ x ′ |) ∂ |⃗ x -⃗ x ′ | = C      1 -q 2 3 (-5q) for q ≤ 2 0 for q > 2 (2.33)
where C = 4 7πh 2 . Figure 2.9 illustrates the kernel function W and its derivative with smoothing radius h = 0.75. Based on the definition of the kernel function and the equation (2.27), we can approximate the value of each field variable by substituting the integration operation with summation operation over finite number of interpolation points

A(⃗ x) a = ∑ b m b ρ b A (⃗ x b )W (h, |⃗ x a -⃗ x b |) (2.34)
Therefore, an approximated value for density as a function of position is

ρ(⃗ x) a = ∑ b m b W (h, |⃗ x a -⃗ x b |) (2.35)
where the indexes a and b represent the particle of interest and the particles in the neighbouring distance, called neighbor particles. m b and ρ b are the mass and density of particle b, respectively.

One of the fascinating features about the SPH method which distinguishes this method from all other mesh-based methods is that in SPH instead of directly calculating the gradient of the field variable f (⃗ x), one can use the differentiation of kernel function [START_REF] Gingold | Smoothed Particle Hydrodynamics -Theory and Application to Non-Spherical Stars[END_REF].

In particular, the formulation of the first derivative of pressure using Monaghan [START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF] reads

ρ a ∇P a = ∑ b m b (P b -P a )∇ a W ab (2.36)
where the kernel function is abbreviated as

W (h, |⃗ x a -⃗ x b |) = W ab .
The equation 2.36 is the pressure gradient estimation between to neighbour particles (a)

and (b).

In SPH modelling, evaluating the pressure at each point depends on the fluid properties, in particular its degree of compressibility. In the literature two main approaches are proposed for each type of fluid.

• Weakly Compressible SPH (WCSPH) where an artificial compressibility is used to evaluate pressure by solving an equation of state Shao and Lo [START_REF] Shao | Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface[END_REF].

• Incompressible SPH (ISPH) where the Poisson equation is solved to impose fluid's incompressibility.

WCSPH has been mostly used to simulate wave propagation Chakraborty and Balachandran [23] and flow through porous media Tartakovsky et al. [START_REF] Tartakovsky | Smoothed Particle Hydrodynamics and Its Applications for Multiphase Flow and Reactive Transport in Porous Media[END_REF]. In WCSPH, the incompressiblity in enforced using a stiff equation of state to adjust pressure oscillations by density variations. The main disadvantage of WCSPH comes from the fact that in this approach the speed of sound is used to satisfy the stability (CFL) condition. High sound speed leads to much smaller time steps making this method computationally costly.

Interactions in EHD-TC Phenomena

At the last paragraph of his work on thermocapillary motion of droplets, Fedosov [41] and after obtaining the validity limit for his analytical solution that is temperature change along the drop perimeter to be less than 50 [K], he mentions that the thermocapillary drop velocity can overcome (and sometimes, for some cases does so significantly) the electrocapillary drop velocity for weakly conductive drops.

In recent decades, the combined interactions of electric and thermal fields applied to multiphase fluid systems have received strong attention. From the theoretical physics point of view, a comprehensive description of the essence of the emergence of electrothermo-hydro-dynamics instabilities has motivated many theoretical and experimental studies [4,[START_REF] Takashima | Elektrohydrodynamic instability in a dielectric fluid between two coaxial cylinders[END_REF][START_REF] Sitte | Influence of the dielectrophoretic force on thermal convection[END_REF][START_REF] Shooshtari | Experimental and numerical analysis of electrohydrodynamic enhancement of heat transfer in air laminar channel flow[END_REF]. One the other hand, several industrial processes have applied combined electric and thermal fields. For instance, heat exchange devices [29,16], solar energy [START_REF] Ghalamchi | Optimizing of solar chimney performance using electrohydrodynamic system based on array geometry[END_REF] and microelectronic devices [START_REF] Yazdani | Heat transfer enhancement of backstep flow by means of ehd conduction pumping[END_REF] have been the motivation behind many studies.

Moreover, the electrohydrodynamics and thermal analysis principles can also be applied in zero to low gravity situations where buoyancy-driven forces are negligible [START_REF] Fylladitakis | Review on the history, research, and applications of electrohydrodynamics[END_REF]31].

The literature of heat transfer enhancement using active techniques [4] such as imposing electric field to the system has been extensively reviewed [START_REF] Laohalertdecha | A review of electrohydrodynamic enhancement of heat transfer[END_REF][START_REF] Léal | An overview of heat transfer enhancement methods and new perspectives: Focus on active methods using electroactive materials[END_REF]. In a numerical study by Hassen [START_REF] Hassen | Electro-thermo-capillary-convection in a square layer of dielectric liquid subjected to a strong unipolar injection[END_REF], the effect of the electric field in controlling thermocapillary instabilities is analysed. They used stream function-vorticity formalism to solve the coupled system of equations including Navier-Stokes, Electro-hydrodynamic (EHD), and heat transfer. It is concluded that depending on the direction and strength of the electric field, thermocapillaty instabilities can increase or decrease.

Another phenomenon where in thin liquid films temperature-induced Marangoni flows and electric field can generate flow instabilities leading to pattern formation is studied by

Corbett [28]. The objective is to use EHD and thermocapillary effect in modifying surface topography by creating smaller structures. Their linear analysis of the lubrication equation on perfect dielectric and leaky dielectric films shows that in the former case, thermal effects tend to overcome the effects of the electric field, while in the latter case, that is the leaky dielectric theory, simultaneous application of the thermal and electric fields can increase the growth rate and reduce the dominant wavelength of the instability. Nonlinear simulations that are performed using a fourth-order accurate finite-difference scheme confirm the results of linear stability analysis and reveal that the coupled thermocapillary and electric forces serve to generate smaller wavelengths than when only one of these sources of instability is used.Nazaripoor [START_REF] Nazaripoor | Thermo-Electrohydrodynamic Patterning in Nanofilms[END_REF] extended previous 1D study to 2D based on lubrication theory. The linear analysis agrees quantitatively with [28] and the numerical simulations support a decreasing trend in the feature size. The mechanism of pattern formation and the shape and size of the structures are investigated by varying Marangoni number and thermal conductivity ratio of layers. Increase in the Marangoni number is found to produce smaller structures and maintaining the same mechanism of pattern formation. The thermal conductivity ratio on the other hand, modifies the mechanism where higher conductivity ratios lead to higher thermocapillary forces in thicker zones that itself prevents the formation of uniform and hexagonal ordered small structures called pillars. The Fig. 2.10 illustrates the process of pillars formation under the effects of thermal and electrical fields.

As mentioned earlier, gradients of surface tension affects the dynamics of a liquid droplet in the presence of the electric field. In a theoretical analysis, Chang and Berg [25] found that surface tension gradient caused by introducing surfactant may produce quasi-stagnant zones along the interface either near the poles or the equator of the droplet. The Marangoni flow may also increase or decrease the droplet terminal velocity compared to when only electric field is applied. Modeling the multiphysics problem of droplet migration and deformation driven by TC and EHD flows has been a challenge for a couple of reasons. First, the physical properties of the system (i) hydrodynamic features such as dynamic viscosities, densities, interfacial tension (ii) electric features such as conductivity and permittivity ratios, and electric field configuration (iii) and thermal features such as thermal conductivities, specific heats, and the equation of surface tension as a function of temperature [START_REF] Scott | An electrically driven multistage countercurrent solvent extraction device: the emulsion-phase contactor[END_REF][START_REF] Mohan | Marangoni effects under electric fields[END_REF][START_REF] Kuipa | Influence of High Voltage Electric Fields Applied across a Horizontal Liquid-Liquid Interface on the Rate of Metal Extraction Using a Rotating Diffusion Cell[END_REF]. Therefore, before introducing the physical model following assumptions about the system need to be made.

• Continuous and disperse phase are incompressible Newtonian fluids;

• Surface tension between the droplet and the bulk is temperature-dependant, while other physical parameters remain independent and constant in time.

• Thermal field is uniformly distributed and is obtained through a constant temperature gradient applied at upper and lower solid boundaries.

• DC electric field with the electric strength lower than the critical value is applied on the upper and lower solid boundaries.

• Electric double-layer theory [START_REF] Mugele | Electrowetting: from basics to applications[END_REF], as well as magnetic effects, are neglected.

A numerical investigation by Traoré [START_REF] Traoré | Heat transfer enhancement induced by electrically generated convection in a plane layer of dielectric liquid[END_REF] on the effects of electric field and thermal gradient in parallel plates immersed in a single-phase dielectric fluid shows that when a strong unipolar injection is applied to the top or bottom boundary, the Nusselt number

Nu = L o ∂ θ ∂ y y=o
is significantly increased resulting in heat transfer enhancement. [START_REF] Lafaurie | Modelling Merging and Fragmentation in Multiphase Flows with SURFER[END_REF].

Conservation of energy

To compute the effect of temperature gradient (∇T ), the conservation of energy should also be considered for each fluid [START_REF] Nazaripoor | Thermo-Electrohydrodynamic Patterning in Nanofilms[END_REF].

ρc p DT Dt = ∇(λ ∇T ), (2.37)
and

κ = λ ρc p , (2.38)
resulting in:

ρc p ( δ T δt + v.∇T ) = κ(∇ 2 T ). (2.39)
where κ is the thermal diffusivity, λ is the thermal conductivity, and c p is the specific heat capacity at constant pressure.

In a recent study by Wakif et al. [START_REF] Wakif | A semi-analytical analysis of electrothermo-hydrodynamic stability in dielectric nanofluids using buongiorno's mathematical model together with more realistic boundary conditions[END_REF], a semi-analytic instability analysis induced by electro-thermo-hydrodynamic (ETHD) effects in dielectric nanofluids is performed by means of power series method where a simultaneous action of AC electric field and imposed negative temperature gradient are considered. The fluid density as a function of temperature is

ρ = ρ p φ * + ρ o (1 -φ * )[1 -β (T * -T 2 )]. (2.40)
Where T * the temperature, T 2 the cold wall temperature, φ * volumetric function of nano particles and β is the thermal expansion coefficient. In addition to density, the effective dielectric constant ε * is also temperature-dependent as .

ε * = ε o [1 -e(T * -T 2 )]. ( 2
These dimensionless parameters destabilize the nanofluidic system. That is to say, in case of a great ETHD stability, to assure the onset of electroconvection and enhanced heat transfer, higher velocity is required.

EHD Applied to Heat Transfer Enhancement

Coupling the electric and thermal fields introduces new variations in the the heat transfer coefficient based on the heat flux, the electrode voltage and the electric field features.

Investigating the mechanism of these processes is crucial due of the extensive use of heat exchangers to ensure large heat fluxes in industrial applications.

The literature of heat transfer enhancement using active techniques [4] such as imposing electric field to the system has been extensively reviewed [START_REF] Laohalertdecha | A review of electrohydrodynamic enhancement of heat transfer[END_REF][START_REF] Léal | An overview of heat transfer enhancement methods and new perspectives: Focus on active methods using electroactive materials[END_REF]. The potential of EHD forces in improving the heat transfer performance of multi-phase systems has been extensively demonstrated. It is therefore necessary to find fluids that are both thermally and electrohydrodynamically efficient. Starting from single phase flows, [START_REF] Herman | Selection of working fluids for electrohydrodynamic experiments in terrestrial conditions and microgravity[END_REF] state that the fluid should have conductive behaviour to facilitate the fluid movement. In two-phase flows, in addtion to the conductivity of the fluids, the permitivity is also an important parameter to to modify an interface or to create an electroconvection phenomenon. In practice, Hydrofluorocarbons such as R123 and R134-a were suggested first, yet due to atmospheric damages they are progressively substituted by non-ozone depleting chemicals such as Hydrofluoroethers such as HFE-7000 and HFE-7100 [START_REF] Nassar | Study of the dielectric behavior of hfe-7000 in function of electric field and temperature variations[END_REF]. Having an immense environmental impact, this area demands more attention in future researches.

Consernig the numerical simulations on the EHD effects on heat transfer, the Lattice Boltzmann numerical method is used to simulate the transition process from laminar to chaotic flow in ETHD convection of a dielectric fluid [START_REF] Li | Transition to chaos in electro-thermoconvection of a dielectric liquid in a square cavity[END_REF]. In a two-dimensional square cavity, the temperature field is uniformly imposed along with a DC electric field on the lateral walls. To keep the injected charge density q o constant, a unipolar, autonomous, and homogeneous injection is used in the model. To describe the charge transport the Nerst-Planck equation is used as follows

∂ q ∂t + ∇.(qKE + qu -D∇q) = 0. (2.42)
Where u = [u x , u y ] the fluid velocity vector, E = [E x , E y ], K the ion mobility and D is the charge diffusion coefficient. The simulations are executed for a range of electric Raleigh number T = ε∆φ µK and based on the results three transition states for the system are identified (i) The quasi-periodic sequence involving four incommensurable frequencies (ii) the intermittency sequence, and (iii) the alternating periodic-chaotic sequence.

In addition to single-phase flows, EHD has been used for decades in heat transfer enhancement in multiphase fluid flows [START_REF] Jones | Electrohydrodynamically enhanced heat transfer in liquids -A review[END_REF]. In particular, heat transfer in melting [START_REF] Luo | Efficient lattice Boltzmann method for electrohydrodynamic solid-liquid phase change[END_REF], boiling [START_REF] Siedel | Electric fields effect on the rise of single bubbles during boiling[END_REF][START_REF] Zu | A numerical investigation of electrohydrodynamic (ehd) effects on bubble deformation under pseudo-nucleate boiling conditions[END_REF][START_REF] Quan | Effects of electric field on microbubble growth in a microchannel under pulse heating[END_REF][START_REF] Leu | Three-dimensional electrohydrodynamic enhanced water evaporation using needle-arrayed electrodes[END_REF][START_REF] Zhang | Effects of pin and wire electrodes on flow boiling heat transfer enhancement in a vertical minichannel heat sink[END_REF] and condensation [17,36] process can be enhanced by means of the Coulomb force associated with the liquid-extraction phenomenon and through enhancement of wettability over the heat transfer surface [18,[START_REF] Zaghdoudi | Electric Field Effects on Pool Boiling[END_REF]. An unfavorable argument against using EHD techniques in such phase change mechanisms is the energy required to feed high voltage DC electric field [START_REF] Mangrulkar | Recent advancement in heat transfer and fluid flow characteristics in cross flow heat exchangers[END_REF].On the other hand, lower pressure drop and higher performance of the EHD-driven methods remain concrete advantages in comparison with passive heat transfer enhancement methods(i.e. chemical composition modification). Table 2.5 summarises some of the main thermal processes in which the EHD is used to enhance heat transfer. In the context of the coupling mechanism when the temperature gradient is imposed on an EHD fluid system, according to the pioneering study by Castellanos [20,19], the temperature equation can be obtained from the entropy equation in moving dielectric fluids. It is assumed that the entropy of the system can be expressed as a function of temperature, mass density, and displacement vector. In polar fluids, the flux of internal energy noted as J' = -k∆T

where k is a function of the temperature, the electrical permittivity, and the electric field.

The other link between the thermal and electric fields lies in the reversible variations of temperature with the electric field. For the case of EHD, the molecular polarization will be in equilibrium with the field, thus the dielectric heating can be neglected. Similarly, when fluids are highly insulating the contribution of Joule heating is not present. The author concludes that imposing temperature difference affects the dynamics of the system through buoyancy and electric forces.

In a recent numerical simulation using the Spectral Element Method (SEM), both uniform and sinusoidal temperature boundary conditions are applied on an enclosure cavity [START_REF] Wang | Spectral element method for numerical simulation of ethd enhanced heat transfer in an enclosure with uniform and sinusoidal temperature boundary conditions[END_REF].

As illustrated in Fig. 2.11, local heat transfer rate enhancement is obtained with greater electric Rayleigh number T = ε∆φ µK where ε the electric permittivity, φ the electric field potential, µ the dynamic viscosity, and K is the ionic mobility.

In an extensive theoretical analysis, Douiebe et al. [35] studied the coupled buoyancy and thermocapillary instability, known as Bénard-Marangoni problem, in an incompressible fluid with small electrical conductivity. Having considered the effects of both the AC electric field and the rotation, they determined the regions where the oscillatory instabilities occur. The system consists of a Boussinesquian fluid phase placed in a horizontal infinite layer for which the density ρ and the dielectric constant ε are varying linearly with temperature and can be expressed as

ρ = ρ o [1 -α 1 (T -T o )] (α 1 > 0), (2.43
)

ε = εo[1 -α(T -T o )].
(2.44)

The upper free surface deforms and has a low heat conductivity leading to linearly decreasing surface tension as a function of temperature.

γ = γ o -τ(T -T o )] (τ > 0). (2.45)
In the same model [35], in addition to the electric force f e , the Coriolis force f c and the centrifugal force f g are also included in the momentum balance [START_REF] Kaddame | Bénard-marangoni convection in a rotating fluid with and without surface deformation[END_REF][START_REF] Vidal | The influence of coriolis force on surface-tensiondriven convection[END_REF].

f c = 2ρV × Ω, (2.46
)

f g = 1 2 ρ∇(|Ω × r| 2 ).
(2.47)

They found that the free surface deformation affects the coupling between electric and capillary forces by means of the ratio as well as the strength of the AC electric field. It is also worth noting that the Biot number (representing the heat flux flow through the interface) h = dh T λ improves the coupling between the buoyancy, thermocapillary and electric forces. In an experimental study performed by [34], the effect of electric field augmenting the heat transfer in the evaporation/boiling process is investigated. For the microchannel evaporator, when they apply a flux of 28.91 W/cm2, the heat exchange coefficient is improved by a factor of 1.45 by imposing an electric field of 1.2 kV/mm. In another configuration with the R141b as the working fluid, under a flux of 167.8 W/cm2, by applying a field of1.3 kV/mm, the heat exchange coefficient is improved by a factor of 1.5. These observations demonstrate that for low fluxes, heat transfer is dominated by vaporization. Under the effect of the electric field, charges accumulate at the interface and create a shear force at the interface, creating an electro-convection movement that promotes heat transfer.In such mechanism, the interface is subjected to the dielectric force which extends the meniscus in the evaporator and thus increases the transfer coefficient.In the case of large heat flows, the transfer is dominated by boiling. Bubbles will stretch along the field lines, increasing the size of the liquid-vapour interface and promoting their departure from the wall.

Many approaches and developments can be explored in this area, from a modeling point of view as well as from an analysis point of view. Moreover, to delineate the governing mechanisms on coupled EHD-TC problems, the role of electrochemistry and its interaction with EHD principles requires more analytical and computational studies [START_REF] Mori | From electrodiffusion theory to the electrohydrodynamics of leaky dielectrics through the weak electrolyte limit[END_REF]11,[START_REF] Schnitzer | The taylor-melcher leaky dielectric model as a macroscale electrokinetic description[END_REF].

Summary

In this chapter we reviewed the existing researches performed on multiphase fluid flows subject to elctrohydrodynamics (EHD), thermocapillary, and the simultaneous effects of these phenomena in coupled systems. It mainly concerns the theoretical, experimental and numerical aspect of the electro-thermo convection as well as the coupled electric field and heat transfer. The following remarks may be taken into account in the future researches

• A comprehensive study on the influence of temperature-dependent variables such as the temperature-dependent model of the electrical permitivity and electrical conductivity when included could enhance the model's reliability and accuracy.

• Since in multiphase flows the simultaneous effects of the hydrodynamic, thermal and electric fields on the stability of the system depends on several dimensionless number, more works on the influence of non-dimensional numbers on the stability of the system is required;

• According to both experimental and theoretical studies,the forces may lead to a reduction in the thermal boundary layer, increased heat convection, gas bubble or liquid droplet migration and deformation, causing interfacial instabilities or phase migration that can result in a flow regime redistribution. Direct Numerical Simulation (DNS) and Large Eddy Simulation (LES) approaches in the turbulent flow problems are therefore helpful in gaining more insight into these phenomena;

• Same as previous point, numerical simulation of some classical instabilities such as Kelvin-Helmholtz and Rose-window instabilities in leaky dielectric fluids could be envisioned.

• Development of 3D numerical simulations capable of dealing with multi-scale problems in spatio-temporal aspects is also an improvement axis in this topic.

Chapter 3

Modelling and simulation of multiphase flows

In this chapter the governing equations of a multiphase system and the numerical methods used to solve them are summarized. The new functions to model temperature-dependent surface tension for coupled EHD-TC simulations are explained.

Numerical methods for multiphase flow

The numerical modeling of multi-phase, here in particular two-phase, flows remains a considerable challenge despite extensive research and the progresses in computer power. Within Computational Fluid Dynamics (CFD), the Smoothed Particle Hydrodynamics (SPH) methods is most used due to its particle-based structure and easy implementation without the need for a mesh generation [START_REF] Shadloo | Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges[END_REF]. [44].

Governing equations

In order to describe the isothermal movement of fluids, the conservation equilibrium for mass and momentum for an infinitesimally small fluid element need to be considered.

The conservation equation for energy is omitted for systems with constant temperature.

Assuming isotropic and incompressible fluids, the Navier-Stokes equations read as surface forces that in 2-D can be solved by the following equations:

∇ • v = 0 (3.1) ρ( ∂ v ∂t + v.∇v) = -∇.p + ∇.Π viscous + ρg + σ κδ s n (3.2) δ ρ δt + ∇.(ρv) = 0 (3.
∂ u ∂t = - ∂ u 2 ∂ x - ∂ (uv) ∂ y + µ ρ ( ∂ 2 u ∂ x 2 + ∂ 2 u ∂ x 2 ) - 1 ρ ∂ ρ ∂ x + f x , (3.4)
∂ v ∂t = - ∂ v 2 ∂ x - ∂ (uv) ∂ x + µ ρ ( ∂ 2 v ∂ x 2 + ∂ 2 v ∂ y 2 ) - 1 ρ ∂ ρ ∂ y + f y (3.5) 0 = ∂ u ∂ x + ∂ u ∂ y (3.6)
It is a system of partial differential equations, which is described as three unknown which contain the velocity in x-or y-direction and the pressure. To solve this, initial conditions and In case of non-isothermal multiphase flows, the conservation of energy yields

ρc p ∂ T ∂t = ∇(λ ∇T ), (3.7)
and

κ = λ ρc p (3.8) results in ∂ T ∂t = κ(∇ 2 T ) (3.9)
where κ is the thermal diffusivity, λ is the thermal conductivity and c p is the specific heat capacity at constant pressure.

wherein the heat capacity c p and the thermal conductivity coefficient λ . In equation 3.2, the viscous dissipation can be neglected since the characteristic flow velocity is much smaller than the speed of sound.

Surface tension

While the equations 3.1,3.2,and 3.2 describe the conservation laws inside the balk flow, to model the jump across the fluid-fluid interface Σ for any quantity

φ = lim ε→0 + [φ (x 0 + εn Σ ) -φ (x 0 -εn Σ )] (3.10)
where the surface normal n Σ is the normal vector at the interface from one phase to the other.

To consider the no-slip condition at the interface, that equivalent to the continuity of the tangential components of the velocity at the interface, we will have

ρ(v -v Σ ) • n Σ = 0, (3.11)
A similar condition must be assumed for the stress continuity at the interface, based on the momentum conservation which reads as

ρv ⊗ (v -v Σ ) + pI -S • v Σ = f Σ , (3.12)
where f Σ is the surface force density and is defined as

f Σ = σ κn Σ + ∇ Σ σ (3.13)
Where the effect of the intrinsic surface viscosity is neglected and the sum of the principal curvatures κ can be defined as

κ = ∇ • (-n Σ ). (3.14)
In the thermocapillary flows the surface tension depends linearly on the temperature, also known as the Eötvös rule, named after the Hungarian physicist Loránd Eötvös, and can be expressed as

σ (T ) = σ (T 0 ) + σ T (T -T 0 ), (3.15)
Where the temperature coefficient is a material dependent property and is mostly negative, i.e. the surface tension mitigates with the increase in the temperature. It is also necessary to include the continuity of the temperature at the interface as

T = 0. (3.16)
In multiphase problems, the sharp interface limit of the capillary stress tensor [START_REF] Hopp-Hirschler | A smoothed particle hydrodynamics approach for thermo-capillary flows[END_REF] can be mathematically modeled as:

∇.Π capillary = (σ κ n + ∇ S σ )δ . (3.17)
where σ is the surface tension coefficient, κ is the curvature, ∇ S is the gradient of the tangential component of surface tension, and n is the normal vector to the interface. The Dirac-Delta function δ is defined to be unity at the interface and zero elsewhere. This implies that the capillary stress tensor is constant inside the balk. The surface tension force can be modeled using the continuum surface force (CSF) model [15].

Considering the total stress tensor balance at the interface in multiphase system, figure 3.1 depicts the normal and tangential components of the total stress tensor, as well the unit vector.

Laplace's law

The Laplace's law which was discovered by the French mathematician Piere-Simon Laplace states that the finite thickness of a sphere and its radius (initial state of a theoretical bubble) determines the tension on the wall when the sphere is filled with a fluid at a specific pressure. This reveals that to reduce the tension on the wall, one needs to increase the thickness of the sphere wall at the same pressure.This law is used in results' chapter to verify the accuracy of the proposed model for pressure jump calculation. Equation 3.19 states the analytical statement of the law of Laplace [108] where p is the pressure at the two sides of the interface. Also, σ and r are the surface tension coefficient and radius of the bubble, respectively.

p in -p out = σ r (3.19)
It is noteworthy to discuss the reason behind the different values obtained for theoretical and numerical pressure jump at the interface. As mentioned before, the pressure and other field variables (velocity, etc) are approximated using the numerical smoothing scheme which converts the sharp values at the interface to smoother ones resulting in a loss of accuracy and introduction of spurious (parasite) oscillations near the wall of the bubble. These currents are generated because of an inappropriate evaluation of curvature of the circular bubble due to unreliable values for unit normal vector (⃗ n = ∇C |∇C| ) in the surface tension force calculation.

Balance laws at the interface

The conservation laws for the interface will be as follows:

Mass conservation :

[ρ(⃗ v -⃗ u)].⃗ n = 0 (3.20)
Linear momentum conservation :

[ρ⃗ v(⃗ v -⃗ u) -T].⃗ n∇ (s) σ + κσ⃗ n (3.21)
where [] sign represents the jump for the quantity f defined on x on the interface Γ such that :

[ f ](x,t) = f + (x,t) -f -(x,t) (3.22) f ± (x,t) = lim h→0 f (x ± hn(x,t),t) (3.23)
assuming the negative and positive sign of the moving discontinuity. Also in 3.21 ⃗ n is the unit normal to the interface, ∇ (s) is defined as the surface gradient operator, σ is the surface tension and κ is the curvature. Moreover in 3.20 and 3.21, the ⃗ u is the velocity of the surface Γ(t) at a point on the surface which may be decomposed into a normal component, u, and a tangential component, ⃗ u tan .⃗ n = 0

We will then have :

⃗ u = u⃗ n + ⃗ u tan (3.24)
3.4 Surface-tension force in SPH

Introduction

Surface tension has been modeled and implemented in SPH by Tartakovsky et al. [START_REF] Tartakovsky | Smoothed Particle Hydrodynamics and Its Applications for Multiphase Flow and Reactive Transport in Porous Media[END_REF].

In this approach the micro-scale interactions between particles are modeled using van der Waals equation of state. This approach is less used in engineering applications as it is dependent to the micro-scale parameters. To address this issue, macro-scale models of surface tension based on Continuum Surface Force are widely used to compute the local curvature of an interface. This section is dedicated to the formulation of this model and the computation of the normal vector to the interface.

Continuum Surface Force model (CSF)

Morris et al. [START_REF] Morris | Modeling Low Reynolds Number Incompressible Flows Using SPH[END_REF] introduced the first surface-tension model based on the Continuum Surface Force (CSF) technique proposed by Brackbill et al. [15]. This method essentially computes the local curvature of an interface. The surface tension can be defined as a force per unit volume

F s = f s δ s (3.25)
where δ s is a surface delta function and

f s = σ κ n + ∇ s σ (3.26)
where f s is the force per unit area, σ is the surface tension coefficient, κ is the local curvature of the interface, and n is the unit normal vector to the interface. The term ∇ s σ links the surface tension to the the Marangoni effect due to change of the surface tension. The unit normal vector can be computed using the color function

n = n |n| = ∇c |∇c| (3.27)

Numerical Simulation Tool

In this thesis, the second publication is done 3Extensible Markup Language (XML) is a meta-language that defines a set of rules for encoding documents in a format that is both human-readable and machine-readable.

4Eclipse is an Integrated Development Environment (IDE) written mostly in Java and is used for developing Java applications, but or applications in other programming languages via plug-ins.

geometry and assigning multiple cores to solve PDE's for 1D and 2D problems. This code uses PETSc5 libraries including MPICH6 and LAPACK.7

SIPER Algorithm

The main algorithm of SIPER code is based on the following steps:

1. Via an input file the necessary simulation data is read in.

2. Initialisation Time Loop starts here 3. Search for Neighboring Particles (After each particle movement) 4. Forces predictor step (in the predictor-corrector scheme, the integration step starts here) 5. Computing v * the temporary velocity 6. Forces corrector step (in the predictor-corrector scheme, the integration step ends here) 7. Solving Pressure-Poisson equation using Stabilized BiConjugate Gradient (BiCGStab) method [9] 8. Computing the velocity v and the pressure p at each grid point 9. Finalisation Time Loop end here 5Portable, Extensible Toolkit for Scientific Computation (PETSc), is a suite of data structures and routines developed by Argonne National Laboratory for the scalable (parallel) solution of scientific applications modeled by PDEs.

6MPICH is a high performance and widely portable implementation of the Message Passing Interface (MPI) standard.

7Linear Algebra PACKage (LAPACK/BLAS) is written in Fortran 90 and provides routines for solving systems of simultaneous linear equations, least-squares solutions of linear systems of equations, eigenvalue problems, and singular value problems.

To improve the condition number (a measure for sensitivity of the output to input variations) of the SPH method the algebraic multigrid preconditioner BOOMERANG from the HYPRE library8 [40] is used.

SPH grid

SIPER uses Smoothed Particle Hydrodynamics methods to solve the governing partial differential equations known as the Navier-Stokes equations. In SPH method, particles are used to discretize the domain. To convert the partial differential equations into algebraic equations using numerical methods, one needs to split the whole domain into multiple points or particles known as grid points, by choosing an appropriate number of grid points across the domain. This decomposition of the domain is known as a SPH grid. SIPER uses Cartesian grids, which are composed of points in one or two dimensions. In case of multiphase problems, each grid point will be assigned an ID that represents the phase of that grid point. This allows us to compute the solutions of differential equations at each point distributed uniformly across the domain.

Field variables

SIPER uses the notion of a Field in relation to the concept of fields in physics. It provides scalar, vector, and tensor fields defined at each grid point. A vector field is a collection of n scalar fields, similarly a tensor field is a collection of n vector fields, where n is the number of dimensions.

8Library of High Performance Preconditioners (HYPER) is a software library for the solution of large, sparse linear systems on massively parallel computers.

According to Morris et al. [START_REF] Morris | Modeling Low Reynolds Number Incompressible Flows Using SPH[END_REF] formulation, if the assumption of incompressibility is applied, and thus the viscosity formulation simplifies to

∇τ = ν∇ 2 v (3.28)
which with SPH discretisation yield to

P i = ∑ j η i + η j m i m j ρ i ρ j 1 r i j ∂W i j ∂ r i j v i j (3.29) U i = -∑ j η i + η j 2 m i m j ρ i ρ j 1 r i j ∂W i j ∂ r i j v 2 i j
(3.30)

Parallel SPH code using MPI

Highly scalable, parallel computer programs have become indispensable tools for the advancement of numerical research. Researchers are more hopeful than ever before, in tackling complex and huge engineering and scientific problems due to the availability of required computational resources in past twenty years. In terms of the hardware resources the progress seems well ahead of the progress in application (software) domain. There are many legacy codes which are still relevant in today, but they lack the modern approach to efficiently use the available hardware resources. Some of the computationally huge problems which cannot be solved without the modern high-performance computing (HPC)

technologies include: microfluidic phenomena, astrophysical analysis, heat transfer analysis, electrohydrodynamics modeling, plasma physics etc. It is unimaginable to work out these problems on single core computers, parallel computing is the only way forward for such problems.

In general, parallel computing refers to solving parts of a problem simultaneously on multi-core computing machines. A problem which can be broken into multiple smaller and discrete parts which can be solved independently, makes a good candidate for parallel computing. The discrete parts of the problem are solved on different computing cores and after finishing they are synchronized to provide the solution of the whole problem. Parallel computing offers several benefits to users: saving time and money, solving complex and large problems, multi-tasking etc. Advancements in HPC have provided another way -'computational science', of doing science along with the classical branches of experimental and theoretical sciences. Computational scientists make use of their simulation methods when they are more advantageous and feasible over the classical approaches of theory and experiments.

Three broad areas of parallel computing are hardware, algorithms and software. In hardware, adding more and more cores and providing efficient inter-communication network among cores has increased the parallel nature of computing machines. In algorithmic terms, scientists seek how a problem can be defined by independent physical mechanisms, and, how it can be solved with independent set of mathematical equations. However, a bigger challenge is posed by the inadequate software, which are not fully able to profit by the progress made in hardware and algorithms. In terms of important characteristics, the modern codes should be optimized, portable and future-proof with every evolving HPC technologies.

A code should make optimal use of the hardware properties such as the cache design, vector registers, multiple cores etc. It should be developed with the standard parallel programming models such as MPI, OpenMP, hybrid models etc.

Surface tension treatment

In multiphase flows, the matter is always classified between the bulk phase, the dispersed phase(s) and the region between two homogeneous phases, known as interfacial state. The physical properties of the interfacial phase being different from those of the bulk and dispersed phases, thus giving rise new variables such as interfacial surface tension. In the equilibrium state, interfaces need to be considered where surface-to-volume ratio is large which leads to large ratio of interfacial free energy to total free energy.

As mentioned earlier on this chapter, in a two-phase system for molecules within the continuous phase the resultant force of all intermolucular interactions gets balanced (see Fig2.5). However, for the molecules at the fluid-fluid interface, the force towards of the continues phase is higher in magnitude, thus pulling the molecules in this direction and creating a curved shape with two radii . The

Color Function

The idea of color [START_REF] Shadloo | Improved Multiphase Smoothed Particle Hydrodynamics[END_REF] function is to assign a value to each phase of the multiphase system so to simplify the value attribution to each field variable and avoid confusion in the calculation.

Color function can be then used specify the interface with some conditions that are to be satisfied.

Here the smoothed color function of the particle i of phase A is defined as :

C i = ∑ j W i j c j ∑ j W i j (3.31) c j =        1, fluid A 0, fluid B
This technique is in particular useful in smoothing the discontinuities between the quantities at each phase using an optional weighting function, W i j .

Furthermore, this approach provides a clear definition for the volume fraction of each fluid i.e. C A i = C i and C B i = 1 -C i defines the volume fraction corresponding to the fluid A and fluid B, respectively such that : ∑ n C n i = 1 for all n phases.

Particle Shifting

The general idea of the particle shifting technique is to slightly modify the equation of evolution of the particle position ( Dx Dt = u) and thus force the particles not to evolve exactly along the trajectories. This allows to artificially break the Lagrangian structures that are classically observed with the SPH method. To this end, it required to add an extra velocity term δ u, called shifting velocity, to the Lagrangian velocity u :

u * = u + δ u, (3.32)
where u * is the modified velocity that will be used in the equation for the evolution of the position of the particles ( Dx Dt = u). Equivalently, equation 3.32 can also be formulated in displacement form by time integration:

r * = r + δ r. (3.33)
It is worth noting that to preserve the Langrangian nature of the SPH method, in case of application of a modified velocity, one should make sure that the shifting velocity δ u is sufficiently small enough compared to the velocity u.

Projection method

The continuity and momentum equations can be discretized according to the first-order projection time integration scheme. For every time step ∆t 1. Evaluate preliminary velocity

u * i = u n i + ∆t(-(u n i .∇)u n i + ν∇ 2 u n i ) (3.34)
2. Using the preliminary velocity, build and solve pressure Poisson equation P n+1

∆P n+1 i = ρ ∆t ∇.u * i (3.35)
3. Using the final pressure and pressure gradients, correct the velocity and project it on to the divergence-free velocity field

u n+1 i = u * i -( 1
ρ ∇p n+1 ) i ∆t (3.36)
This method is relatively stable but it's not so accurate. This is because there is a numerical boundary layer occurrence that makes the velocity order O(∆t) and the pressure order O(∆t 0.5 ). This numerical boundary layer issue could limit the accuracy of the model when mixed with an interpolation scheme that is of spatially fourth-order accuracy. One can also use the second-order projection time integration scheme proposed by [START_REF] Hu | An Incompressible Multi-Phase SPH Method[END_REF] 1. Evaluate preliminary velocity

u * i = (4u n i -u n-1 i + 2∆t(-(u * i .∇)u * i -1 ρ ∇P n i + ν∇ 2 u * i )) 3 (3.37)
2. Using the preliminary velocity, build and solve pressure Poisson equation for pressure increment q n+1

∆q n+1 i = 3ρ 2∆t ∇.u * i (3.38)
3. Evaluate the pressure at step n+1

P n+1 i = q n+1 i + P n i -µ∇.u * i (3.39)
3.8 Particle Shifting 73 4. Using the final pressure and pressure gradients, correct the velocity and project it on to the divergence-free velocity field at step n+1

u n+1 i = u * i -( 1
ρ ∇q n+1 ) i 2∆t 3 (3.40)
The second-order method is even less stable due to pressure-velocity colocation but it's more accurate. The velocity order O(∆t 2 ) and the pressure order O(∆t).

Summary

To sum up, the code which is used in this project is an in-house program including but not limited to many subroutines such as kernel, grid, ,and interface simulations. Nevertheless, there are numerous studies related to the mathematical foundation of SPH including its numerical stability, neighbouring points search algorithm, dissipation and so on and so forth which are not mentioned in this chapter [START_REF] Shadloo | Improved Multiphase Smoothed Particle Hydrodynamics[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF][START_REF] Monaghan | Smoothed particle hydrodynamics[END_REF].

Chapter 4

Results

Introduction

In this chapter two cases of two-phase systems: EHD and Thermocapillary simulation and a coupled case: EHD-TC simulations, are presented. The EHD and Thermocapillary simulations are demonstrated using the published or submitted papers. These cases are analyzed qualitatively and quantitatively, namely analysis of the effect of the numerical resolution and the surface tension on EHD deformation and thermocapillary droplet migration is described. Finally, the evolution of the droplet subject to coupled electric and thermocapillary forces is demonstrated.

Results: Numerical Simulation of Multi-phase Electro-Hydrodynamics Systems

The numerical simulations of complex multi-physics electrohydrodynamics (EHD) problems using ISPH method was first studied by Shadloo [START_REF] Shadloo | A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid[END_REF] were the first group to suggest and Results implement a model for such problems. The code was first validated with simple EHD deformation of droplets suspended in a neutrally buoyant Newtonian fluid. We demonstrated a step-by-step validation of the SPH method for each individual force using the same methodology to answer the existing gap in the literature of this topic. The following paper extends the limits of understanding about the simulations of multi-phase EHD flows using a simple ISPH method for different cases. In this work, a numerical simulation of multi-phase flows under the effects of an external electric field using various scenarios ranging from low to high deformations, droplet migration, and the effect of shear flow on the droplet's deformation is performed. This study answers the existed research gap, provides a broader perspective into the capabilities of the SPH method for such applications.

Introduction

Predicting the behavior of multi-phase flow systems has attracted for decades the attention of many industries due to their wide ranges of applications in the chemical engineering, aerospace engineering and renewable energy sectors, among others [1][2][3]. In multi-phase systems two or more fluids share interfaces which can deform/ migrate as outcome of exerted forces and constitutional laws. Some applications of multi-phase systems include boiling, condensation, water purification and petroleum refinement processes where these phenomena have been investigated mostly experimentally and theoretically. However, with the ever-increasing power of Computational Fluid Dynamics (CFD) methods, numerical simulations of these systems became of great interest among researchers [4][5][6].

Smoothed Particle Hydrodynamics (SPH) is a relatively recent and promising mesh-less Lagrangian method which discretizes the domain into a set of nodes, known as material particles. These particles can freely move inside the computational domain subject to an external force or particle-particle interactions. Initially introduced by Gingold and Monaghan [7], and Lucy [8] for astrophysics applications, SPH was soon found to be suitable for fluid dynamics problems, where complex geometries [9][10][11], large deformations [12][13][14], multi-phase [15][16][17] and multi-physics problems [18][19][20] are involved. A recent overview for the application of SPH can be found in [21].

One of the most important engineering problems which involves many of above cases is the Electrohydrodynamics (EHD) one, where hydrodynamics of a fluid system is coupled with its response to an external electric field. In EHD problems, one may control the interface between the two fluids (here, the droplet and the bulk fluid) by controlling the flow conditions and fluid properties [22,23]. In such problems, the coupling may lead to a large interfacial deformation (i.e. merge/breakup) or migration. Indeed, EHD is a very complex problem including multi-phase, multi-physics and multiscale phenomena with strong topological changes of the interface shape [24,25]. Although, there are many experimental and theoretical studies available in the literature on the coupled modeling of EHD problems [26][27][28]. Nevertheless, some discrepancies between experiments and analytical data still exist [28]. As such, numerical simulations have been developed to tackle these difficulties and provide insight into EHD problems.

Considering the numerical simulations of EHD using SPH method, Shadloo et al. [29] were the first group to provide a model for such problems. They validated their code with the simple EHD deformation of droplets suspended in a neutrally buoyant Newtonian fluid. Rahmat et al. have proposed a multi phase ISPH method based on the lubrication theory and the drainage model to simulate droplet electro-coalescence for wide ranges of simulation conditions [30,31]. Rahmat et al. [32] also provided the first simulation results for the Rayleigh-Taylor instability under the combined effect of electric field and gravitational forces. Yet, step-by-step validation of the SPH method for each individual force using the same methodology is not well-documented. Additionally, numerical simulation of a multi-phase flow under the effects of an electric field using various scenarios ranging from low to high deformations, droplet migration, and effect of shear flow on the droplet's deformation would provide a broader perspective into the capabilities of the SPH method for such applications. To this end, this article aims at introducing a mesh-less numerical methodology, i.e. Incompressible SPH (ISPH) approach, to deal with such complex problems. Additionally, we verify the applicability of some of the used algorithms for a range of problems including hydrodynamic, capillary, gravity, shear and EHD forces.

This article is organized as follows: First, we introduce the mathematical formulation of the SPH method as well as the numerical discretization scheme. Then, we incorporate the governing equations of the multi-phase system including the conservation equations for mass, momentum and electrical charges in a Lagrangian form. Thereafter, a code validation and numerical convergence study is asserted in the absence of electric field. Numerical results cover solutions with and without electric field sections. Additionally, the effect of surface tension through Laplace law, the effect of gravitational force, and the Couette flow for a multi-phase system are examined and validated against analytical solution and available numerical data in the literature. Finally, conclusions are provided in the last section.

Mathematical formulation of SPH

The idea of SPH comes from the fact that any field variable f (x) can be calculated by an exact mathematical relation as

f (x) = ∫ Ω f (x ′ )δ(x ′ )dx ′ . (1)
Upon approximating δ(x ′ ) by an interpolation function W (xx ′ , h), this equation can be formulated as

f (x) = ∫ Ω f (x ′ )W (x -x ′ , h)dx ′ , (2)
where x and x ′ are the position vectors and h is the smoothing length. In our case, h = ζ dx where ζ = 1.6 is a constant value, and dx is the initial particle spacing. The interpolation function, also known as smoothing function or kernel function, should have, among others, the following properties [33] -Normalized over the domain

∫ Ω W (x -x ′ , h)dx = 1.
(3) -Produces δ function for a small enough smoothing length

lim h→0 W (x -x ′ , h) = δ(x ′ ). ( 4
)
-Remains monotonically decreasing throughout the entire domain.

-Has a compact support, meaning that for |x -

x ′ | > kh W (x -x ′ , h) = 0. (5)
-Is a symmetric function.

Initially, the kernel functions were defined such that each particle should have interactions with all others [7]. By introducing the concept of Neighboring particles, 1 kernel function affect only a compact support around it were substituted (see Eq. ( 5)). Depending on the smoothing length parameter h, only a few number of particles in the entire space affect the approximated value of the kernel function (around 25 to 35 in 2D). In the current work, a cubic spline kernel function is used both for the bulk fluid and the interface modeling while taking harmonic average.

W ij = A ⎧ ⎪ ⎨ ⎪ ⎩ 2/3 -(r/h) 2 + 1/2(r/h) 3 r/h ∈ [0, 1) 1/6(2 -r/h) 3 r/h ∈ [1, 2). 0 r/h ≥ 2 (6)
Hereafter, W (x-x ′ , h), will be shown by W ij and A = 15 7π h 2 . Also, i, j, and r represent the index of the particle of interest, the index of its neighbors, and the smoothing radius.

To calculate the SPH gradients, one can show that it is sufficient to differentiate the kernel function W (x-x ′ , h). In other words, in SPH there is no need to differentiate the field function f (x); instead one can differentiate the kernel function. The latter is one of the fascinating features of the SPH method which distinguishes this method from other mesh-based techniques. In this work, we use an improved version of the first derivative, presented in [34] as

∂f m i ∂x k i a kl i = ∑ j 1 ψ j (f m j -f m i ) ∂W ij ∂x l i . (7)
Also, the derivatives for vectorial and scalar quantities are calculated, respectively, as follows:

∂ 2 f m i ∂x k i ∂x k i a ml i = 8 ∑ j 1 ψ j (f m i -f m j ) ∂W ij ∂x l i r m ij r 2 ij , (8)
and

∂ 2 f i ∂x k i ∂x k i (2 + a kk i ) = 8 ∑ j 1 ψ j (f i -f j ) ∂W ij ∂x k i r k ij r 2 ij , ( 9
)
where ψ is the particle number density and a kl i represents a corrective second rank tensor to avoid particle inconsistencies [9].

Governing equations

Assuming an immiscible two-phase Newtonian, viscous, incompressible, isothermal fluid system, the corresponding mass and momentum conservations in a Lagrangian formulation are given as follows

Dρ Dt = -ρ∇. ⃗ V , (10)
and

ρ D ⃗ V Dt = ∇.T + ⃗ f b + ⃗ f s + ⃗ f e , (11)
where, ρ is the fluid density, D Dt is the material time derivative operator, 2

∇. ⃗

V is the divergence of the velocity vector,

T is the total stress tensor which is defined as T = -pI + τ where p is the static pressure, I is the identity matrix and

τ = µ(∇ ⃗ V + (∇ ⃗ V ) T
) is the viscous dissipation term for µ being the dynamic viscosity. Additionally, ⃗ f b = ρ ⃗ g is the body force due to gravity and ⃗ f s is the surface tension which can be calculated using the volumetric force proposed by Brackbill et al. [35], so called the Continuum Surface Force (CSF) method, as

⃗ f s = γ κ ⃗ nδ s . ( 12
)
Here, γ is the surface tension coefficient, κ = -∇.⃗ n is the interface curvature, ⃗ n = ∇C |∇C| is the unit vector normal to the interface, and δ s = |∇C| is surface Dirac-delta function, and finally, ⃗ f e is the electric field force.

To avoid sharp discontinuities at the interface, the smoothed color function of the particle i is defined as

C i = ∑ j W ij c j ∑ j W ij , ( 13
)
1 Particles that are located within the range of the kernel function with respect to the particle of interest. Outside of this range, the kernel function has already dropped to zero. 2 The material time derivative is a directional time derivative for a fixed point.

where the color function c assigns a unit value to one phase and zero to the other phase in a two-phase system such that As mentioned before, in this study the electrostatics and the hydrodynamics are coupled together. This coupling is achieved through the Maxwell stress tensor. Maxwell equations provide a mathematical framework for the interaction and the connection between the electric and the magnetic fields [36]. Here, the EHD part of the system can be regarded as quasi-static model, and dynamic currents values are so low, hence the induced magnetic field effects are negligible. Therefore, the contribution from the induced magnetic field is neglected. Consequently, the volumetric electric force can be written as

⃗ f e = ∇.T E . ( 14
)
In case of an application of the external electric field on a multi-phase fluid flow, this new term for the electric force, will be added to the right hand side of the momentum equation (see Eq. ( 11)), where the Maxwell's stress tensor defines as

T E = ⃗ D ⃗ E - 1 2 ( ⃗ D. ⃗ E)I, ( 15
)
where ⃗ E is an external electric field, ⃗ D = ϵ ⃗ E is the dielectric displacement vector, and ϵ is the electrical permittivity. Also, based on the Gauss' law [36] ∇. ⃗

D = q v , (16)
where q v is the free electric charge density. Application of Eqs. ( 15) and ( 14) will result in

⃗ f e = q v ⃗ E - 1 2 ⃗ E. ⃗ E∇ϵ. ( 17
)
In this work both fluids are considered to be leaky dielectric, (i.e. electric relaxation time is much shorter compared to its viscous counterpart or t e ≪ t v ).

Time integration

We apply a predictor-corrector scheme to advance the governing flow equations in time considering a first-order Euler approach. The time-step is selected based on Courant-Friedrichs-Lewy (CFL) condition in which ∆t = C CFL h/V max , with V max being the largest magnitude of particle velocity and the C CFL is the constant taken as 0.25. During the predictor step, we first advance all the variables to an intermediate value denoted by (*), from the variables' value at the nth time-step denoted by superscript (n), as

⃗ r * i = ⃗ r (n) i + ⃗ V (n) i ∆t + δ⃗ r (n) i , (18)
⃗

V * i = ⃗ V (n) i + RHS ρ (n) i ∆t, (19)
ψ ( * ) i = ψ (n) i -∆tψ (n) i (∇. ⃗ V * i ). ( 20
)
RHS denotes the right hand side of Eq. ( 11), ψ i = ∑ j W ij , is the number density associated with the particle of interest i, which is calculated from the summation of kernel function at all neighboring particles j, δr i is the artificial particle displacement, defined as

δr k i = α ∑ N j (r k ij /r 3 ij )r 2
i,o V max ∆t, and its constant α is set to 0.05 according to [34].

These intermediate values will then be used to solve the Poisson equation which gives the pressure value at the next time-step (n + 1). Using this pressure, new velocity and displacement vectors are updated as following ∇.

(

1 ρ * i ∇p (n+1) i ) = ∇. ⃗ V * i ∆t , (21)
∇. ⃗

V (n+1) i = ⃗ V * i - 1 ρ * i ∇p (n+1) i ∆t, ( 22
) ⃗ r (n+1) i = ⃗ r (n) i + 0.5( ⃗ V (n) i + ⃗ V (n+1) i )∆t + δ⃗ r (n) i . ( 23
)
Fig. 1. Schematic of the test case for validation and numerical convergence test, bubble rising as well as bubble deformation under the effect of electrohydrodynamics (EHD). For the first test case ⃗ g = 0 and ⃗ E = 0, for the second one ⃗ g ̸ = 0 and ⃗ E = 0 while for the third one ⃗ g = 0 and ⃗ E ̸ = 0.

Fig. 2.

Comparison of (left) the pressure jump across the droplet interface for three particle resolutions and (right) its comparison with the theoretical pressure jump, i.e. Laplace's law, for different cases.

Results

Validation and convergence

To ensure a suitable particle resolution based on the numerically computed pressure jump across the interface, in to capture the physical jump across the interface. Once more for the reported simulations, the relative error is less than 1% when the numerical results are compared to the Laplace's law.

It is noteworthy to discuss the reason behind the different values obtained for the theoretical and numerical pressure jump at the interface. As mentioned before, the pressure and other flow field variables (velocity, color function, etc.) are approximated using the numerical smoothing scheme which converts the sharp values at the interface to smoother ones resulting in a loss of accuracy and introduction of spurious oscillations near the surface of the droplet. Furthermore, it is found that these spurious currents are generated because of an inappropriate evaluation of the curvature of the circular droplet due to unreliable values for the unit normal vector (⃗ n = ∇C |∇C| ) in the surface tension force calculation [37].

Bubble rising

In this section, the ISPH method is applied to test and to validate the simulation of the bubble rising problem due to the gravitational force. The computational geometry for this test case is similar to the one shown in Fig. 1 except that the domain size is increased in the normal direction (orthogonal), i.e. x/r = 6 and y/r = 12, and the bubble is initially placed such that x o /r = 3 and y o /r = 2.4. The grid resolution is set to 240 × 480 in x and y direction, respectively. The velocity boundary conditions are set to be free slip for BC -X 1 and BC -X 2, and no slip for BC -Y 1 and BC -Y 2. Also, pressure boundary conditions are set to be Dirichlet with a constant value at BC -Y 2 and Neumann for the other three boundaries (∇p • ⃗ n = 0) where ⃗ n is normal direction to the given boundary. Here, both bubble and bulk phases are set to have stationary conditions at initial time step. The bubble starts to rise during the simulation due to the gravitational forces. This problem can be characterized by Reynolds and Bond numbers defined as following:

Re = ρ 2 g 0.5 (2r o ) 1.5 µ 2 , (24)
and

Bo = ρ 2 g(2r o ) 2 γ . ( 25
)
respectively.

For the first simulation, a case with low density and viscosity ratios and high surface tension is considered, where Additionally, the gravity is set to be ⃗ g = -1 [m/s 2 ] in y direction such that it produces Reynolds number Re = 35, and the Bond number Bo = 10 for this case. The time snapshot of this test case for dimensionless times of t * = 0.5, t * = 18.5, t * = 37 and t * = 64.8 are shown in Fig. 3. As observed in this figure, the bubble starts to rise straight upwards due to the gravity, while its velocity increases from zero to 0.36 [m/s] and remains constant until it feels the pressure coming from the stationary upper-wall. Additionally, the bubble shape is changing from circular shape to a quasi elliptical one due to hydrodynamic pressures on its tip. The final shape comes from the competition among surface tension, gravitational, and viscous forces. As it is observed, the shape remains unchanged after some time steps which is the main reason for the bubble's almost constant terminal velocity. Fig. 4 shows the mean migration velocity and the position of the droplet with respect to dimensionless time t * = t √ (g/D), which are in agreement with the results of [31] and [38]. Here, the velocity gradient near the stationary wall starts to deviate at the final times which could be due to the confinement effects.

To show the applicability of the proposed algorithm for capturing larger deformations and breakups, we perform a second bubble rising test case. This time, the computational domain size is x/r = 6 and y/r = 10 with the particle resolution of 240 × 400. The bubble is initially placed at x o /r = 3 and y o /r = 2. Here, all four boundaries have no slip boundary condition for velocity. However, the pressure boundary conditions are kept the same as before. The gravity is selected as ⃗ g = -1 [m/s 2 ] in y direction, while the surface tension coefficient is set to be γ = 20 [N/m]. Additionally, the density and viscosity ratios are ρ 2 /ρ 1 = 1000 and µ 2 /µ 1 = 2.828/10 with ρ 1 = 1 [kg/m 3 ] and µ 1 = 1 [Pa s], respectively. These choices are for mimicking the test case presented in [15] in order to produce the fluid flow system with Reynolds and Bond numbers of Re = 1000 and Bo = 200, respectively.

The snapshots of our current simulations are illustrated (middle) in Fig. 5 for the dimensionless times between t * = 3.2 and t * = 5.6 with a time increment of ∆t * = 0.4. These snapshots are compared to their experimental (top) and Volume of Fluid method (bottom) counterparts, presented in [39] and [40], respectively. As can be seen, the proposed ISPH approach can predict the large deformation and bubble breakup as accurate as its well established mesh based method. The presented Volume of fluid (VoF) method uses a hybrid VoF-level-set method [40] to accurately capture the interface.

EHD droplet deformation

In this section, we consider a suspended circular droplet under the effect of an external applied electric field. The schematic of the computational domain is similar to what is presented in Fig. 1 with an increment in the size of the domain. Here, we double the domain size in each direction, i.e. x/r = 8 and y/r = 8, in order to reduce the confinement effect. In Fig. 6-left, a particle resolution of 240 × 240 is used with a circular droplet initially placed at the center of Comparison of bubble rising and break up due to gravitational and hydrodynamic forces using experiment [39], SPH method and Volume of Fluid method [40], respectively from top to bottom rows. Here, the density ratio ρ 2 /ρ 1 = 1000 and viscosity ratio µ 2 /µ 1 = 10 are applied. Reynolds and Bond numbers, as previously defined, are Re = 1000 and Bo = 200, respectively. The dimensionless time is defined as t * = t √ g/r o , starting from t * = 3.2 at the very first frame on the left, up to t * = 5.6, with a time increment of 0.4 per frame. The deformation of a suspended circular droplet under such conditions is a commonly utilized test case for validation of a EHD solver, where two theories are available in the literature. Taylor et al. [41] estimate the droplet deformation D T as

D T = 9f dT E 2 o ϵ 2 r o 8(2 + R) 2 γ , (26)
where f dT is the discriminating function evaluated as For the same problem, Feng [42] suggests the following relation

f dT = R 2 + 1 -3.5S + 1.5R. (27)
D F = f dF E 2 o ϵ 1 r o 3(1 + R) 2 Sγ , (28)
where f dF is estimated from

f dF = R 2 + 1 -3S + R. (29)
In Eqs. ( 26) and ( 28), r o is the initial droplet radius before its deformation and E o is the electric field intensity in vertical direction which is calculated from E o = (ϕ + -ϕ -)/H, H being the domain height. Additionally, the permittivity ratio and the conductivity ratio of droplet to the bulk are called S and R, presented as

S = ϵ 1 ϵ 2 , R = σ 1 σ 2 , (30)
where ϵ and σ are the electrical permittivity and conductivity, respectively. Also the subscripts 1 and 2 show, droplet and bulk medium properties, respectively. Another point in the theory of droplet deformation is to investigate the velocity recirculation vectors inside and outside of the droplet when a vertical electric field is applied. The relative differences in the electric permittivity and conductivity of both constituent phases define the direction of the flow rotation in either phase. This is shown for two cases in Fig. 6.

On the left side test case is adopted for S = 0.5 and R = 0.05 with the electrical permittivity of the droplet (ϵ 1 ) being 0.5 [F/m] which is half of that of the bulk fluid. Also, the electrical conductivity of the droplet (σ 1 ) is set to 150 [S/m] which is three-times more than the background fluid. On the right side of Fig. 6, the test case has S = 0.5 and R = 3 with ϵ 1 = 0.5 and σ 1 = 1. As can be seen in the left sub-figure, the re-circulation zone in the first quarter (i.e. the top-right quarter inside the droplet) of the droplet orients clockwise. This should be the case for S > R and is consistent with the results of [29]. The opposite flow circulation pattern should be expected for the case of S < R as it also presented on the right side of the same figure.

Additionally, Eqs. ( 27) and ( 29) define the sign of Eqs. ( 26) and ( 28), respectively. The positive sign, so called prolate deformation, indicates that the droplet is elongated in the direction parallel to the electric field. The positive sign, so called oblate deformation, shows the droplet elongation in the opposite direction. For the comparison with the simulation results, we define the numerical deformation as

D N = A -B A + B , (31)
where A and B are the elliptic droplet diameters, parallel and perpendicular to the direction of the external electric field, respectively, at the steady state condition. When this parameter is equal to zero, the droplet keeps its initial circular shape. On the other hand, more deviation from zero indicates more deformation from its initial shape. Fig. 7 provides comparison between the current ISPH results and the two aforementioned theories for multiple cases. As can be seen, numerical data reasonably follow the available theories. However, in most of the cases, an over-prediction is reported by the simulations. Some might be some possible reasons for such behavior can be mentioned: (i) in the theory it is assumed that the droplet remains circular even after applying the electric field. This means the change in the curvature is not integrated in the deformation equation, but only the surface tension coefficient. (ii) Another reason might be due to the confinement effect. In theories, it is assumed that the droplet is suspended in an unbounded domain for simplicity. However, providing such domain numerically, or even with twice larger computational domain, is very expensive computationally. Finally, the hydrodynamical properties of droplets such as density and viscosity are not taken into account theoretically and the problem is considered to be inviscid.

Couette flow

This section investigates the deformation of a droplet suspended between two parallel plates subjected to a constant shear. The flow configuration, known as Couette flow, is presented in Fig. 8. In this case, the flow is driven by viscous forces or pressure gradients [43]. Different cases with and without an external electric field with the magnitude of unity perpendicular to the flow direction are simulated. The computational domain consists of a rectangle box with the size of x/r = 16 and y/r = 4 discretized by a set of initially equally spaced 400 × 100 particles, arranged in a Cartesian grid. The two-phase system contains a droplet with Additionally, particles inside the droplet are initialized to be at rest, while background fluid particles having undisturbed Couette flow velocity. The simulations are performed for a range of electrical permittivity and electrical conductivity ratios shown as (S, R), while neglecting the gravity and keeping the surface tension coefficient constant γ = 0.02 [N/m]. The droplet radius is a quarter of the distance between two parallel plates (i.e. r = H/4). The dimensionless Reynolds, Weber and Electroinertial numbers, respectively, are as follows

Re = ρ 2 U o r o µ 2 , We = ρ 2 U 2 o r o γ , Ei = ρ 2 U 2 o ϵ 2 E 2 o , (32)
where Re = 1, We = 0.2, and Ei = 50. Fig. 9 shows the time evolution of the droplets' deformation under the same external electric field and the shear stress conditions, but for different working fluids having different electrical permittivity and conductivity ratios (i.e. only a change in S and R is considered here). In this figure, the dimensionless time is defined as t * = tU o /r o , while the droplet deformation is calculated from

D f = L max -L min L max + L min , (33)
where L max and L min are major and minor droplet diameters, respectively. As can be seen in this figure, the droplet deformation increases when an external electric field is applied regardless of S and R. Additionally, at the constant electrical permittivity ratio, larger deformations can be achieved by an increment in the electrical conductivity ratio when S > R. However, for the similar condition (i.e. S = cte.), smaller deformations are seen by an increment in the electrical conductivity ratio when S < R. Similarly, at constant electrical conductivity ratio larger electrical permittivities result in larger droplet deformations for R < S, while decreases the same for S < R. It is noted that the (5.0, 0.2) and (5.0, 0.5) test problems did not reach a steady profile during the simulation time, here fixed at t * = 1 due to large computational costs.

Fig. 10 provides a comparison on the interface shape in the absence (in blue) and the presence of the electric field (in black) at this time. Following the previous figure, it can be seen that the droplets are more slender for larger deformation factors.

It is notable that the angle between the major axis of the elliptic droplet and the stream-wise direction becomes smaller with an increment in the conductivity ratio, while the cases of (5.0, 0.2) and (5.0, 0.5) are immediately distinguishable due to their larger deformations. This elongation is due to the suppression of surface tension forces which is the result of the larger electric field force at higher electrical conductivity ratios. It is also noted that the droplet is no longer elliptical and is suspected to have the breakup at larger time steps.

Finally, the time snapshot of an extreme test case with the electrical permittivity ratio of S = 10 and the electrical conductivity ratio of R = 0.2 is illustrated in Fig. 11 to show the ability of the presented method to capture very large deformations with the interfacial topological change. Here, the interface is represented by red color, while the velocity streamlines and electrical vectors are represented by blue and black arrows, respectively. In this case, the circular droplet becomes elliptical soon after the start of the simulation. The elliptical interface is elongated in the stream-wise axis direction due to the four re-circulation zones in the bulk flow just next to the interface. Soon after that, pairs of recirculation merge with each other at the both tips of the droplet and cause the creation of the third re-circulation in its center. By time, the droplet gets folded in four different places and new re-circulation zones appear close to the interface which makes the droplet very susceptible to breakup. As can be seen, this is a promising test problem to show the ability of the present ISPH code to treat the complex multi-phase fluid behavior under an extreme EHD conditions with large interfacial deformation.

Summary

In this work, we presented an effective multi-phase Incompressible Smoothed Particle Hydrodynamics (ISPH) approach to simulate complex multi-physics electrohydrodynamics (EHD) problems. We showed a step-by-step validation of the multiphase code for surface tensions, hydrodynamic forces and electric forces, respectively, by solving Laplace's law, bubble rising and buoyant droplet deformation under an applied electric field problems. Results are validated either against available analytical or numerical results. An overall satisfactory agreement was found. Finally, we presented, for the first time, results of droplet deformation under sheared Couette flow with external electric field for a range of simulation parameters. Different parameters such as time resolved topological changes, droplet deformation magnitudes as well as velocity field and electrical potential vectors were presented and compared with each other. It was shown that the current ISPH approach is able to be easily adopted for different multi-physics problems. It is also capable of predicting large interfacial topological changes such as folding and breakup. In future, our strategy would be to include more complex transport and multi-physics phenomena.

Results: Coupled EHD-TC Droplet Migration

Having that interfacial flows are important phenomena in many engineering processes and to extend the application of SPH method for complex transport and multi-physics phenomena, Hopp-hirschler [START_REF] Hopp-Hirschler | A smoothed particle hydrodynamics approach for thermo-capillary flows[END_REF] used an ISPH model for thermo-capillary driven flows. Such flows are driven by the gradient of surface tension. The model is based on CSF approach and includes Marangoni forces. To increase the accuracy of the results, density-invariant divergence-free (DIDF) [START_REF] Hu | An Incompressible Multi-Phase SPH Method[END_REF], corrected SPH [14,[START_REF] Shepard | A Two-Dimensional Interpolation Function for Irregularly-Spaced Data[END_REF], and particle shifting [START_REF] Shadloo | A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid[END_REF][START_REF] Shadloo | Improved Incompressible Smoothed Particle Hydrodynamics method for simulating flow around bluff bodies[END_REF] techniques are used. We have shown that in mutiphase systems, momentum equation is affected by a non-isotropic pressure field, due to the fluid-fluid interface which give rise to the capillary stress tensor. For immiscible systems, this stress is referred to as the surface tension. The normal component of this stress can be validated using the Young-Laplace law. Here we investigate a two-dimensional problem at the steady-state.

Introduction

Suspended bubbles or liquid droplets deform in fluid flows and in doing so demonstrate a host of phenomena with high industrial importance. For instance, in petroleum industry, there is a major need to accurately drive and control the demulsification process of crude oil [1,2]. In doing so, conventional experimental techniques, such as heat treatment, electrical field, and membrane separation, require complex and expensive setups to enable significant insight into these complex multi-physics problems [3]. In terms of physical modelling, the simulations of the dynamics of bubble rising including deformation and possible merging or break-down, require a correct treatment of the sharp liquid-gas interfaces with a fine modeling of the surface tension that can lead to large deformations of the interfaces. In this article, we focus on the numerical study of a single bubble deformation under different flow conditions. Since the numerical simulation of two-phase flows is inherently a multiscale problem that needs sophisticated strategies for time and space integration schemes, both small and large scale deformations should be treated simultaneously. Moreover, in a two phase system of leaky dielectric fluids inclusion of the temperature response to electric and hydrodynamic response of the system requires special treatment. The application of an external electric field to a droplet can result in large topology and velocity changes with possible break-down and coalescence. Understanding of the underlying principles of Electrodynamics (EHD) can be applied to better control and predict the motion and deformation of droplets.

From a fundamental point of view, most of the simplest multi-physics models cannot be theoretically solved in their closed form, hence the interest in numerical solution arises [4,5]. The classical methods such as Finite Difference [6], Finite Element [7], Finite Volume [8], Level Set Method [9], and combined methods [10] cover the majority of the mesh-based numerical strategies in which the accuracy of the results is highly dependent on the numerical aspects (i.e., numerical methods, time step size and mesh refinement) [11]. Numerical methods can be classified with regards to their modeling approach. In the Eulerian multi-fluid approach, each phase is considered as an interpenetrating continua from which an ensemble averaging of the multi-phase Navier-Stokes equations is calculated [12]. On the other hand, the Lagrangian mesh-free methods notably Smoothed Particle Hydrodynamics (SPH) [13] offers a promising and flexible framework in modeling complex coupled multi-phase fluid problems. As listed in reference [14], among several mesh-free methods, SPH inherently provide notable efficiency in calculating partial derivatives [15] by considering particles which remove the necessity for mesh generation and refinement.

From physical point of view, thermocapillary instability (TC) is caused by inhomogenities in interfacial tension in multi-phase systems. This in-homogeneity is a result of a thermal gradient at the interface introduced by surfactants or temperature variations. Here, we consider the variation of temperature on the surface of the micro-droplet. Temperature variation creates a non-uniform surface tension, which results in interruption in the balance of forces and the introduction of a new shear stress on the surface of the droplet. The imposed strains by the continuous phase alter the structure of the particles of the disperse phase. In particular, coupling the temperature gradient with an electric field leads to a bubble destabilization and deformation. The electrostatic pressure enhances the instability since the electrostatic force on the droplet surface is higher than the capillary pressure. On the other hand, the capillary pressure affects the fluctuations of the free surface, which lead to diminishing instabilities. The instability grows when the electrostatic pressure is larger than the capillary pressure. As a result, the bubble starts to deform or migrate.

There are few experimental and theoretical studies available in the literature of EHD-TC coupled problem. In principle, a droplet in such system evolves in order to reduce the total free energy of the system. The total energy is defined as the sum of internal and the kinetic energies. At steady state, the final shape of the droplet is reached at the lowest interfacial energy level. The linear stability analysis shows a negative correlation between the temperature and the interfacial tension. Regarding the one-dimension thin liquid films, the EHD-TC forces lead to the creation of smaller structures (eddies) [16,17]. Nevertheless, some inconsistencies between the experimental and theoritical approaches are reported because of the electric breakdown effects, when a sufficiently high voltage is applied [18]. Therefore, the current paper concerns the study of complex flow physics including thermocapillary phenomena (Marangoni forces) and electrohydrodynamics by means of SPH simulation. The current work is an extension of our previous numerical study in which an electric field is coupled [19].

This article is organised as follows; the governing equations of a multi-phase system with thermal gradient and electric field are first presented, before introducing the SPH method and the associated space and time discretization schemes. Afterwards, numerical convergence studies are presented. Then we discuss the surface tension, the electrohydrodynamics and the thermo-capillary effects, separately. Finally, the results for electrohydrodynamics-thermo-capillary bubble deformations are presented.

Governing Equations

Mathematical formulations of governing equations of the coupled EHD-TC problem consist of a set of continuity, momentum and energy balance equations.

Continuity Equation

The continuity equation can be derived based on volume flux conservation laws as

Dρ Dt + ρ∇ • v = 0 (1)
where ρ, v and Dρ/Dt represent the density, the averaged fluid velocity and the material derivative operator, respectively. For incompressible flows, the continuity equation reduces to ∇ • v = 0.

Momentum Equation

Assuming that the total stress tensor is symmetric, the momentum equation for a two-phase Newtonian, immiscible, and non-reactive, fluid flow with constant electrical permitivity and electrical conductivity can be written as

ρ Dv Dt = -∇p + ∇ • Π v + ∇ • Π c + F g + F e ( 2
)
where p is the pressure, Π v = µ[∇v + (∇v) T ] is the viscous stress tensor with the dynamic viscosity µ and transpose operator T, Π c is the capillary stress tensor, F g is the body force due to gravity, and F e is the electric force. Note that the sharp interface limit of the capillary stress tensor is defined as ∇ • Π c = (σχn + ∇ S σ)δ where σ is the surface tension, χ is the curvature, ∇ S is the gradient of the surface tension in tangential direction with respect to the surface, and n is the unit normal vector of the interface. The Dirac-delta function δ is defined to be unity at the interface and zero elsewhere. As a result, the capillary stress tensor is constant inside the bulk. The surface tension force is calculated using the continuum surface force (CSF) model [19,20]. Thus, Equation ( 2) can be rewritten as

ρ Dv Dt = -∇p + µ∇ 2 v + (σχn + ∇ S σ)δ + ρg + F e (3)
To calculate the electric force F e , we apply the electrohydrodynamics theory for leaky dielectric fluids [21,22]. The electric force is obtained from

F e = ∇ • T ( 4
)
where T is the Maxwell stress tensor. For weak electric currents, the magnetic field is negligible because the electric field is assumed to be irrotational (∇ × E = 0). The Maxwell stress tensor T reads

T = DE - 1 2 (D • E) Ī ( 5
)
where E is an external electric field, D = E is the dielectric displacement vector, and is the electrical permittivity. Ī denotes the identity tensor. Using Gauss's law,

∇ • D = q v (6)
where q v is the free electric charge density. Knowing that the gradient of the electric field vector is symmetric and by application of the product rule in differentiation, by taking the divergence of the Maxwell stress tensor and combined use of Equations ( 4)-( 6), the electric field force per unit volume F e can be obtained as [23]

F e = q v E - 1 2 E • E∇ (7)

Conservation of Energy

In the energy balance equation, we only include the heat transfer induced by the temperature gradient (∇T) and we neglect the effect of viscous heating and other source terms. However, the heat transfer from one phase to the other is accounted for. The resulting energy equation can be written as :

DT Dt = κ∇ 2 T ( 8
)
where κ is the thermal diffusivity.

Boundary Conditions

This study is based on impermeable boundary conditions. For temperature and pressure, Dirichlet or Neumann boundary conditions are used. For the gradient of the color only Neumann boundary conditions are applied to avoid influence of the fictitious boundary conditions. The velocity boundary conditions are set to either slip or no-slip conditions. Based on the features of our proposed SPH model (see Section 3), there is no need to consider normal and surface tension boundary conditions. Boundary treatment is done using the implementation of Cummins and Rudman [24] approach at straight walls, also known as mirror particles. According to this approach, the fluid particles are mirrored in every time step across the wall. The image of the particles, distanced to a certain value from the fluid domain, represent the discretized wall particles. Thereafter, the properties of the wall particles are chosen to apply Dirichlet and Neumann boundary conditions. More specifically, any property of an image particle Ψ i and a particle Ψ i (such as velocity, pressure, temperature, color) can be defied such that

Ψ i = 2Ψ wall -δ wall Ψ i (9)
where the Ψ wall is the value of the wall and the sign function is defined as

δ wall = +1
for Dirichlet boundary conditions. -1 for Neumann boundary conditions with Ψ wall = 0 (10) More details can be found in [25,26]. It is more suitable for straight wall conditions (i.e., our currant work) compared to curved walls.

Smoothed Particle Hydrodynamics (SPH) Method

Mathematical Formulation of SPH

Smoothed Particle Hydrodynamics (SPH) method was developed independently by Gingold and Monaghan [27], and Lucy [28] as a truly Lagrangian particle-based method with a superior ability in modeling complex geometries and large fluid flow deformations. Several studies show the wide range of applications of SPH method as presented for example in [15,29]. In SPH, the continuum system is discretized into interpolation points, called particles, that can move freely and carry physical properties such as mass, momentum and temperature. An arbitrary function f (x) can be exactly reformulated as

f (x) = f (x )W(x -x , h)dx (11)
Here, the smoothing (kernel) function W is introduced. The smoothing length h measures the radius of the kernel, while x and x are the position vectors of two different particles. According to Monaghan [30], the kernel function must satisfy certain conditions. It has to be normalized over the whole domain Ω, such as Ω W(xx , h)dx = 1 and contracts to Dirac-delta function δ, so that lim h→0 W(xx , h) = δ(x ).

Moreover, a suitable kernel function should have a compact support; i.e., for every

k ∈ IR + if |x -x | > kh then W(x -x , h) = 0.
This last condition ensures the numerical efficiency in approximating physical variables such as velocity, density or temperature. The discrete form of Equation (11) reads

f (x i ) = ∑ j f (x j )W(x i -x j , h)V j (12)
Here, we use the 4-th order Wendland C2 kernel function [31] given by

W(x -x , h) = C w h d (1 - q 2 ) 4 (1 + 2q) q ∈ [0, 2]. 0 otherwise. ( 13
) W (x -x , h) = C w h d (1 - q 2 ) 3 (-5q) q ∈ [0, 2]. 0 otherwise. ( 14
)
where q is the dimensionless smoothing length, q =

xx h . The normalization constant C w at each dimension d is

C w =                        3 4 d = 1 7 4π d = 2 21 16π d = 3 (15)
Figure 1 illustrates the kernel function W and its derivative with smoothing radius h = 0.75. 13) and (14).

One of the main advantages of SPH over other mesh-less methods is that by starting from Equation (11), one can calculate the derivatives of the function f by means of the gradient of the kernel function. Here, we use two formulations for the first derivative for different conditions. More details on the derivatives in SPH can be found in the following References [32,33].

The first formulation, called the negative formulation, Monaghan [34], is often used for the divergence free velocity since it guaranties that the gradient of a constant field (

f i = f j ) is exactly zero. ∇ f (x) i = ∑ j m j ρ j f j -f i ∇W ij ( 16
)
where W ij is an abbreviated notation for W(x i -x j , h) and i is the particle of interest, among j neighboring particles. The gradient of the kernel function

∇W ij = x i -x j |x i -x j |
∂W ∂x is hence calculated using the analytical expression of W shown in Equation (13).

The second formulation of the first derivative, known as positive formulation, proposed by Monaghan [30], is often used to calculate the pressure gradient. As mentioned in [30] the advantage of this formulation over Equation ( 16) is that due to the symmetric term f i + f j , the conservation of both linear momentum and angular momentum hold when calculating the pressure forces.

1 ρ ∇ f (x) i = ∑ j m j f j + f i ρ i • ρ j ∇W ij (17)

Application of the SPH to the Governing Equations

Present section concerns the discretized formulation of the set of governing equations in the context of SPH method.

As for the continuity equation Hu and Adams [35] proposed a conservative formulation of the discrete mass conservation as following

ρ i = m i ∑ j W ij (18)
In this formulation, ρ i is calculated directly from m i V i with V i being the volume of particle i. The error of the total volume of a wall-bounded system of particles is therefore bounded, since the neighboring particles contribute to the particle density only by affecting the specific volume of particle i [35].

Using the presented first and second derivatives of SPH, the discrete momentum balance can be written as

Dv Dt = -∑ j m j ρ i • ρ j (p i + p j )∇W ij + ∑ j m j ρ i • ρ j µ i + µ j x i -x j |x i -x j | 2 ∇ f (x) i • (v i -v j ) + F g,i + δ[σ i χ i n i + (∇ S σ) i ] + F e,i (19)
where F e,i corresponds to the electric field force and is the source term in the momentum equation calculated from Equation ( 7) that is obtained from the discrete form of the electric field density

E = - ρ i m i ∑ j m i ρ i 2 + m j ρ j 2 φ ij ∇W ij (20)
with φ ij being the inter-particle average of the electric potential. The electric body force links the electric field equations to the momentum balance.

To close the system, the energy equation is discretized as

DT Dt i = ∑ j m j ρ j λ i + λ j x ij x 2 ij ∇ i W ij • T i -T j (21)
where T i is the ith particle temperature and λ i is the ith particle thermal conductivity. It must be noted that any type of phase change in the system is neglected here due to the low temperature gradient compared to the fluid properties. If not negligible, the effect of phase change on the energy equation can be represented in the form of additional terms.

Pressure Calculation and Temporal Discretization

One way to calculate the pressure term in the momentum balance is by enforcing the incompressiblity condition. It is therefore called the Incompressible SPH (ISPH) [24] method, which consists of the following steps. First, the Helmholtz-Hodge decomposition of the momentum balance is obtained. This decomposition divides the momentum balance into a part including pressure (divergence-free) and a pressure-free part (curl-free) shown in Equation (19). From the curl-free contribution, one can estimate an intermediate velocity and density, which represents a predictor step. From a Pressure-Poisson Equation (PPE), one can compute the pressure term. Pressure is applied to correct the estimated velocity (corrector step). This implies that in this method an estimation of velocity is projected on a divergence-free space.In SPH, this projection method was introduced by Cummins and Rudman [24]. The discrete equations used in the current study are reviewed in the following. Having that the curl-fee part of the acceleration in the intermediate step (indicated with asterisk * ), a a, * in the momentum balance can be defined as

a i, * = ∑ j m j ρ i ρ j µ i + µ j r ij r 2 ij ∇ i W ij • v i -v j + F e,i + F g,i + n i ρ i (σ i χ i n i + (∇ S σ) i ) (22)
Next, we obtain the intermediate velocity using an explicit Euler scheme

v i, * = v t i + a i, * ∆t (23)
where ∆t is the time step. Again, sing the intermediate velocity and an explicit Euler scheme, we obtain the new particle position

x i, * = x t i + v i, * ∆t (24)
This position is then used to calculate the intermediate density according to Equation (18). Given the calculated intermediate density, the particles can be put back to their prior position x t i . The pressure can therefore be calculated according to the PPE

∇ • ∇p ρ * = ∇ • v * ∆t (25)
More details about projection procedure can be found in reference [19] .

In the present model, the temporal discretization is obtained by the fractional step method [36]. The details of the application of this method for velocity and position advancement using the so-called intermediate velocity calculation is elaborated at references [19,26]. In the current study, two time-step integration schemes are used, namely a predictor-corrector scheme to solve the momentum equation and an explicit Euler scheme to sequentially couple the mass fraction to the velocity in the predictor step. A stability criterion, based on the CFL condition is needed following the minimum time scales between convection and diffusion schemes. In the present study, we consider fixed boundary conditions for all conducted cases. For each test mentioned further we adopt the method of Morris et al. [37], where the boundary conditions are projected to fixed boundary particles to fulfill either Neumann or Dirichlet boundary conditions [25,38,39].

Results and Discussion

This section presents the results of the SPH model implemented in the 2D studies. First, in order to verify the SPH model in the pressure prediction, a Young-Laplace problem is solved with SPH model in a 2D simulation domain. Second, a linear thermal profile is imposed on a 2D simulation domain to verify the SPH model against analytical solution of the Marangoni force in the static case followed by a dynamic case to validate the 2D thermoacapillary bubble rising and comparison of the SPH results with direct numerical simulation results. Third, the effects of EHD on a single droplet immersed in continues phase are calculated with SPH model and compared with analytical results. This includes the flow orientation inside and outside of the droplet, the bubble deformation and the velocity compared with analytical results. In the last section, the SPH model for two-phase flow subject to coupled EHD and thermocapillary forces is presented and the evolution of a 2D droplet is predicted for multiple fluid properties.

Validation of the interfacial forces begins with the investigation of a static pressure jump. Initially, a quiescent system with a velocity field equal to zero is assumed. Since all time derivatives are zero, the mass is conserved. Young-Laplace law relates the droplet curvature and the pressure gradient at the interface. Based on this law, the pressure gradient at the interface of two phases would be equal to the product of the mean curvature of the interface and the surface tension such as

p d -p b = σ 1 R 1 + 1 R 2 (26)
where R 1 and R 2 are the radii of the curvature of a curved surface. Note that in case of a circle

R 1 = R 2 = R resulting in p d -p b = σ R .
The pressure jump at phase boundary follows the Young-Laplace equation, which describes the relationship between the pressure, the surface tension and the radius of the droplet. The schematic of the pressure inside and outside of the bubble are shown in Figure 2-left. To investigate the accuracy of the numerical results compared to analytical results when droplet radius R are taken 0.25 [m], 0.3125 [m], 0.375 [m] and 0.625 [m]. The L/R ratio is set to 8 in these simulations so that the same confinement effect will be applied for all. It is shown in Figure 2-right that as expected, a linear relationship between ∆p and 1 R is observed. The slope of each straight ling indicates the surface tension and are equal to σ = 0.1 and σ = 0.2. This indicates that the simulation results obtained from SPH multiphase model confirm well with Young-Laplace law. According to Equation (26), surface tension tends to minimize the surface, whereas the pressure difference tends to increase the surface curvature. When considering free-surface flows and based on the geometry, the fluid-fluid interface is flat while in the problems stated here, one phase is fully surrounded by the other and is under full tension from all directions, leading to droplet circular shape (when no external forces are applied).

To study the effect of particle resolution, the simulation setup is designed with the domain size of 4 [cm], sufficiently low confinement effects with L/R = 4 and a surface tension of σ = 0.00837 [N/m]. The following relative error norms are defined.

L 2 (p Exact ) ≡ 1 p Exact ∑ i (p SPH,i -p Exact ) 2 ∑ i L ∞ (p Exact ) ≡ 1 p Exact max i (|p SPH,i -p Exact |) (27)
Table 1 shows that when the particle resolutions increase, both L 2 and L ∞ norms decrease. We can conclude that the pressure gradient inside the bubble converges to the analytical solution ∆p = 0.837 [Pa].

σ = σ re f + σ T (T -T re f ) (28)
where σ re f and T re f are the reference values of the surface tension and the temperature, respectively. The negative surface tension coefficient σ T = -0.002 2. The gravity and the heat dissipation are neglected, therefore the only acting force is the Marangoni force for which the sampled profile along the horizontal line is depicted and one can observe that the higher the particle spacing, the wider the range of volumetric Marangoni force at the interface. These results are in agreement with results reported in [19,40]. In thermocapillary droplet motion, having that surface tension depends on temperature, assuming hot wall temperature T H and cold wall temperature T C on parallel boundaries (T H > T C ) leads to the introduction of a surface tension gradient along the interface. Thermocapillarity consists of applying a temperature gradient along an interface to induce a surface tension gradient. If ∇T and ∂σ denote the temperature and the surface gradient operator, respectively, the thermocapillary tangential stress writes as

∇ s σ = ∂σ ∂T ∇ s T (29)
The surface tension gradients (i.e., Marangoni effect) can be used to control the dynamics of the bubble. To this end, a square box (5.76 [mm] Marangoni stress that is the tangential component of the surface tension gradient acts in opposition of the surface motion and hence results in less flexibility of the interface and droplet mobility restriction if no additional force is applied. As time passes, the timeevolution of the non-dimensional droplet velocity which, and after some oscillations, meets the velocity associated direct numerical simulation Ma and Bothe [41], hence the initial unsteady leading to a steady droplet motion.

Hereafter, the EHD solver is validated by setting up a similar test to first part of the Section 4. The direction of these streamlines is determined as the mutual relationship between the electrical conductivity and electrical permitivity ratios of the phases. They are defined as S = d / b and R = ς d /ς b for the electrical permitivity and electrical conductivity ratios, respectively. Note that the subscript d and b refer to the droplet and bulk fluid properties.

The recirculation zones of the fluid inside and outside of the droplet are theoretically predicted by Taylor et al. [42] as depicted in Figure 5 with Pole-to-Equator and Equator-to-Pole flow directions. A specific numerical simulation for the case S D4 with R < S and S D3 for R > S from Table 3 are compared with Taylor's results in Figure 6.

A qualitative agreement of the flow orientation inside and outside of the droplet is found between the theoretical prediction and the SPH capability in correctly capturing [START_REF] Rahmat | Numerical simulation of the electrohydrodynamic effects on bubble rising using the sph method[END_REF] Results

the direction of the recirculation zones. For R > S, the streamlines start from the equator towards the pole, while R > S situation, the opposite direction is observed. Investigation of circular droplet deformation given small deformations subject to electric field is presented here using two theories from literature. Taylor [42] estimates the droplet deformation D T as

D T = 9 f dT E 2 o d R 8(2 + R) 2 σ ( 30
)
where f dT is the discriminating function evaluated as

f dT = R 2 + 1 -3.5S + 1.5R (31)
For the same problem, Feng [44] suggests the following relation

D F = f dF E 2 o b R 3(1 + R) 2 Sσ ( 32
)
where f dF is estimated from

f dF = R 2 + 1 -3S + R (33)
In Equations ( 30) and ( 32), R is the initial droplet radius before its deformation and E o is the electric field magnitude in the vertical direction deduced from the electric potential difference E o = (ϕ +ϕ -)/h, with h being the height of the domain. Numerically, the droplet deformation parameter D can be defined based on the droplet's deviation from the circular shape,

D N = A -B A + B ( 34
)
where A and B are the elliptic droplet diameters at the steady-state condition, parallel and perpendicular to the direction of the external electric field, respectively. When D N = 0, the droplet is at its initial circular shape. On the other hand, more deviation from zero indicates more deformation from its initial shape. Positive and negative values of D N refer to deformation in the direction and perpendicular to the electric field direction, respectively.

For our study, we consider a circular droplet with R = 0. 4 represents the obtained along with the theoretical values of Feng and Taylor, and the droplet deformed shape at the steady-state. These setups are selected from the simulations proposed by Shadloo et al. [43] and show good agreement with this study. Feng [44] has also proposed an analytical solution for the droplet velocity subject to DC electric field, assuming that the droplet deformation is negligible, i.e., the final shape of the droplet is assumed circular. The velocities inside and outside the droplet can be theoretically calculated as where r is the radial position, v θ and v r are the tangential velocity and the radial velocity, respectively. The characteristic velocity U can be evaluated as

U = R -S 2S (1 + R) 2 d E 2 o R µ d + µ b ( 39
)
Our numerical results are in agreement with the analytical solutions for both the velocity θ = 0 • and θ = 45 • at which one of the velocity components is maximized as shown in Figure 8. As can be observed in this figure, theoretical velocity profile for tangential and radial components are shown with dashed-lines and solid lines in the given order. The simulation data for tangential and radial velocity components are respectively depicted with filled circles and unfilled circles. The simulation input parameters correspond to S D4 as shown in Table 3 with (R = 0.5 [m], E o = 1 [V/m]). When the droplet deformation is sufficiently large such that the terminal shape can no longer be assumed as circular in two dimension, the conformity of the numerical and analytical results tends to reduce. Based on the Equations ( 35) and ( 37) the radial velocity needs to be zero at the droplet interface (r = R). However, because of the small deformation of the droplet (r = R) and causes small deviation between the numerical and analytical results. Based on the equations of tangential components (36) and (38), having (sin 2θ = 0) leading to v θ = 0 [m/s] and maximum values of v r that is in agreement with the observations in Figure 8. In this section the simultaneous effects of thermocapillarity and electrohydrodynamics on a single suspended droplet are studied. When a multiphase system is solely subject to [START_REF] Safdari Shadloo | Special topic on turbulent and multiphase flows[END_REF] Results electric field, the motion of the dispersed phase is forced by the electric field and damped by the viscous forces. The instabilities caused by EHD-driven flows will occur when the viscous drag force is much smaller than the electric force. Hence, the droplet losses its balance and starts to migrate and/or deform. The electric force is more sensitive to the droplet size than the viscous counterpart, augmenting more instability for larger L R ratios,

where R is the initial droplet radius and L is the size of the domain.

As mentioned before, the Marangoni effect roots into the tangential component of the surface tension gradient which if aligned with the tangential component of the electric force, drives the flow to a more unstable configuration. Subsequently, in the next section we will discuss the scaling parameters that lead to observation of both convective Marangoni effects and electrohydrodynamic effects. To capture electrohydrodynamics driven droplet deformations at the equilibrium state, the electric and the hydrodynamic time scales should match properly based on the geometry of the system namely the droplet radius and the channel width. and ς are electric permittivity and electric conductivity. This accumulation of the bulk free charges at the phase boundaries will further result in creation of dipoles on the droplet. The electric field acting on these induced free charges, will generate shear stress at the fluid-fluid interface that can be balanced by shear viscous stress and, if thermal gradient is applied, the tangential component of the surface tension gradient. Another aspect to take into account regarding the instability of an EHD-TC multiphase system, is the difference between prolate and oblate droplet deformation. Assuming a vertical potential difference across the domain, and given that weakly conductive droplets embedded in the weakly conductive fluids tend to flip such that the dipole aligns with the direction of the electric field, one can expect that the oblate deformations are more prone to potential instabilities as their dipole is in the opposite direction of the electric field. 5. The time evolution of the droplet subject to thermocapillary flow and EHD forces is illustrated in Figure 9. 5.

Based on the results shown in Figure 9, a combination of R and S ratios and the direction of thermal and electric potential gradients affect the deformation and migration of the droplet.Thermocapillary induced motion due to surface tension gradient can be decomposed into two components.The perpendicular temperature gradient component to the fluid layer that generates Benard-Marangoni circulations inside the droplet. The tangential temperature gradient component generate surface tension gradient along the surface of the fluid and induce surface flow towards the regions with higher surface tension. In our simulation, we combine these two components to obtain surface force. Furthermore, when a droplet is suspended in an imposed flow field generated by the EHD forces, which itself depends on the R and S relation as explained in the EHD section; both surface force and the imposed EHD force are responsible for the deformation and migration of the droplet. For the case 1 (top row) and case 2 (middle row), vertical and horizontal electric field are applied, respectively, while the thermal gradient is kept vertical in both cases. All the other parameters are kept the same. In case 1, where the electric field is oriented the same as thermal gradient, in the vertical direction, it is observed that the droplet forms a prolate shape. Because the viscosity of the droplet (µ d = 0.12 [Pa•s]) is chosen close to continues phase viscosity (µ b = 0.24 [Pa•s]), similar to those in the thermocapillary induced motion cases, the resistance due to the presence of the droplet is relatively low. However, the non-uniform distribution of the electric charges on the surface of the droplet, generates a shear force (from equator-to-pole since R > S) which, in addition to the Marangoni stress caused by the variation of surface tension on the droplet surface, deforms the droplet. This deformation modifies the effective viscosity of the droplet compared to its initial state with respect to the continues field. Because the surface tension coefficient is negative σ T = -0.002 [N/m], the droplet tends to move in the opposite direction of the thermal gradient, that is from top to bottom. But the effective viscosity modification, retards this migration. The top side of the droplet, closer to hot wall, has lower surface tension, and it is as a result more prone to strong deformation caused by both internal EHD-induced circulations and surface-force-induced circulations. In the case 2, where horizontal electric field in applied with higher electric potential set to be at the right side of the setup, our results show a symmetric elongation along the electric field direction while the droplet is moving downwards. Finally, case 3 (bottom row) shows the evolution of the droplet towards a break-up when a vertical electric field is applied taking R < S. As can be observed, the droplet starts to deform from its center while keeping symmetrical oblate structure thought the break-up process. Comparison between case 3 in coupled EHD-TC cases and case 4 in Table 4 where both systems are characterised by R < S, we observe a different orientation for the droplet due to the presence of Marangoni force. It is important to take into account the difference between the scale of effectiveness in EHD phenomena and Marangoni flows due to thermal gradient. The former has a really small time scale and therefor affects the system much faster than the Marangoni flows. This highlights the reason for-which numerical simulation of coupled physics require many trials as time-scale and length scale are not in the same order for each physics involved.

Conclusions

•

After independently validating the effect of thermocapillary (TC) and electrohydrodynamic (EHD) forces on two-phase flows, the coupled EHD-TC effects on a single droplet migration is studied.

•

The two-phase system is analysed first at each isolated condition through which multiple field parameters such as system's velocity and pressure are validated.

•

The behaviour of the droplet subject to EHD forces and TC phenomena is characterized using the deformation values for different material properties. • By choosing the appropriate ratios of electrical permitivity, electrical conductivity, thermal gradient and electric potential gradient, it is now possible to use simulations to predict and control the migration and deformation of a droplet using SPH method.

•

The results show strong agreement with previous literature. Simulation results for a single droplet subject to direct current electric field and Marangoni flow induced by thermal gradient in the absence of gravity reveal an insight to the evolution of the droplet migration and deformation.

•

Results show that the main parameters that influence the droplet topology include the R and S corresponding to electrical permitivitty and electrical conductivity, respectively as well as the orientation of the applied electric field with respect to thermal gradient.

•

For R > S and thermal and electric gradient applied both in the vertical direction, a prolate deformation of the droplet is observed.

•

Having R > S and an electric field perpendicular to the thermal gradient results in symmetric and oblate deformation of the droplet while displacing the droplet in the opposing direction of the thermal gradient.

•

The investigation of the droplet break-up process under the coupled effect of the thermal and electric field demonstrates that in the case of R < S with both electric and thermal gradients applied vertically, the droplet's topology evolves towards a break-up state. • Nevertheless, the different time-scales of each phenomenon motivates further researches to provide better understanding of EHD-TC instabilities where gravity effects are neglected.

•

Regarding the limitations of this research work, it is worth mentioning that despite having some similarities with topics such as EHD Heat transfer enhancement, it is yet to be found sufficiently related numerical or experimental studies with the close physics, dimensions and parameters to compare with. The approach that we decided to choose is to validate our model extensively, using available related recent studies, and propose our solution to a novel coupled multi-physics problem which we believe is valid based on the validated components. We expect that new experimental and numerical studies emerge to discuss this configuration following to the pioneer appearance of this work.

•

For future works, we believe that uncertainty quantification and reliability analysis of this model can be subject of future directions. In that scope, a meta-model is used to couple the mechanical and probabilistic models.

•

In this study, given that the foundations of the mathematical model, in its general form, to explain the contribution of each force and the general mechanism affecting the deformation of the droplet are established in this work, the parametric and sensitivity analysis can be subject of further studies.

R Electrical conductivity of the droplet to bulk ratio E o

Electric field vector magnitude (V/m) set to be E o = Chapter 5

Conclusion and Future work 5.1 Synthesis

In this project a 2D numerical simulation of the coupled problem concerning multi-phase fluid system has been studied. Several test cases have been carried out to find more about the deformation and rising of bubbles subject to an electric field. The influence of some parameters has also been studied through benchmark problems. This chapter is dedicated to provide an initial conclusion to the project and some recommendations for the continuation of this research.

• Paradoxically, as computational power increases, we become aware of the finer scale effects and limitations that our physical model poses. Furthermore, we are more sensitive to the spread of error and uncertainty. Therefore, despite the extended limitations on computing power, in the near future we will continue to be hampered in our efforts to understand more complex coupled multiscale multiphysics phenomena.

For this purpose, mathematical and computational modeling will remain the main supporting technologies that need to be further developed and exploited.

• Multi-phase fluid flow modelling

The numerical modelling and simulation of multi-phase fluid flows especially when coupled with electric or magnetic field or with interaction with solid boundary has numerous industrial applications including energy, aerospace and bio-medical sectors.

To solve the governed system of equations many grid-based and mesh-less methods have been proposed.

• SPH method SPH is a mesh-less Lagrangian method that solves the N-S equations on freely moving particles of the incompressible fluid system. In the multi-phase systems, each phase is described by its own set of material particles, thus, interface problems are inherently simpler when it comes to distinguishing the phases, large deformations, and complex geometries.

• Boundary conditions

As for the square or rectangular computational domain of the problem geometry, each side is attributed to have a "tag". Assigning one of the pre-defined values to the corresponding field variable of each tag could be interpreted as no-slip, free-slip and "to be computed". The periodic boundary condition is also as option for some cases.

• Time and space discretization

The computational domain consists of a 2D array of initially arranged particles in Cartesian grid for the background fluid (oil) and circular particles for the droplet or bubble ( water or gas). The time step used in the computation is relatively small ( order of E-5 for some cases) which is set to avoid instabilities. For a writing interval varying between 200 to 1000, one should consider a significant computation time (order of few days) considering that the code is developed in serial. The latter emphasizes the interest to paralleling and distributing the computation.

2. 4 1 Introduction 1 . 1

 4111 Three distinguished flow regimes can be observed by increasing temperature gradient (left) steady flow, (middle) oscillatory flow, and (right) turbulent flow reproduced sketch from: https://iss.jaxa.jp/ . . . . . . . . . . . . . . xiv List of Figures 2.5 A molecule at the fluid-fluid interface (here, magnified an colored in red) is exposed to different tensions compared to a molecule inside the bulk. . . 2.6 Comparison of ISPH and VOF method [159]. . . . . . . . . . . . . . . . 2.7 Dispersed phases or particle laden flows in (left) air bubble inside water and (right) air-oil-in-water. reproduced from [97]. . . . . . . . . . . . . . . . 2.8 Dispersed Flow classification from low to high volume fraction. reproduced from [97] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.9 Kernel function W (x, h) and its first derivative ∂ ∂ x (W (xx ′ , h)) according to equation 2.32-2.33. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2.10 Interface height change versus time (a) and spatiotemporal evolution of interface b(i-iv) in combined EHD-Marangoni instabilities. Applied voltage, ψ up = 5V, initial film thickness, h 0 = 26nm, relative thermal conductivity, k r = 2 and modified Marangoni number, M = 0.22. Nondimensional times, T =(i) 1.2 × 10 4 ,(ii) 1.5 × 10 4 , (iii) 1.7 × 10 4 , (iv) 2.0 × 10 4 and (v) 7.0 × 10 4 .λ = 5.34µm [96]. . . . . . . . . . . . . . . . . . . . . . . . . 2.11 Profiles of local Nusselt numbers on the hot wall at different electric Rayleigh number[151] . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3.1 The components of the total stress tensor and the surface tension gradient at the interface . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .Chapter Multi-phase modeling for multi-physics problemsMultiphase flows are used to refer to systems formed by contacting more than one partially miscible or immiscible fluids where the components possess different chemical and/or mechanical properties. Multiphase fluids are observed in almost every process. The energy generation in gas turbines and diesel engines or the production of nano particles with specific properties are examples of industrial applications of multiphase flows.

Fig. 1 . 1

 11 Fig. 1.1 Examples of multiphase flows

)

  With ρ density, p pressure, η dynamic viscosity, D Dt the material time derivative operator1, ∇.⃗ u the divergence of the velocity field, κ the thermal conductivity and c p the specific heat capacity at constant pressure. Herein, ⃗ F body = ρ⃗ g is the body force due to gravitational acceleration ⃗ g and Π capillary the capillary stress tensor and ⃗ F e is the volumetric electric force. The first three terms in the momentum balance are common between single-phase and multi-phase flows, while the forth terms accounts for the fluid-fluid interface and the last term includes the external electric force contribution.

Faraday's law

  of induction states that electric fields exist even without the presence of charges. Changing magnetic flux through a surface induces an electromotive force (EMF) in any boundary path of that surface. A changing magnetic field induces a circulating electric field. The voltage accumulated around a closed circuit is proportional to the time rate of change of the magnetic flux it encloses.This is in comparison with electric fields, which are generated by charges (sources) and are not rotation-free. It is a so-called vortex field and is therefore source-free, i.e. ∇.E = 0. The integral and differential form of this law are shown in 2.3.

Fig. 2 . 1

 21 Fig. 2.1 (left) Velocity and (right) temperature distributions in a buoyancy-driven cavity for Ra=10000 reproduced from [58].

Fig. 2 . 3

 23 Fig. 2.3 Surface Marangoni flow induced due to surface temperature gradient from hot to cold region.The flow propagates inside due to drag and reverse flow forms to hold mass conservation (a ). Droplet on a solid surface with temperature gradient leans toward the cold region due to higher surface tension at the inception of movement and Marangoni circulation forms inside it (b ). In the case of channel flows, the carrier liquid obeys a bulk flow from hot to cold region due to surface tension forces and the confined droplet moves to the opposite direction (toward the hot region) to conserve the mass (c ). Schematic reproduced form[START_REF] Karbalaei | Thermocapillarity in microfluidics-a review[END_REF]

Figure 2 . 3 illustrates

 23 the Marangoni flow generated in response to the tangential component of the surface tension gradient that leads to formation of circulation zones inside the dispersed phase.

  Navier-Stokes equations describe the motion of fluid by which a relationship between velocity field and stress field (shear stress, normal stress, and external forces) is established to track fluids' temporal and spatial dynamics.

  surface tension is considered to balance the normal and tangential components of the stress boundary condition at the fluid-fluid interface. The hydrodynamic normal stress balance at the interface reads for every n unit normal vector to the fluids surface n.T.nn. T .n = σ (∇.n) (2.23)

Fig. 2 . 4

 24 Fig. 2.4 Three distinguished flow regimes can be observed by increasing temperature gradient (left) steady flow, (middle) oscillatory flow, and (right) turbulent flow reproduced sketch from: https://iss.jaxa.jp/

Fig. 2 . 5 A

 25 Fig. 2.5 A molecule at the fluid-fluid interface (here, magnified an colored in red) is exposed to different tensions compared to a molecule inside the bulk.

Fig. 2 . 7

 27 Fig. 2.7 Dispersed phases or particle laden flows in (left) air bubble inside water and (right) air-oilin-water. reproduced from [97].

Fig. 2 . 8

 28 Fig. 2.8 Dispersed Flow classification from low to high volume fraction. reproduced from[START_REF] Pachpute | Multi-Phase Flow and its CFD Modeling[END_REF]

  )where ⃗ x ′ is the position of any point in the volume Ω. The Dirac-Delta function δ (|⃗ x -⃗ x ′ |) is defined to be unity for ⃗ x = ⃗ x ′ and zero otherwise. Thereafter, an interpolation function, also known as kernel (smoothing) function, W (h, |⃗ x -⃗ x ′ |), is defined to replace the Dirac-Delta function where h is the radius of the support volume, so-called smoothing length and its value depends on type of the kernel function, dimension of the problem as well as type of the problem.

Fig. 2 . 9

 29 Fig. 2.9 Kernel function W (x, h) and its first derivative ∂ ∂ x (W (xx ′ , h)) according to equation 2.32-2.33.

Fig. 2 .

 2 Fig. 2.10 Interface height change versus time (a) and spatiotemporal evolution of interface b(i-iv) in combined EHD-Marangoni instabilities. Applied voltage, ψ up = 5V, initial film thickness, h 0 = 26nm, relative thermal conductivity, k r = 2 and modified Marangoni number, M = 0.22. Nondimensional times, T =(i) 1.2 × 10 4 ,(ii) 1.5 × 10 4 , (iii) 1.7 × 10 4 , (iv) 2.0 × 10 4 and (v) 7.0 × 10 4 .λ = 5.34µm [96].

  .41) where e has a small positive value representing the thermal expansion coefficient. Assuming zero nanoparticle mass flux, They argued that stability of the system mainly depends on the thermal Rayleigh number R a = ρ o β ∆T gL 3 µα , the Lewis number L e = α D B , the modified diffusivity ratio N a = D T ∆T (φ o D B T 2 ) and the nanoparticle Rayleigh number R N = (ρ p -ρ o )φ o gL 3 µα

Fig. 2 .

 2 Fig. 2.11Profiles of local Nusselt numbers on the hot wall at different electric Rayleigh number[START_REF] Wang | Spectral element method for numerical simulation of ethd enhanced heat transfer in an enclosure with uniform and sinusoidal temperature boundary conditions[END_REF]

3 ) 1

 31 where v = (u, v, w) the fluid flow velocity, density ρ, pressure p, the gravitational acceleration g and the total viscous stress tensor Π viscous = µ[∇v + (∇v) T ], and the Boussinesq approximation Here, eq. 3.1 is the continuity equation, which describes the mass conservation of the fluid flow under consideration. Equation 3.2 is is about maintaining momentum balance.The convective force term (ρv -∇v) on the right side of eq. 3.2 describes the convective transport within the fluid. The viscous force term (η∇ 2 v) presents frictional effects. The pressure term indicates the pressure the forces acting on the fluid again. One can define f = ( f x, f y, f z) as external forces (source terms). For instance, the gravitational force g or 1Density variation is only important in the buoyancy term, and can be neglected in the rest of the equation are assumed.

  boundary conditions are required. The mainly used boundary conditions are the Dirichlet and Neumann conditions. A Dirichlet boundary condition will specify the value of the variable on the edge. For the Neumann boundary condition defines the normal gradient of a variable. Periodic boundary conditions are used to combine an infinitely large area which is regenerated in time.

  By defining the interfacial tension σ and the stress exerted by one fluid on the other by n.T and n. T , respectively, we can derive the total stress balance at the fluid-fluid interface.n.Tn. T = σ n(∇.n) -∇.σ(3.18) Interface between two phases with their respective normal unit vector n,-n, stress components uAn element of the interface in subject to a normal stress T nn and a tangential component T tt , the gradient of the surface tension ∇σ .

Fig. 3 . 1

 31 Fig. 3.1 The components of the total stress tensor and the surface tension gradient at the interface

  using Smoothed Particle Hydrodynamics in Process Engineering (SIPER ©2016) which is developed by a team of engineers and researchers, mainly Manuel Hopp-Hirschler, Philip Kunz and Christian Zander, at the Institute of Chemical Process Engineering (ICVT), Stuttgart University. It is an in-house code for the description of multiphase flows. For this purpose, the ICVT relies on a particle-based, grid-free method called "Smoothed Particle Hydrodynamics" (SPH). Due to the Lagrangian property, this approach is suitable on the one hand for the dynamic description of interfaces and for problems with geometrically changing boundary conditions.This program is written in C language, organized with doxygen2, uses XML-PlugIn 3 components and combines them into a deployable format under eclipse IDE 4. Massage Passing Interface (MPI) has been implemented to parallelize the code by dividing the 2Doxygen is a free software used as documentation generator and static analysis tool for software source trees.

  approach provides a clear definition for the volume fraction of each fluid, i.e. C A i = C i and C B i = 1 -C i define the volume fraction corresponding to the fluid A and fluid B, respectively, such that ∑ n C n i = 1 for all n phases, here n = 2.

Fig. 2 Fig. 2 -

 22 Fig.2the data is represented for 60 × 60, 100 × 100 and 140 × 140 grids. To study numerical convergence, a droplet with the radius of r = 0.01 [m] is situated at the center of a square domain, i.e. x o /r = y o /r = 2, with the side lengths of x/r = y/r = 4 (see Fig.1). While the Dirichlet (no-slip) boundary condition is set for the velocity at all four boundaries, namely, BC-X1, BC-X2, BC-Y1, and BC-Y2, the Neumann boundary condition is applied for the pressure field. As for the hydrodynamics properties, we keep both viscosity and density ratios equal to unity such that ρ 1 = ρ 2 = 1000 [kg/m 3 ] and µ 1 = µ 2 = 0.1 [Pa s] and set the surface tension to γ = 0.01 [N/m], given neither electrical nor gravitational force.This problem, known as Young-Laplace problem, has an analytical solution which is ∆p = p i -p o = γ /r = 1. As can be seen in Fig.2-Left, the pressure oscillations decrease by increasing the resolution from the coarser to the finest and the results, converging towards the analytical solution. It is noted that the relative error is less than 1% for the intermediate particle resolution. Therefore, we chose the 100 × 100 resolution for our simulations as it provides accurate results with reasonable computation cost. Similar simulations, with different surface tension coefficient are tested, while putting r = 0.5 [m] and keeping x/r and y/r ratios constant to validate the accuracy of the used method. As can be seen in Fig.2-Right, the pressure jump increases by an increment in the surface tension confirming the capability of the method

Fig. 3 .

 3 Fig. 3. Time evolution of bubble rising problem with the density ratio of ρ 1 /ρ 2 = 0.1 and viscosity ratio of µ 1 /µ 2 = 0.1 at different dimensionless times t * = 0.5, t * = 18.5, t * = 37 and t * = 64.8. Here the dimensionless time is defined as t * = t √ (g/D) and Reynolds and Bond numbers are Re = 35 and Bo = 10, respectively. The left half of each snapshot shows the velocity streamlines in black and the droplet interface in red, while the right half shows the velocity magnitude's 10 highest levels' contour in the range of [0, 4.5] [m/s]. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

ρ 2 /ρ 1 =

 21 10 with ρ 1 = 100 [kg/m 3 ], µ 2 /µ 1 = 10 with µ 1 = 1 [Pa s], and surface tension coefficient is γ = 24.5 [N/m].

Fig. 4 .

 4 Fig. 4. The average normalized central position of the droplet, ȳ (top) and its average vertical migration velocity, v (bottom) as a function of dimensionless time, t * .

  computational domain. The initial zero velocity are assigned to both fluids and wall particles. Density ratio and viscosity ratio are set to unity with values of ρ 1 = ρ 2 = 1000 [kg/m 3 ] and µ 1 = µ 2 = 1 [Pa s]. The surface tension coefficient γ = 1 [N/m]. The velocity and pressure boundary conditions are exactly the same as those imposed in Section 5.1. The electrical boundary conditions are Dirichlet (ϕ = cte.) and Neumann boundary (∇ϕ • ⃗ n = 0) conditions for horizontal (i.e. BC -Y 1 and BC -Y 2) and vertical walls (i.e. BC -X 1 and BC -X 2), respectively, where ⃗ n is normal direction to the given boundary.

Fig. 5 .

 5 Fig. 5. Comparison of bubble rising and break up due to gravitational and hydrodynamic forces using experiment [39], SPH method and Volume of

Fig. 6 .

 6 Fig. 6. Deformation of a suspended droplet in response to an external vertical electric field in the steady-state simulation for the cases with (left) S = 0.5, R = 0.05 and (right) S = 0.5, R = 3.

Fig. 7 .

 7 Fig. 7. Comparison of deformation for all cases. permittivity ratio, conductivity ratio and surface tension coefficient of each simulation, S, R, γ , respectively, mentioned below or above the corresponding case.

Fig. 8 .

 8 Fig. 8. Schematic of the Couette flow test case.

Fig. 9 .

 9 Fig. 9. Numerical deformation for twelve cases (Left), a close-up look at the rate of deformation at the steady state (Right). Cases without electric field are denoted by a black + sign. The pair number on the legend box, corresponds to the electrical permittivity and electrical conductivity (S, R), respectively.

Fig. 10 .

 10 Fig. 10. Bubble interface at t * = 1. The pair number above each case corresponds to the electrical permittivity and electrical conductivity (S, R), respectively. If zero, the electric field is not applied. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Fig. 11 .

 11 Fig. 11. Deformation of a suspended droplet in Couette flow with S = 10 and R = 0.2 subject to electric field at t * = 0, t * = 0.4, t * = 0.8, t * = 1.2 and t * = 1.6, respectively, from top to bottom where dimensionless time is defined by t * = tU o /r. The velocity streamlines (in blue), the electric field vectors (in black) and the droplet interface (in red) are shown at five moments. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Figure 1 .

 1 Figure 1. Kernel functionW(x, h) and its first derivative ∂ ∂x (W(x -x , h)) according to Equations (13) and(14).

Figure 2 .

 2 Figure 2. (left) Schematic of the pressure inside (p d ) and pressure outside (p b ) the bubble used to calculate theoretical pressure according to Young-Laplace law, (right) The validation of Young-Laplace law. Straight lines indicate analytical results and scattered points indicate numerical results from SPH method.

Table 1 .

 1 Compariosn between the exact analytical results and SPH results in the pressure computation in the case of 2D droplet simulation with 60, 100, and 140 particles per direction. validate the Marangoni force separately and investigate convergence of the numerical method. First, the static capillary stress tangential to the interface under thermocapillary effect at three grid resolutions is studied. A linear thermal profile is imposed on a two-layer square (5.76 [mm] × 5.76 [mm]) domain where droplet is characterised by ρ d = 250 [kg/m 3 ], λ d = 96 × 10 -6 [m 2 /s] and the background fluid by ρ b = 500 [kg/m 3 ], λ b = 48 × 10 -6 [m 2 /s] for density and thermal conductivity, respectively. The surface tension evolves based on the Equation (28) showing a linear relationship between the surface tension and the temperature,

  [N/mK] implies the decrease of the surface tension σ with respect to the temperature. The linear thermal gradient 200 [K/m] is imposed from bottom wall T C towards the upper wall T H . According to the dependence of the surface tension to the temperature with σ T = -0.002 [N/mK], the Marangoni force acts vertically on the interface. Note that the lateral walls are adiabatic while the top and the bottom walls are subject to T H = [291.15] K and T C = 290 [K].Figure3shows the profile of interfacial Marangoni force along a horizontal line at the center of the simulation domain, at three different particle spacing. The values form SPH model are the normal component of the Marangoni force perpendicular to the fluid-fluid interface which are calculated based on the Continuum Surface Force (CSF) method and as a consequence are volumetric force calculated per particle volume. Lower particle spacing correspond to sharper interface. If an infinite resolution was possible numerically, the Marangoni force profile would tend to Dirac function where its magnitude is equal to Marangoni force. From Figure3one can also observe that the magnitude of the Marangoni force increases with particle resolution augmentation. Since the σ T = -0.002 [N/mK], the interfacial Marangoni force direction is downwards. To compare the magnitude numerical results with exact solution (σ T × ∇T = 0.4 [N/m 2 ] ). For the numerical case, magnitude of the Marangoni force predicted by SPH method is computed by integrating the Marangoni force part of the surface tension (not the CSF part). The comparative table of the results in shown in Table

Figure 3 .

 3 Figure 3. Profile of the volumentric Marangoni force along the normal direction to the interface using a resolution of Lo = 180 [µm], Lo = 90 [µm] and Lo = 45 [µm].

Figure 4 .

 4 Figure 4. Comparison between droplet migration velocity in this work and Ma and Bothe [41].

Figure 5 .Figure 6 .

 56 Figure 5. Schematics representation of two types of induced flows based on the Taylor's theory: (a) R < S and (b) R > S. Reproduced from Shadloo et al. [43].

Figure 7 .

 7 Figure 7. Pressure distribution across the left side of the droplet and the vector velocity contours inside and outside of the droplet on the right side. The red to blue colors indicate high to low pressure region.

  5 [m] placed at the center of a squared domain of L = 4 [m]. The domain contains 240 particles per direction. The droplet properties are ρ d = 1000 [kg/m 3 ] and µ d = 1 [Pa•s] while the bulk fluid has identical density and viscosity. All four boundaries are set to no-slip velocity and Neumann pressure boundary conditions except for the top wall where a Dirichlet pressure boundary condition is used. An electrical field is imposed on the system by E = 1 [V/m] resulting in unique electric force value directed towards the bottom wall.

Figure 8 .

 8 Figure 8. Components of the velocity profile at (a) θ = 0 • and (b) θ = 45 • compared to the formulae given by Feng [44].

  The EHD time scale can be characterised by Maxwell-Wagner polarization time t = d + 2 b ς d + 2ς b . It is worth reminding that electric charge accumulation at the fluid-fluid interface happens when each phase has a different charge relaxation time τ = ς where

  The coupled EHD-TC problem setup consists of a square domain of size L = 0.04 [m] with L R = 8 where R is the radius of the centralised droplet at x o = 0.02 [m] and y o = 0.02 [m] with 120 × 120 particle resolution in x and y directions, respectively. An electric potential of φ = 0.04 [V] is imposed on the top boundary while other boundaries are set to φ = 0 [V]. A linear thermal profile is imposed on the top and bottom walls with respective temperatures of T top = 300 [K] and T bottom = 290 [K] while the reference temperature and the initial temperature are also kept at 290 [K] throughout the simulation. The velocity boundary conditions are set to be free slip at the lateral walls, and no slip at the top and bottom walls. Also, pressure boundary conditions are set to be Dirichlet with a constant value at top wall and Neumann for the other three boundaries (∇p • n = 0) where n is the normal direction to the given boundary. The density and viscosity of are chosen to be ρ d = 250 [kg/m 3 ], µ d = 0.12 [Pa•s] for the bubble and ρ b = 500 [kg/m 3 ], µ b = 0.24 [Pa•s] for the bulk phase. Both phases are set to have stationary conditions at initial time step. The fluid electrical properties of the three cases are given in Table

Table 5 .tFigure 9 .

 59 Figure 9. Time evolution of the droplet interface for EHD-TC coupled simulations; case 1 (top row), case 2 (middle row), case 3 (bottom row) with vertical temperature gradient set to 10 [K]. The physical simulation time t = 0 [s] at the very first frame on the left, up to t = 4 [s], with a time increment of 1 per frame. For the simulation parameters see Table 5.

Table 2 .

 2

1 Gauss's law for static electric fields expressed in two forms.

Table 2 . 2

 22 Maxwell-Ampèr's law for electric field changing with time producing a magnetic field expressed in two forms.

Table 2 . 3

 23 Faraday's Law of Induction expressed in two forms.

Table 2 .4

 2 Gaussian law for electric fields expressed in two forms.

Table 2 . 5

 25 Summary of approaches in multiphase modelling.

	Issues	In contrast to Eulerian-	Lagrangian approach, Eule-	rian models does not mini-	mize numerical diffusion in	dispersed-phase fields e.g.	volume fraction and mean	velocity and are not well	posed for successive grid	refinements.	Computationally expensive	for complex geometries	with large deformations and	strong interface instabili-	ties.				behavior of phases arises	from interaction of the par-	ticles alone and is irrelevant	to the position of interface
	Applicability	In problems with great	number of particles within	each cell.							Free-surface flows and	droplets migration and de-	formation.						capable of dealing with any	interface deformations and	fragmentations.
	Modelling	Transport equations are	needed and computed for	each cell							Mass conservation is used	to track the interface and	momentum conservation is	shared by the fluids while	the volume fraction of each	phase in the control volume	is tracked throughout the	whole domain	fluid into a set of discrete	elements, referred as parti-	cles and is based on integral	representation of quantities	and spatial derivatives.	Application of Navier-	Stokes equations for all	phases.
	Approach Comments	Eulerian Each phase is modeled sep-	arately using fluid particle	groups.							V-O-F Surface-tracking technique	is applied to a fixed Eule-	rian mesh.						Lagrangian Each phase is modeled sep-	arately and the interphase	is inherently distinguished.	Eulerian-Velocity, pressure and tem-	Lagrangian perature fields for each	phase can be accurately pre-	dicted as continuous func-	tions using the Eulerian ap-	proach, while the droplet	properties is evaluated us-	ing Lagrangian trajectories.

Table 2 .

 2 Comparison between SPH results of the integrated Marangoni force and analytical solution.

	Particle Spacing [µm]	Numerical Result [N]	Analytical Solution [N]	Relative Error Percentage
	45 90 180	0.394421 0.394411 0.394425	0.4 0.4 0.4	0.013947 0.013971 0.013975

Table 3 .

 3 Simulation parameters for EHD cases using small deformation theory.

Table 3

 3

summarizes different EHD test cases and Table

Table 4 .

 4 Comparison

	Case	D N	D T	D F	Bubble Shape
	S D1	0.077	0.065	0.061

of numerical (N) and theoretical results (Taylor et al.

[42]

: Equation (

30

) and Feng

[44]

: Equation (32) of the droplet deformation (D) for different combinations of conductivity and permittivities.

4The material time derivative is a directional time derivative for a fixed point

4.2 Results: Numerical Simulation of Multi-phase Electro-Hydrodynamics Systems

Energies 2022, 15, 2576. https://doi.org/10.3390/en15072576 https://www.mdpi.com/journal/energies

4.3 Results: Coupled EHD-TC Droplet Migration

Acknowledgements

Acknowledgments

This work has benefited from the financial support of TOTAL SA within the framework of a call for projects with exploratory projects coordinated by TOTAL, France SA and the CNRS, France. The authors also acknowledge the access to French HPC resources provided by the French regional computing center of Normandy CRIANN (2017002). The first author acknowledges the financial supports provided by Normandy region, France. The second author acknowledges the support provided by Alexander von Humboldt Foundation, Germany through the project FRA-1204799-HFST-E for the experienced researcher. The last author gratefully acknowledge the financial support provided by the Scientific and Technological Research Council of Turkey (TUBITAK) for grant number 112M721 and by the Natural Sciences and Engineering Research Council of Canada (NSERC) for grant number 05862, the Canada Research Chair program and the Canada Foundation for Innovation.

This study has been carried out at the COmplexe de Recherche Interprofessionnel en Aérothermochimie (CORIA UMR-6614) and National Institute of Applied Sciences (for his precious advice, inspiration and scientific discussions.

Funding: This research received no external funding.

Acknowledgments:

The authors acknowledge the access to French HPC resources provided by the French regional computing center of Normandy CRIANN (2022006). The first author acknowledges the financial supports provided by Normandy region council. The third author acknowledge the support provided by Alexander von Humboldt Foundation through the project FRA-1204799-HFST-E for the experienced researcher.

Conflicts of Interest:

The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Nomenclature

Future works

• Due to the fat that SPH is a fully Lagrangian method and its efficiency in modeling interfaces, it has been greatly developed for multiphse and free-surface problems.

To take advantage of the Eulerian methods in the continuous phase where complex deformations are non-existent, coupling SPH with mesh-based methods such as Finite Volume Methods is presented for simulation of incompressible interfacial flows with large density difference. In case of coupled multiphase-multi-physics problems, hybrid numerical approaches, such as coupled Level Set-SPH and FEM-SPH, can be used to improve the efficiency of the method and accuracy of the results.

• New approaches to simulate multiphase -multiphysics problems are conquering the research summits. Recent advances in HPC and Machine Learning application in thermal systems has been investigated in special edition with topics such as thermal parameters prediction and flow modelling [START_REF] Safdari Shadloo | High-performance computing and machine learning applied in thermal systems analysis preface[END_REF]. Liu et al. [START_REF] Liu | The effect of alcohol-gasoline fuel blends on the engines' performances and emissions[END_REF] proposed an artificial neural network (ANN) modeling to predict the exhaust emissions and engine performance by investigating different alcohol-gasoline fuel blends. Such new approaches on multi-phase problems could be applied even for complicated and nonlinear cases without the need for expensive, complex, and time-consuming experimental studies. In an investigation to predict the viscosity of biodiesel blends in temperatures ranging from 268.15 to 373.15K Zheng et al. [START_REF] Zheng | Prediction of viscosity of biodiesel blends using various artificial model and comparison with empirical correlations[END_REF] showed excellent agreement between ANN predictions and experimental data. The ANN models have also been used in the pressure drop prediction in two-phase long horizontal pipes, Shadloo et al. [START_REF] Shadloo | Estimation of Pressure Drop of Two-Phase Flow in Horizontal Long Pipes Using Artificial Neural Networks[END_REF], Haghighi et al. [START_REF] Haghighi | Using Committee Neural Network for Prediction of Pressure Drop in Two-Phase Microchannels[END_REF] with ANN models presenting good agreement with experimental data. Artificial intelligence applied to heat transfer exchangers [START_REF] Ghalandari | Applications of intelligent methods in various types of heat exchangers: a review[END_REF] is also a new promising field of research in multi-phase multi-physics phenomena.