Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method

: This paper concerns the study of coupled effects of electrohydrodynamic (EHD) and thermocapillary (TC) on the dynamic behaviour of a single liquid droplet. An incompressible Smoothed Particle Hydrodynamic (ISPH) multiphase model is used to simulate EHD-TC driven ﬂows. The complex hydrodynamic interactions are modeled using the continuum surface force (CSF) method, in which the gradient of the interfacial tension and the Marangoni forces are calculated with an approximated error or 0.014% in the calculation of Marangoni force compared to the analytical solutions which is a signiﬁcant improvement in comparison with previous SPH simulation studies, under the assumption that the thermocapillarity generates sufﬁciently large stress to allow droplet migration, while the electrohydrodynamic phenomena inﬂuences the droplet morphology depending on the electrical and thermal ratios of the droplet and the ambient ﬂuid. This study shows that, when applying a vertical electric ﬁeld and thermal gradient, the droplet starts to stretch horizontally towards a break-up condition at a high rate of electrical permitivity. The combined effect of thermal gradient and electric ﬁeld tends to push further the droplet towards the break-up regime. When the thermal gradient and the electric ﬁeld vector are orthogonal, results show that the droplet deformation would take place more slowly and the Marangoni forces cause the droplet to migrate, while the stretching in the direction of the electric ﬁeld is not seen to be as strong as in the ﬁrst case.


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 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 where p is the pressure, 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 To calculate the electric force F e , we apply the electrohydrodynamics theory for leaky dielectric fluids [21,22].The electric force is obtained from 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 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, 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]

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 : 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 where the Ψ wall is the value of the wall and the sign function is defined as 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.

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 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(x − x , h)dx = 1 and contracts to Dirac-delta function δ, so that lim h→0 W(x − x , h) = δ(x ).
Moreover, a suitable kernel function should have a compact support; i.e., for every This last condition ensures the numerical efficiency in approximating physical variables such as velocity, density or temperature.The discrete form of Equation (11) reads Here, we use the 4-th order Wendland C2 kernel function [31] given by where q is the dimensionless smoothing length, q = x − x h .The normalization constant C w at each dimension d is 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 ( 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 ∂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.

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 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 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 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 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 Next, we obtain the intermediate velocity using an explicit Euler scheme where ∆t is the time step.Again, sing the intermediate velocity and an explicit Euler scheme, we obtain the new particle position 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
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 where R 1 and R 2 are the radii of the curvature of a curved surface.Note that in case of a circle 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.
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].Next, we 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.  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, where Figure 3 shows 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 Figure 3 one 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 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 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 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.Figure 7 depicts the velocity profile on the right-half and the pressure distribution on the left-half.Note that the recirculation zones of velocity streamlines inside (in blue) and outside (in red) of the droplet are clearly visible.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 where f dT is the discriminating function evaluated as For the same problem, Feng [44] suggests the following relation where f dF is estimated from 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, 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.  3 summarizes different EHD test cases and Table 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 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 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 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.The EHD time scale can be characterised by Maxwell-Wagner polarization It is worth reminding that electric charge accumulation at the fluid-fluid interface happens when each phase has a different charge relaxation time τ = ς where 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.
Table 5. Simulation parameters of the coupled EHD-TC cases where vertical (V) and horizontal (H) electric fields are imposed.5.

Case Electric Field Direction
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.

•
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 =

Figure 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.
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 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[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].
× 5.76 [mm]) is discretized using 32 particles in each direction.The droplet is initially placed at the center of the domain and has a radius R = 1.4 [mm].The velocity boundary conditions are set to be free slip at the lateral walls, and no-slip at the top and bottom walls.Neumann Pressure boundary conditions are used on all walls except for the left wall, where a Dirichlet Pressure boundary condition is used due to bootstrap condition of the Pressure Poisson Equation.The temperature is fixed at 290 [K] at the bottom wall and linearly increases to 291.15 [K] at the top wall.Assuming that both droplet and the ambient fluid are initially at the stationary state with ρ d = 250 [kg/m 3 ] and µ d = 0.012 [Pa•s] for the droplet and ρ b = 500 [kg/m 3 ] and µ b = 0.024 [Pa•s]) for the background fluid, we consider that the heat conductivity of the droplet (1.2 × 10 −6 [W/mK]) is half of the background fluid.With constant surface tension σ re f = 0.01 [N/m] and the rate of the change of the surface tension with temperature σ T = 0.002 [N/mK] in Equation (28), Figure 4 shows the time evolution of the droplet migration velocity subject to Marangoni force compared to results obtained by Ma and Bothe [41].The dimensionless velocity U * = U U c where the characteristic velocity U c = σ T |∇T|R µ b and the dimensionless time t * = t •

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

Figure 5 .Figure 6 .
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 .Table 3 .
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.Table 3. Simulation parameters for EHD cases using small deformation theory.Case d [F/m] ς d [S/m] S R σ [N/m] S D1 0.3 40 0.5 2 0.01 S D2 00.5 40 0.5 2 0.01 S D3 0.5 150 0.5 3 0.01 S D4 0.5 1 0.5 0.05 0.01 S D5 3 10 5 0.5 0.03 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.Table
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

Figure 9 .
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 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.

Table 2 .
Comparison between SPH results of the integrated Marangoni force and analytical solution.
1 unless stated otherwise S Electrical permittivity of the droplet to bulk ratio φ