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

Fatemeh Almasi, Manuel Hopp-Hirschler, Abdellah Hadjadj, Ulrich Nieken, Mostafa Safdari Shadloo

► To cite this version:

Fatemeh Almasi, Manuel Hopp-Hirschler, Abdellah Hadjadj, Ulrich Nieken, Mostafa Safdari Shadloo. Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method. Energies, 2022, 15 (7), pp.2576. 10.3390/en15072576. [hal-03824997](https://hal.science/hal-03824997)

HAL Id: hal-03824997

https://hal.science/hal-03824997

Submitted on 21 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method

Fatemeh Almasi, Manuel Hopp-Hirschler, Abdellah Hadjadj, Ulrich Nieken, Mostafa Safdari Shadloo

► To cite this version:

Fatemeh Almasi, Manuel Hopp-Hirschler, Abdellah Hadjadj, Ulrich Nieken, Mostafa Safdari Shadloo. Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method. Energies, 2022, 15 (7), pp.2576. 10.3390/en15072576. [hal-03824997](https://hal.science/hal-03824997v1)

HAL Id: hal-03824997

https://hal.science/hal-03824997v1

Submitted on 21 May 2024

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
Article

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

Fatemeh Almasi<fnref n="1" />*, Manuel Hopp-Hirschler<fnref n="2" />, Abdellah Hadjadj<fnref n="1" />, Ulrich Nieken<fnref n="2" /> and Mostafa Safdari Shadloo<fnref n="1" />,<fnref n="2" />

**Abstract:** 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 flows. 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 significant improvement in comparison with previous SPH simulation studies, under the assumption that the thermocapillarity generates sufficiently large stress to allow droplet migration, while the electrohydrodynamic phenomena influences the droplet morphology depending on the electrical and thermal ratios of the droplet and the ambient fluid. This study shows that, when applying a vertical electric field and thermal gradient, the droplet starts to stretch horizontally towards a break-up condition at a high rate of electrical permittivity. The combined effect of thermal gradient and electric field tends to push further the droplet towards the break-up regime. When the thermal gradient and the electric field 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 field is not seen to be as strong as in the first case.

**Keywords:** electrohydrodynamics; thermocapillary; multiphase fluid flows

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Almasi, F.; Hopp-Hirschler, M.; Hadjadj, A.; Nieken, U.; Safdari Shadloo, M. Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method. *Energies* **2022**, *15*, 2576. https://doi.org/10.3390/en15072576

Academic Editors: Gabriela Huminic and Fabrizio Ascione

Received: 13 January 2022

Accepted: 24 March 2022

Published: 1 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. 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

<footnote n="1">¹ National Institute for Applied Sciences, INSA Rouen Normandie, Clean Combustion Laboratory, CORIA UMR 6614 CNRS, 76000 Rouen, France; abdellah.hadjadj@insa-rouen.fr (A.H.); msshadloo@coria.fr (M.S.S.)</footnote>

<footnote n="2">² Institute of Chemical Process Engineering, University of Stuttgart, 70199 Stuttgart, Germany; manuel.hirschler@icvt.uni-stuttgart.de (M.H.-H.); ulrich.nieken@icvt.uni-stuttgart.de (U.N.)</footnote>

<footnote>* Correspondence: almasif@coria.fr</footnote>
Article

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

Fatemeh Almasi<fnref n="1" />*, Manuel Hopp-Hirschler<fnref n="2" />, Abdellah Hadjadj<fnref n="1" />, Ulrich Nieken<fnref n="2" /> and Mostafa Safdari Shadloo<fnref n="1" />,<fnref n="2" />

**Abstract:** 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 flows. 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 significant improvement in comparison with previous SPH simulation studies, under the assumption that the thermocapillarity generates sufficiently large stress to allow droplet migration, while the electrohydrodynamic phenomena influences the droplet morphology depending on the electrical and thermal ratios of the droplet and the ambient fluid. This study shows that, when applying a vertical electric field and thermal gradient, the droplet starts to stretch horizontally towards a break-up condition at a high rate of electrical permittivity. The combined effect of thermal gradient and electric field tends to push further the droplet towards the break-up regime. When the thermal gradient and the electric field 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 field is not seen to be as strong as in the first case.

**Keywords:** electrohydrodynamics; thermocapillary; multiphase fluid flows

<figure><img src="image_2.png" /></figure>

check for
updates

**Citation:** Almasi, F.; Hopp-Hirschler, M.; Hadjadj, A.; Nieken, U.; Safdari Shadloo, M. Coupled Electrohydrodynamic and Thermocapillary Instability of Multi-Phase Flows Using an Incompressible Smoothed Particle Hydrodynamics Method. *Energies* **2022**, *15*, 2576. https://doi.org/10.3390/en15072576

Academic Editors: Gabriela Huminic and Fabrizio Ascione

Received: 13 January 2022

Accepted: 24 March 2022

Published: 1 April 2022

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

<figure><img src="image_3.png" /></figure>

**Copyright:** © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/).

## 1. 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

<footnote n="1">¹ National Institute for Applied Sciences, INSA Rouen Normandie, Clean Combustion Laboratory, CORIA UMR 6614 CNRS, 76000 Rouen, France; abdellah.hadjadj@insa-rouen.fr (A.H.); msshadloo@coria.fr (M.S.S.)</footnote>

<footnote n="2">² Institute of Chemical Process Engineering, University of Stuttgart, 70199 Stuttgart, Germany; manuel.hirschler@icvt.uni-stuttgart.de (M.H.-H.); ulrich.nieken@icvt.uni-stuttgart.de (U.N.)</footnote>

<footnote>* Correspondence: almasif@coria.fr</footnote>
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 can-
not 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 Eule-
rian 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 Hy-
drodynamics (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 inhomogen-
ities 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 sys-
tem with thermal gradient and electric field are first presented, before introducing the
SPH method and the associated space and time discretization schemes. Afterwards, nu-
merical 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.

## 2. Governing Equations

Mathematical formulations of governing equations of the coupled EHD-TC problem consist of a set of continuity, momentum and energy balance equations.
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 can-
not 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 Eule-
rian 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 Hy-
drodynamics (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 inhomogen-
ities 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 sys-
tem with thermal gradient and electric field are first presented, before introducing the
SPH method and the associated space and time discretization schemes. Afterwards, nu-
merical 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.

## 2. Governing Equations

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

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

$$ \frac{D\rho}{Dt} + \rho\nabla \cdot \mathbf{v} = 0 \qquad (1) $$

where $\rho$, $\mathbf{v}$ and $D\rho/Dt$ represent the density, the averaged fluid velocity and the material derivative operator, respectively.

For incompressible flows, the continuity equation reduces to $\nabla \cdot \mathbf{v} = 0$.

## 2.2. 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 permittivity and electrical conductivity can be written as

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \nabla \cdot \bar{\Pi}_v + \nabla \cdot \bar{\Pi}_c + \mathbf{F}_g + \mathbf{F}_e \qquad (2) $$

where $p$ is the pressure, $\bar{\Pi}_v = \mu[\nabla\mathbf{v} + (\nabla\mathbf{v})^T]$ is the viscous stress tensor with the dynamic viscosity $\mu$ and transpose operator $T$, $\bar{\Pi}_c$ is the capillary stress tensor, $\mathbf{F}_g$ is the body force due to gravity, and $\mathbf{F}_e$ is the electric force. Note that the sharp interface limit of the capillary stress tensor is defined as $\nabla \cdot \bar{\Pi}_c = (\sigma\chi\mathbf{n} + \nabla_S\sigma)\delta$ where $\sigma$ is the surface tension, $\chi$ is the curvature, $\nabla_S$ is the gradient of the surface tension in tangential direction with respect to the surface, and $\mathbf{n}$ is the unit normal vector of the interface. The Dirac-delta function $\delta$ 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

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \mu \nabla^2 \mathbf{v} + (\sigma\chi\mathbf{n} + \nabla_S\sigma)\delta + \rho g + \mathbf{F}_e \qquad (3) $$

To calculate the electric force $\mathbf{F}_e$, we apply the electrohydrodynamics theory for leaky dielectric fluids [21,22]. The electric force is obtained from

$$ \mathbf{F}_e = \nabla \cdot \bar{T} \qquad (4) $$

where $\bar{T}$ is the Maxwell stress tensor. For weak electric currents, the magnetic field is negligible because the electric field is assumed to be irrotational ($\nabla \times \mathbf{E} = 0$). The Maxwell stress tensor $\bar{T}$ reads

$$ \bar{T} = \mathbf{D}\mathbf{E} - \frac{1}{2}(\mathbf{D} \cdot \mathbf{E})\bar{I} \qquad (5) $$

where $\mathbf{E}$ is an external electric field, $\mathbf{D} = \epsilon\mathbf{E}$ is the dielectric displacement vector, and $\epsilon$ is the electrical permittivity. $\bar{I}$ denotes the identity tensor. Using Gauss's law,

$$ \nabla \cdot \mathbf{D} = q^v \qquad (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 $\mathbf{F}_e$ can be obtained as [23]

$$ \mathbf{F}_e = q^v \mathbf{E} - \frac{1}{2} \mathbf{E} \cdot \nabla \mathbf{E} \qquad (7) $$

## 2.3. Conservation of Energy

In the energy balance equation, we only include the heat transfer induced by the temperature gradient ($\nabla T$) and we neglect the effect of viscous heating and other source
## 2.1. Continuity Equation

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

$$ \frac{D\rho}{Dt} + \rho\nabla \cdot \mathbf{v} = 0 \qquad (1) $$

where $\rho$, $\mathbf{v}$ and $D\rho/Dt$ represent the density, the averaged fluid velocity and the material derivative operator, respectively.

For incompressible flows, the continuity equation reduces to $\nabla \cdot \mathbf{v} = 0$.

## 2.2. 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 permittivity and electrical conductivity can be written as

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \nabla \cdot \bar{\Pi}_v + \nabla \cdot \bar{\Pi}_c + \mathbf{F}_g + \mathbf{F}_e \qquad (2) $$

where $p$ is the pressure, $\bar{\Pi}_v = \mu[\nabla\mathbf{v} + (\nabla\mathbf{v})^T]$ is the viscous stress tensor with the dynamic viscosity $\mu$ and transpose operator $T$, $\bar{\Pi}_c$ is the capillary stress tensor, $\mathbf{F}_g$ is the body force due to gravity, and $\mathbf{F}_e$ is the electric force. Note that the sharp interface limit of the capillary stress tensor is defined as $\nabla \cdot \bar{\Pi}_c = (\sigma\chi\mathbf{n} + \nabla_S\sigma)\delta$ where $\sigma$ is the surface tension, $\chi$ is the curvature, $\nabla_S$ is the gradient of the surface tension in tangential direction with respect to the surface, and $\mathbf{n}$ is the unit normal vector of the interface. The Dirac-delta function $\delta$ 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

$$ \rho \frac{D\mathbf{v}}{Dt} = -\nabla p + \mu \nabla^2 \mathbf{v} + (\sigma\chi\mathbf{n} + \nabla_S\sigma)\delta + \rho g + \mathbf{F}_e \qquad (3) $$

To calculate the electric force $\mathbf{F}_e$, we apply the electrohydrodynamics theory for leaky dielectric fluids [21,22]. The electric force is obtained from

$$ \mathbf{F}_e = \nabla \cdot \bar{T} \qquad (4) $$

where $\bar{T}$ is the Maxwell stress tensor. For weak electric currents, the magnetic field is negligible because the electric field is assumed to be irrotational ($\nabla \times \mathbf{E} = 0$). The Maxwell stress tensor $\bar{T}$ reads

$$ \bar{T} = \mathbf{D}\mathbf{E} - \frac{1}{2}(\mathbf{D} \cdot \mathbf{E})\bar{I} \qquad (5) $$

where $\mathbf{E}$ is an external electric field, $\mathbf{D} = \epsilon\mathbf{E}$ is the dielectric displacement vector, and $\epsilon$ is the electrical permittivity. $\bar{I}$ denotes the identity tensor. Using Gauss's law,

$$ \nabla \cdot \mathbf{D} = q^v \qquad (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 $\mathbf{F}_e$ can be obtained as [23]

$$ \mathbf{F}_e = q^v \mathbf{E} - \frac{1}{2} \mathbf{E} \cdot \nabla \epsilon \qquad (7) $$

## 2.3. Conservation of Energy

In the energy balance equation, we only include the heat transfer induced by the temperature gradient ($\nabla 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 :

$$ \frac{DT}{Dt} = \kappa \nabla^2 T \qquad (8) $$

where $κ$ is the thermal diffusivity.

## 2.4. 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 $\Psi'_i$ and a particle $\Psi_i$ (such as velocity, pressure, temperature, color) can be defied such that

$$ \Psi'_{i} = 2\Psi_{wall} - \delta_{wall}\Psi_{i} \qquad (9) $$

where the $\Psi_{wall}$ is the value of the wall and the sign function is defined as

$$ \delta_{wall} = \begin{cases} +1 & \text{for Dirichlet boundary conditions.} \\ -1 & \text{for Neumann boundary conditions with } \bar{\Psi}_{wall} = 0 \end{cases} \qquad (10) $$

More details can be found in [25,26]. It is more suitable for straight wall conditions (i.e., our current work) compared to curved walls.

# 3. Smoothed Particle Hydrodynamics (SPH) Method

## 3.1. 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) = \int f(x')W(x - x', h)dx' \qquad (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 $\Omega$, such as $\int_{\Omega} W(x - x', h)dx = 1$ and contracts to Dirac-delta function $\delta$, so that $\lim_{h \to 0} W(x - x', h) = \delta(x')$.

Moreover, a suitable kernel function should have a compact support; i.e., for every $k \in \mathbb{R}^+$ if $|x - x'| > kh$ then $W(x - x', h) = 0$.
terms. However, the heat transfer from one phase to the other is accounted for. The resulting energy equation can be written as :

$$ \frac{DT}{Dt} = \kappa \nabla^2 T \qquad (8) $$

where $\kappa$ is the thermal diffusivity.

## 2.4. 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 $\Psi'_i$ and a particle $\Psi_i$ (such as velocity, pressure, temperature, color) can be defied such that

$$ \Psi'_{i} = 2\Psi_{wall} - \delta_{wall}\Psi_{i} \qquad (9) $$

where the $\Psi_{wall}$ is the value of the wall and the sign function is defined as

$$ \delta_{wall} = \begin{cases} +1 & \text{for Dirichlet boundary conditions.} \\ -1 & \text{for Neumann boundary conditions with } \bar{\Psi}_{wall} = 0 \end{cases} \qquad (10) $$

More details can be found in [25,26]. It is more suitable for straight wall conditions (i.e., our current work) compared to curved walls.

# 3. Smoothed Particle Hydrodynamics (SPH) Method

## 3.1. 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) = \int f(x')W(x - x', h)dx' \qquad (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 $\Omega$, such as $\int_{\Omega} W(x - x', h)dx = 1$ and contracts to Dirac-delta function $\delta$, so that $\lim_{h \to 0} W(x - x', h) = \delta(x')$.

Moreover, a suitable kernel function should have a compact support; i.e., for every $k \in \mathbb{R}^+$ 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) = \sum_j f(x_j) W(x_i - x_j, h) V_j \quad (12)$$

Here, we use the 4-th order Wendland C2 kernel function [31] given by

$$W(x - x', h) = \frac{C_w}{h^d} \begin{cases} (1 - \frac{q}{2})^4 (1 + 2q) & q \in [0, 2]. \\ 0 & \text{otherwise.} \end{cases} \quad (13)$$

$$W'(x - x', h) = \frac{C_w}{h^d} \begin{cases} (1 - \frac{q}{2})^3 (-5q) & q \in [0, 2]. \\ 0 & \text{otherwise.} \end{cases} \quad (14)$$

where $q$ is the dimensionless smoothing length, $q = \frac{x - x'}{h}$. The normalization constant $C_w$ at each dimension $d$ is

$$C_w = \begin{cases} \frac{3}{4} & d = 1 \\ \frac{7}{4\pi} & d = 2 \\ \frac{21}{16\pi} & d = 3 \end{cases} \quad (15)$$

Figure 1 illustrates the kernel function $W$ and its derivative with smoothing radius $h = 0.75$.

<figure><img src="image_4.png" /><figcaption>Figure 1. Kernel function $W(x, h)$ and its first derivative $\frac{\partial}{\partial x}(W(x - x', h))$ according to Equations (13) and (14).</figcaption></figure>

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].
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) = \sum_j f(x_j) W(x_i - x_j, h) V_j \quad (12)$$

Here, we use the 4-th order Wendland C2 kernel function [31] given by

$$W(x - x', h) = \frac{C_w}{h^d} \begin{cases} (1 - \frac{q}{2})^4 (1 + 2q) & q \in [0, 2]. \\ 0 & \text{otherwise.} \end{cases} \quad (13)$$

$$W'(x - x', h) = \frac{C_w}{h^d} \begin{cases} (1 - \frac{q}{2})^3 (-5q) & q \in [0, 2]. \\ 0 & \text{otherwise.} \end{cases} \quad (14)$$

where $q$ is the dimensionless smoothing length, $q = \frac{x - x'}{h}$. The normalization constant $C_w$ at each dimension $d$ is

$$C_w = \begin{cases} \frac{3}{4} & d = 1 \\ \frac{7}{4\pi} & d = 2 \\ \frac{21}{16\pi} & d = 3 \end{cases} \quad (15)$$

Figure 1 illustrates the kernel function $W$ and its derivative with smoothing radius $h = 0.75$.

<figure><img src="image_4.png" /><figcaption>Figure 1. Kernel function $W(x, h)$ and its first derivative $\frac{\partial}{\partial x}(W(x - x', h))$ according to Equations (13) and (14).</figcaption></figure>

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 guarantees that the gradient of a constant field ($f_i = f_j$) is exactly zero.

$$ \nabla f(x)_i = \sum_j \frac{m_j}{\rho_j} (f_j - f_i) \nabla W_{ij} \qquad (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 $\nabla W_{ij} = \frac{x_i - x_j}{|x_i - x_j|} \frac{\partial W}{\partial 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.

$$ \left( \frac{1}{\rho} \nabla f(x) \right)_i = \sum_j m_j \left( \frac{f_j + f_i}{\rho_i \cdot \rho_j} \right) \nabla W_{ij} \qquad (17) $$

## 3.2. 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

$$ \rho_i = m_i \sum_j W_{ij} \qquad (18) $$

In this formulation, $\rho_i$ is calculated directly from $\frac{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

$$ \begin{aligned} \frac{D\mathbf{v}}{Dt} = & -\sum_j \frac{m_j}{\rho_i \cdot \rho_j} (p_i + p_j) \nabla W_{ij} \\ & + \sum_j \frac{m_j}{\rho_i \cdot \rho_j} \left( \mu_i + \mu_j \right) \frac{x_i - x_j}{|x_i - x_j|^2} \nabla f(x)_i \cdot (\mathbf{v}_i - \mathbf{v}_j) \\ & + \mathbf{F}_{g,i} + \delta[\sigma_i \chi_i \mathbf{n}_i + (\nabla_S \sigma)_i] + \mathbf{F}_{e,i} \end{aligned} \qquad (19) $$

where $\mathbf{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

$$ \mathbf{E} = -\frac{\rho_i}{m_i} \sum_j \left[ \left( \frac{m_i}{\rho_i} \right)^2 + \left( \frac{m_j}{\rho_j} \right)^2 \right] \phi_{ij} \nabla W_{ij} \qquad (20) $$

with $\phi_{ij}$ being the inter-particle average of the electric potential. The electric body force links the electric field equations to the momentum balance.
The first formulation, called the negative formulation, Monaghan [34], is often used for the divergence free velocity since it guarantees that the gradient of a constant field ($f_i = f_j$) is exactly zero.

$$ \nabla f(x)_i = \sum_j \frac{m_j}{\rho_j} (f_j - f_i) \nabla W_{ij} \quad (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 $\nabla W_{ij} = \frac{x_i - x_j}{|x_i - x_j|} \frac{\partial W}{\partial 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.

$$ \left( \frac{1}{\rho} \nabla f(x) \right)_i = \sum_j m_j \left( \frac{f_j + f_i}{\rho_i \cdot \rho_j} \right) \nabla W_{ij} \quad (17) $$

## 3.2. 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

$$ \rho_i = m_i \sum_j W_{ij} \quad (18) $$

In this formulation, $\rho_i$ is calculated directly from $\frac{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

$$ \begin{aligned} \frac{D\mathbf{v}}{Dt} = & -\sum_j \frac{m_j}{\rho_i \cdot \rho_j} (p_i + p_j) \nabla W_{ij} \\ & + \sum_j \frac{m_j}{\rho_i \cdot \rho_j} \left( \mu_i + \mu_j \right) \frac{x_i - x_j}{|x_i - x_j|^2} \nabla f(x)_i \cdot (\mathbf{v}_i - \mathbf{v}_j) \\ & + \mathbf{F}_{g,i} + \delta[\sigma_i \chi_i \mathbf{n}_i + (\nabla_S \sigma)_i] + \mathbf{F}_{e,i} \end{aligned} \quad (19) $$

where $\mathbf{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

$$ \mathbf{E} = -\frac{\rho_i}{m_i} \sum_j \left[ \left( \frac{m_i}{\rho_i} \right)^2 + \left( \frac{m_j}{\rho_j} \right)^2 \right] \phi_{ij} \nabla W_{ij} \quad (20) $$

with $\phi_{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

$$ \left(\frac{DT}{Dt}\right)_i = \sum_j \frac{m_j}{\rho_j} (\lambda_i + \lambda_j) \frac{\mathbf{x}_{ij}}{x_{ij}^2} \nabla_i W_{ij} \cdot (T_i - T_j) \quad (21) $$

where $T_i$ is the ith particle temperature and $\lambda_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.

### 3.3. Pressure Calculation and Temporal Discretization

One way to calculate the pressure term in the momentum balance is by enforcing the incompressibility 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-free part of the acceleration in the intermediate step (indicated with asterisk *), $\vec{a}_{i,*}$ in the momentum balance can be defined as

$$ \mathbf{a}_{i,*} = \sum_j \frac{m_j}{\rho_i \rho_j} (\mu_i + \mu_j) \frac{\mathbf{r}_{ij}}{r_{ij}^2} \nabla_i W_{ij} \cdot (\mathbf{v}_i - \mathbf{v}_j) + \mathbf{F}_{e,i} + \mathbf{F}_{g,i} \\ + \frac{\mathbf{n}_i}{\rho_i} (\sigma_i \chi_i \mathbf{n}_i + (\nabla_S \sigma)_i) \quad (22) $$

Next, we obtain the intermediate velocity using an explicit Euler scheme

$$ \mathbf{v}_{i,*} = \mathbf{v}_i^t + \mathbf{a}_{i,*} \Delta t \quad (23) $$

where $\Delta t$ is the time step. Again, sing the intermediate velocity and an explicit Euler scheme, we obtain the new particle position

$$ \mathbf{x}_{i,*} = \mathbf{x}_i^t + \mathbf{v}_{i,*} \Delta t \quad (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 $\mathbf{x}_i^t$. The pressure can therefore be calculated according to the PPE

$$ \nabla \cdot \left( \frac{\nabla p}{\rho_*} \right) = \frac{\nabla \cdot \mathbf{v}_*}{\Delta t} \quad (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
To close the system, the energy equation is discretized as

$$ \left(\frac{DT}{Dt}\right)_i = \sum_j \frac{m_j}{\rho_j} (\lambda_i + \lambda_j) \frac{\mathbf{x}_{ij}}{x_{ij}^2} \nabla_i W_{ij} \cdot (T_i - T_j) \quad (21) $$

where $T_i$ is the ith particle temperature and $\lambda_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.

### 3.3. Pressure Calculation and Temporal Discretization

One way to calculate the pressure term in the momentum balance is by enforcing the incompressibility 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-free part of the acceleration in the intermediate step (indicated with asterisk *), $\vec{a}_{i,*}$ in the momentum balance can be defined as

$$ \mathbf{a}_{i,*} = \sum_j \frac{m_j}{\rho_i \rho_j} (\mu_i + \mu_j) \frac{\mathbf{r}_{ij}}{r_{ij}^2} \nabla_i W_{ij} \cdot (\mathbf{v}_i - \mathbf{v}_j) + \mathbf{F}_{e,i} + \mathbf{F}_{g,i} \\ + \frac{\mathbf{n}_i}{\rho_i} (\sigma_i \chi_i \mathbf{n}_i + (\nabla_S \sigma)_i) \quad (22) $$

Next, we obtain the intermediate velocity using an explicit Euler scheme

$$ \mathbf{v}_{i,*} = \mathbf{v}_i^t + \mathbf{a}_{i,*} \Delta t \quad (23) $$

where $\Delta t$ is the time step. Again, sing the intermediate velocity and an explicit Euler scheme, we obtain the new particle position

$$ \mathbf{x}_{i,*} = \mathbf{x}_i^t + \mathbf{v}_{i,*} \Delta t \quad (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 $\mathbf{x}_i^t$. The pressure can therefore be calculated according to the PPE

$$ \nabla \cdot \left( \frac{\nabla p}{\rho_*} \right) = \frac{\nabla \cdot \mathbf{v}_*}{\Delta t} \quad (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].

## 4. 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 = \sigma \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \qquad (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 = \frac{\sigma}{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 $\Delta p$ and $\frac{1}{R}$ is observed. The slope of each straight line indicates the surface tension and are equal to $\sigma = 0.1$ and $\sigma = 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 $\sigma = 0.00837$ [N/m]. The following relative error norms are defined.

$$L_2(p_{\text{Exact}}) \equiv \frac{1}{p_{\text{Exact}}} \sqrt{\frac{\sum_i (p_{\text{SPH},i} - p_{\text{Exact}})^2}{\sum_i}} \qquad (27) \\ L_\infty(p_{\text{Exact}}) \equiv \frac{1}{p_{\text{Exact}}} \max_i (|p_{\text{SPH},i} - p_{\text{Exact}}|)$$

Table 1 shows that when the particle resolutions increase, both $L_2$ and $L_\infty$ norms decrease. We can conclude that the pressure gradient inside the bubble converges to the analytical solution $\Delta p = 0.837$ [Pa].
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].

## 4. 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 = \sigma \left( \frac{1}{R_1} + \frac{1}{R_2} \right) \qquad (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 = \frac{\sigma}{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 $\Delta p$ and $\frac{1}{R}$ is observed. The slope of each straight line indicates the surface tension and are equal to $\sigma = 0.1$ and $\sigma = 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 $\sigma = 0.00837$ [N/m]. The following relative error norms are defined.

$$L_2(p_{\text{Exact}}) \equiv \frac{1}{p_{\text{Exact}}} \sqrt{\frac{\sum_i (p_{\text{SPH},i} - p_{\text{Exact}})^2}{\sum_i}} \qquad (27) \\ L_\infty(p_{\text{Exact}}) \equiv \frac{1}{p_{\text{Exact}}} \max_i (|p_{\text{SPH},i} - p_{\text{Exact}}|)$$

Table 1 shows that when the particle resolutions increase, both $L_2$ and $L_\infty$ norms decrease. We can conclude that the pressure gradient inside the bubble converges to the analytical solution $\Delta p = 0.837$ [Pa].
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 1.** Comparisons 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.</caption><thead><tr><td>Mesh Resolution</td><td>Absolute Value of L<sub>2</sub> Norm</td><td>Maximum L<sub>∞</sub> Norm</td></tr></thead><tbody><tr><td>60</td><td>0.049</td><td>0.937</td></tr><tr><td>100</td><td>0.044</td><td>0.854</td></tr><tr><td>140</td><td>0.025</td><td>0.821</td></tr></tbody></table>

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.76 [mm] × 5.76 [mm]) domain where droplet is characterised by $ρ_d = 250$ [kg/m³], $λ_d = 96 × 10^{-6}$ [m²/s] and the background fluid by $ρ_b = 500$ [kg/m³], $λ_b = 48 × 10^{-6}$ [m²/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,

$$ \sigma = \sigma_{ref} + \sigma_{T}(T - T_{ref}) \qquad (28) $$

where $σ_{ref}$ and $T_{ref}$ 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]. 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
<figure><img src="image_5.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 1.** Comparisons 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.</caption><thead><tr><td>Mesh Resolution</td><td>Absolute Value of L<sub>2</sub> Norm</td><td>Maximum L<sub>∞</sub> Norm</td></tr></thead><tbody><tr><td>60</td><td>0.049</td><td>0.937</td></tr><tr><td>100</td><td>0.044</td><td>0.854</td></tr><tr><td>140</td><td>0.025</td><td>0.821</td></tr></tbody></table>

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.76 [mm] × 5.76 [mm]) domain where droplet is characterised by $ρ_d = 250$ [kg/m³], $λ_d = 96 × 10^{-6}$ [m²/s] and the background fluid by $ρ_b = 500$ [kg/m³], $λ_b = 48 × 10^{-6}$ [m²/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,

$$ \sigma = \sigma_{ref} + \sigma_{T}(T - T_{ref}) \qquad (28) $$

where $σ_{ref}$ and $T_{ref}$ 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]. 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 ($\sigma_T \times \nabla T = 0.4$ [N/m²]). 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 is 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].

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Profile of the volumetric Marangoni force along the normal direction to the interface using a resolution of $Lo = 180$ [µm], $Lo = 90$ [µm] and $Lo = 45$ [µm].</figcaption></figure>

<table><caption>**Table 2.** Comparison between SPH results of the integrated Marangoni force and analytical solution.</caption><thead><tr><th>Particle Spacing [µm]</th><th>Numerical Result [N]</th><th>Analytical Solution [N]</th><th>Relative Error Percentage</th></tr></thead><tbody><tr><td>45</td><td>0.394421</td><td>0.4</td><td>0.013947</td></tr><tr><td>90</td><td>0.394411</td><td>0.4</td><td>0.013971</td></tr><tr><td>180</td><td>0.394425</td><td>0.4</td><td>0.013975</td></tr></tbody></table>

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 $\nabla T$ and $\partial\sigma$ denote the temperature and the surface gradient operator, respectively, the thermocapillary tangential stress writes as

$$ \nabla_s \sigma = \frac{\partial \sigma}{\partial T} \nabla_s T \qquad (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] × 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
interfacial Marangoni force direction is downwards. To compare the magnitude numerical results with exact solution ($\sigma_T \times \nabla T = 0.4$ [N/m²]). 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].

<figure><img src="image_6.png" /><figcaption>**Figure 3.** Profile of the volumetric Marangoni force along the normal direction to the interface using a resolution of $Lo = 180$ [µm], $Lo = 90$ [µm] and $Lo = 45$ [µm].</figcaption></figure>

<table><caption>**Table 2.** Comparison between SPH results of the integrated Marangoni force and analytical solution.</caption><thead><tr><td>Particle Spacing [µm]</td><td>Numerical Result [N]</td><td>Analytical Solution [N]</td><td>Relative Error Percentage</td></tr></thead><tbody><tr><td>45</td><td>0.394421</td><td>0.4</td><td>0.013947</td></tr><tr><td>90</td><td>0.394411</td><td>0.4</td><td>0.013971</td></tr><tr><td>180</td><td>0.394425</td><td>0.4</td><td>0.013975</td></tr></tbody></table>

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 $\nabla T$ and $\partial\sigma$ denote the temperature and the surface gradient operator, respectively, the thermocapillary tangential stress writes as

$$ \nabla_s \sigma = \frac{\partial \sigma}{\partial T} \nabla_s T \qquad (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] × 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 $\rho_d = 250$ [kg/m³] and $\mu_d = 0.012$ [Pa·s] for the droplet and $\rho_b = 500$ [kg/m³] and $\mu_b = 0.024$ [Pa·s]) for the background fluid, we consider that the heat conductivity of the droplet ($1.2 \times 10^{-6}$ [W/mK]) is half of the background fluid. With constant surface tension $\sigma_{ref} = 0.01$ [N/m] and the rate of the change of the surface tension with temperature $\sigma_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^* = \frac{U}{U_c}$ where the characteristic velocity $U_c = \frac{\sigma_T |\nabla T|R}{\mu_b}$ and the dimensionless time $t^* = \frac{t \cdot U_c}{R}$.

<figure><img src="image_7.png" /><figcaption>Figure 4. Comparison between droplet migration velocity in this work and Ma and Bothe [41].</figcaption></figure>

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 time-evolution 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 permittivity ratios of the phases. They are defined as $S = \epsilon_d/\epsilon_b$ and $R = \zeta_d/\zeta_b$ for the electrical permittivity 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
at the stationary state with $\rho_d = 250$ [kg/m³] and $\mu_d = 0.012$ [Pa·s] for the droplet and $\rho_b = 500$ [kg/m³] and $\mu_b = 0.024$ [Pa·s]) for the background fluid, we consider that the heat conductivity of the droplet ($1.2 \times 10^{-6}$ [W/mK]) is half of the background fluid. With constant surface tension $\sigma_{ref} = 0.01$ [N/m] and the rate of the change of the surface tension with temperature $\sigma_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^* = \frac{U}{U_c}$ where the characteristic velocity $U_c = \frac{\sigma_T |\nabla T|R}{\mu_b}$ and the dimensionless time $t^* = \frac{t \cdot U_c}{R}$.

<figure><img src="image_7.png" /><figcaption>Figure 4. Comparison between droplet migration velocity in this work and Ma and Bothe [41].</figcaption></figure>

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 time-evolution 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 permittivity ratios of the phases. They are defined as $S = \epsilon_d/\epsilon_b$ and $R = \zeta_d/\zeta_b$ for the electrical permittivity 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.

<figure><img src="image_8.png" /><figcaption>**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].</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between the orientation of induced flows inside and outside of the droplet with $S > R$: (a) Numerical results (b) Taylor's theory [42].</figcaption></figure>
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.

<figure><img src="image_8.png" /><figcaption>**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].</figcaption></figure>

<figure><img src="image_9.png" /><figcaption>**Figure 6.** Comparison between the orientation of induced flows inside and outside of the droplet with $S > R$: (a) Numerical results (b) Taylor's theory [42].</figcaption></figure>
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 3.** Simulation parameters for EHD cases using small deformation theory.</caption><thead><tr><th>Case</th><th>ε<sub>d</sub> [F/m]</th><th>ζ<sub>d</sub> [S/m]</th><th>S</th><th>R</th><th>σ [N/m]</th></tr></thead><tbody><tr><td>S<sub>D1</sub></td><td>0.3</td><td>40</td><td>0.5</td><td>2</td><td>0.01</td></tr><tr><td>S<sub>D2</sub></td><td>00.5</td><td>40</td><td>0.5</td><td>2</td><td>0.01</td></tr><tr><td>S<sub>D3</sub></td><td>0.5</td><td>150</td><td>0.5</td><td>3</td><td>0.01</td></tr><tr><td>S<sub>D4</sub></td><td>0.5</td><td>1</td><td>0.5</td><td>0.05</td><td>0.01</td></tr><tr><td>S<sub>D5</sub></td><td>3</td><td>10</td><td>5</td><td>0.5</td><td>0.03</td></tr></tbody></table>

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 = \frac{9 f_{dT} E_o^2 \epsilon_d R}{8(2 + \mathcal{R})^2 \sigma} \qquad (30)$$

where $f_{dT}$ is the discriminating function evaluated as

$$f_{dT} = \mathcal{R}^2 + 1 - 3.5\mathcal{S} + 1.5\mathcal{R} \qquad (31)$$

For the same problem, Feng [44] suggests the following relation

$$D_F = \frac{f_{dF} E_o^2 \epsilon_b R}{3(1 + \mathcal{R})^2 \mathcal{S} \sigma} \qquad (32)$$

where $f_{dF}$ is estimated from

$$f_{dF} = \mathcal{R}^2 + 1 - 3\mathcal{S} + \mathcal{R} \qquad (33)$$

In Equations (30) and (32), $\mathcal{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 = (\varphi^+ - \varphi^-)/h$, with $h$ being the height of the domain. Numerically, the
<figure><img src="image_10.png" /><figcaption>**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.</figcaption></figure>

<table><caption>**Table 3.** Simulation parameters for EHD cases using small deformation theory.</caption><thead><tr><th>Case</th><th>ε<sub>d</sub> [F/m]</th><th>ζ<sub>d</sub> [S/m]</th><th>S</th><th>R</th><th>σ [N/m]</th></tr></thead><tbody><tr><td>S<sub>D1</sub></td><td>0.3</td><td>40</td><td>0.5</td><td>2</td><td>0.01</td></tr><tr><td>S<sub>D2</sub></td><td>00.5</td><td>40</td><td>0.5</td><td>2</td><td>0.01</td></tr><tr><td>S<sub>D3</sub></td><td>0.5</td><td>150</td><td>0.5</td><td>3</td><td>0.01</td></tr><tr><td>S<sub>D4</sub></td><td>0.5</td><td>1</td><td>0.5</td><td>0.05</td><td>0.01</td></tr><tr><td>S<sub>D5</sub></td><td>3</td><td>10</td><td>5</td><td>0.5</td><td>0.03</td></tr></tbody></table>

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 = \frac{9 f_{dT} E_o^2 \epsilon_d R}{8(2 + \mathcal{R})^2 \sigma} \qquad (30)$$

where $f_{dT}$ is the discriminating function evaluated as

$$f_{dT} = \mathcal{R}^2 + 1 - 3.5\mathcal{S} + 1.5\mathcal{R} \qquad (31)$$

For the same problem, Feng [44] suggests the following relation

$$D_F = \frac{f_{dF} E_o^2 \epsilon_b R}{3(1 + \mathcal{R})^2 \mathcal{S} \sigma} \qquad (32)$$

where $f_{dF}$ is estimated from

$$f_{dF} = \mathcal{R}^2 + 1 - 3\mathcal{S} + \mathcal{R} \qquad (33)$$

In Equations (30) and (32), $\mathcal{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 = (\varphi^+ - \varphi^-)/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 = \frac{A - B}{A + B} \qquad (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.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³] 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 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.

<table><caption>**Table 4.** Comparison 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.</caption><thead><tr><th>Case</th><th>D<sub>N</sub></th><th>D<sub>T</sub></th><th>D<sub>F</sub></th><th>Bubble Shape</th></tr></thead><tbody><tr><td>S<sub>D1</sub></td><td>0.077</td><td>0.065</td><td>0.061</td><td></td></tr><tr><td>S<sub>D2</sub></td><td>0.121</td><td>0.109</td><td>0.101</td><td></td></tr><tr><td>S<sub>D3</sub></td><td>0.157</td><td>0.143</td><td>0.119</td><td></td></tr><tr><td>S<sub>D4</sub></td><td>-0.040</td><td>-0.045</td><td>-0.063</td><td></td></tr><tr><td>S<sub>D5</sub></td><td>-0.097</td><td>-0.139</td><td>-0.196</td><td></td></tr></tbody></table>

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

$$v_{r,in} = U \left[ \left( \frac{r}{R} \right)^3 - \left( \frac{r}{R} \right) \right] \cos 2\theta \qquad (35)$$

$$v_{\theta,in} = U \left[ \left( \frac{r}{R} \right) - 2 \left( \frac{r}{R} \right)^3 \right] \sin 2\theta \qquad (36)$$
droplet deformation parameter *D* can be defined based on the droplet's deviation from the circular shape,

$$D_N = \frac{A - B}{A + B} \qquad (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*<sub>*N*</sub> = 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*<sub>*N*</sub> 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.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 ρ<sub>*d*</sub> = 1000 [kg/m<sup>3</sup>] and μ<sub>*d*</sub> = 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 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.

<table><caption>**Table 4.** Comparison 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.</caption><thead><tr><th>Case</th><th>D<sub>N</sub></th><th>D<sub>T</sub></th><th>D<sub>F</sub></th><th>Bubble Shape</th></tr></thead><tbody><tr><td>S<sub>D1</sub></td><td>0.077</td><td>0.065</td><td>0.061</td><td></td></tr><tr><td>S<sub>D2</sub></td><td>0.121</td><td>0.109</td><td>0.101</td><td></td></tr><tr><td>S<sub>D3</sub></td><td>0.157</td><td>0.143</td><td>0.119</td><td></td></tr><tr><td>S<sub>D4</sub></td><td>-0.040</td><td>-0.045</td><td>-0.063</td><td></td></tr><tr><td>S<sub>D5</sub></td><td>-0.097</td><td>-0.139</td><td>-0.196</td><td></td></tr></tbody></table>

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

$$v_{r,in} = U \left[ \left( \frac{r}{R} \right)^3 - \left( \frac{r}{R} \right) \right] \cos 2\theta \qquad (35)$$

$$v_{\theta,in} = U \left[ \left( \frac{r}{R} \right) - 2 \left( \frac{r}{R} \right)^3 \right] \sin 2\theta \qquad (36)$$
$$v_{r,out} = U \left[ \left( \frac{R}{r} \right) - \left( \frac{R}{r} \right)^3 \right] \cos 2\theta \quad (37)$$

$$v_{\theta,out} = -U \left[ \left( \frac{R}{r} \right)^3 \right] \sin 2\theta \quad (38)$$

where $r$ is the radial position, $v_\theta$ and $v_r$ are the tangential velocity and the radial velocity, respectively. The characteristic velocity $U$ can be evaluated as

$$U = \frac{\mathcal{R} - S}{2S(1 + \mathcal{R})^2} \frac{\epsilon_d E_0^2 R}{\mu_d + \mu_b} \quad (39)$$

Our numerical results are in agreement with the analytical solutions for both the velocity $\theta = 0^\circ$ and $\theta = 45^\circ$ 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_0 = 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 \neq R$) and causes small deviation between the numerical and analytical results. Based on the equations of tangential components (36) and (38), having ($\sin 2\theta = 0$) leading to $v_\theta = 0$ [m/s] and maximum values of $v_r$ that is in agreement with the observations in Figure 8.

<figure><img src="image_11.png" /><figcaption>Figure 8. Components of the velocity profile at (a) $\theta = 0^\circ$ and (b) $\theta = 45^\circ$ compared to the formulae given by Feng [44].</figcaption></figure>

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
$$v_{r,out} = U \left[ \left( \frac{R}{r} \right) - \left( \frac{R}{r} \right)^3 \right] \cos 2\theta \quad (37)$$

$$v_{\theta,out} = -U \left[ \left( \frac{R}{r} \right)^3 \right] \sin 2\theta \quad (38)$$

where $r$ is the radial position, $v_\theta$ and $v_r$ are the tangential velocity and the radial velocity, respectively. The characteristic velocity $U$ can be evaluated as

$$U = \frac{\mathcal{R} - \mathcal{S}}{2\mathcal{S}(1 + \mathcal{R})^2} \frac{\epsilon_d E_0^2 R}{\mu_d + \mu_b} \quad (39)$$

Our numerical results are in agreement with the analytical solutions for both the velocity $\theta = 0^\circ$ and $\theta = 45^\circ$ 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_0 = 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 \neq R$) and causes small deviation between the numerical and analytical results. Based on the equations of tangential components (36) and (38), having ($\sin 2\theta = 0$) leading to $v_\theta = 0$ [m/s] and maximum values of $v_r$ that is in agreement with the observations in Figure 8.

<figure><img src="image_11.png" /><figcaption>Figure 8. Components of the velocity profile at (a) $\theta = 0^\circ$ and (b) $\theta = 45^\circ$ compared to the formulae given by Feng [44].</figcaption></figure>

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 loses 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 $\frac{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 time $t = \frac{\epsilon_d + 2\epsilon_b}{\zeta_d + 2\zeta_b}$. It is worth reminding that electric charge accumulation at the fluid-fluid interface happens when each phase has a different charge relaxation time $\tau = \frac{\epsilon}{\zeta}$ where $\epsilon$ and $\zeta$ 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.

The coupled EHD-TC problem setup consists of a square domain of size $L = 0.04$ [m] with $\frac{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 $\phi = 0.04$ [V] is imposed on the top boundary while other boundaries are set to $\phi = 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 ($\nabla p \cdot \mathbf{n} = 0$) where $\mathbf{n}$ is the normal direction to the given boundary. The density and viscosity of are chosen to be $\rho_d = 250$ [kg/m³], $\mu_d = 0.12$ [Pa·s] for the bubble and $\rho_b = 500$ [kg/m³], $\mu_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 5. The time evolution of the droplet subject to thermocapillary flow and EHD forces is illustrated in Figure 9.

<table><caption>**Table 5.** Simulation parameters of the coupled EHD-TC cases where vertical (V) and horizontal (H) electric fields are imposed.</caption><thead><tr><th>Case</th><th>Electric Field Direction</th><th>ε<sub>b</sub> [F/m]</th><th>ζ<sub>b</sub> [S/m]</th><th>ε<sub>d</sub> [F/m]</th><th>ζ<sub>d</sub> [S/m]</th></tr></thead><tbody><tr><td>1</td><td>V</td><td>1</td><td>50</td><td>0.5</td><td>150</td></tr><tr><td>2</td><td>H</td><td>1</td><td>50</td><td>0.5</td><td>150</td></tr><tr><td>3</td><td>V</td><td>0.5</td><td>150</td><td>1</td><td>50</td></tr></tbody></table>
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 loses 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 $\frac{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 time $t = \frac{\epsilon_d + 2\epsilon_b}{\zeta_d + 2\zeta_b}$. It is worth reminding that electric charge accumulation at the fluid-fluid interface happens when each phase has a different charge relaxation time $\tau = \frac{\epsilon}{\zeta}$ where $\epsilon$ and $\zeta$ 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.

The coupled EHD-TC problem setup consists of a square domain of size $L = 0.04$ [m] with $\frac{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 $\phi = 0.04$ [V] is imposed on the top boundary while other boundaries are set to $\phi = 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 ($\nabla p \cdot \mathbf{n} = 0$) where $\mathbf{n}$ is the normal direction to the given boundary. The density and viscosity of are chosen to be $\rho_d = 250$ [kg/m³], $\mu_d = 0.12$ [Pa·s] for the bubble and $\rho_b = 500$ [kg/m³], $\mu_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 5. The time evolution of the droplet subject to thermocapillary flow and EHD forces is illustrated in Figure 9.

<table><caption>**Table 5.** Simulation parameters of the coupled EHD-TC cases where vertical (V) and horizontal (H) electric fields are imposed.</caption><thead><tr><th>Case</th><th>Electric Field Direction</th><th>ε<sub>b</sub> [F/m]</th><th>ζ<sub>b</sub> [S/m]</th><th>ε<sub>d</sub> [F/m]</th><th>ζ<sub>d</sub> [S/m]</th></tr></thead><tbody><tr><td>1</td><td>V</td><td>1</td><td>50</td><td>0.5</td><td>150</td></tr><tr><td>2</td><td>H</td><td>1</td><td>50</td><td>0.5</td><td>150</td></tr><tr><td>3</td><td>V</td><td>0.5</td><td>150</td><td>1</td><td>50</td></tr></tbody></table>
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 generates 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 ($\mu_d$ = 0.12 [Pa·s]) is chosen close to continues phase viscosity ($\mu_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 $\sigma_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
<figure><img src="image_12.png" /><figcaption>**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.</figcaption></figure>

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 generates 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 ($\mu_d$ = 0.12 [Pa·s]) is chosen close to continues phase viscosity ($\mu_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 $\sigma_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 $\mathcal{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 $\mathcal{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.

**5. 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 permittivity, 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 $\mathcal{R}$ and $S$ corresponding to electrical permittivity and electrical conductivity, respectively as well as the orientation of the applied electric field with respect to thermal gradient.

* For $\mathcal{R} > S$ and thermal and electric gradient applied both in the vertical direction, a prolate deformation of the droplet is observed.

* Having $\mathcal{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 $\mathcal{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
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 $\mathcal{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 $\mathcal{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.

**5. 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 permittivity, 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 $\mathcal{R}$ and $S$ corresponding to electrical permittivity and electrical conductivity, respectively as well as the orientation of the applied electric field with respect to thermal gradient.

* For $\mathcal{R} > S$ and thermal and electric gradient applied both in the vertical direction, a prolate deformation of the droplet is observed.

* Having $\mathcal{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 $\mathcal{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.

**Author Contributions:** Conceptualization, M.S.S.; Methodology, M.S.S. and M.H.-H.; Software, F.A. and M.H.-H.; Investigation: F.A.; Validation, F.A.; Formal Analysis, F.A. and M.H.-H.; Resources, A.H., M.S.S. and U.N.; Writing—Original Draft Preparation, F.A. and M.H.-H.; Writing—Review & Editing, F.A.; Visualization, F.A.; Supervision, M.S.S. and A.H.; Project Administration, M.S.S., A.H. and U.N.; Funding Acquisition, M.S.S. and A.H. All authors have read and agreed to the published version of the manuscript.

**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**

<table><tr><td>χ</td><td>Curvature</td></tr><tr><td>κ</td><td>Thermal diffusivity (m²/s)</td></tr><tr><td>μ</td><td>Dynamic viscosity (Pa·s)</td></tr><tr><td>σ</td><td>Surface tension (N/m)</td></tr><tr><td>R</td><td>Initial droplet radius (m)</td></tr><tr><td>F<sub>e</sub></td><td>Electric force (N)</td></tr><tr><td>g</td><td>Gravity acceleration (m/s²)</td></tr><tr><td>c<sub>p</sub></td><td>Specific heat capacity at constant pressure (J/kg K)</td></tr><tr><td>ρ</td><td>Density (kg/m³)</td></tr><tr><td>ε</td><td>Electric permittivity (F/m)</td></tr><tr><td>ζ</td><td>Electric conductivity (S/m)</td></tr><tr><td>p</td><td>Static pressure (kg/m s²)</td></tr><tr><td>T</td><td>Temperature (K)</td></tr><tr><td>T&#773;</td><td>Maxwell stress tensor</td></tr><tr><td>M<sub>a</sub></td><td>Marangoni number</td></tr><tr><td>v</td><td>Local average fluid velocity</td></tr><tr><td>Π</td><td>Stress tensor</td></tr><tr><td>q<sup>v</sup></td><td>Free electric charge density (C/m³)</td></tr><tr><td>Ω</td><td>Domain of integration of the kernel function</td></tr><tr><td>δ</td><td>Dirac-delta function</td></tr><tr><td>C<sub>w</sub></td><td>Normalisation constant of the kernel function</td></tr><tr><td>W(x, h)</td><td>Kernel function</td></tr><tr><td>h</td><td>Smoothing length of the kernel function</td></tr><tr><td>V<sub>i</sub></td><td>Volume of the particle i</td></tr><tr><td>Φ&#773;<sub>ij</sub></td><td>Inter-particle average of the electric potential</td></tr><tr><td>λ<sub>i</sub></td><td>Thermal conductivity (W/m K)</td></tr><tr><td>L<sub>o</sub></td><td>Characteristic length (m)</td></tr><tr><td>τ</td><td>Bulk relaxation time (s)</td></tr></table>
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.

**Author Contributions:** Conceptualization, M.S.S.; Methodology, M.S.S. and M.H.-H.; Software, F.A. and M.H.-H.; Investigation: F.A.; Validation, F.A.; Formal Analysis, F.A. and M.H.-H.; Resources, A.H., M.S.S. and U.N.; Writing—Original Draft Preparation, F.A. and M.H.-H.; Writing—Review & Editing, F.A.; Visualization, F.A.; Supervision, M.S.S. and A.H.; Project Administration, M.S.S., A.H. and U.N.; Funding Acquisition, M.S.S. and A.H. All authors have read and agreed to the published version of the manuscript.

**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

<table><tr><td>χ</td><td>Curvature</td></tr><tr><td>κ</td><td>Thermal diffusivity (m²/s)</td></tr><tr><td>μ</td><td>Dynamic viscosity (Pa·s)</td></tr><tr><td>σ</td><td>Surface tension (N/m)</td></tr><tr><td>R</td><td>Initial droplet radius (m)</td></tr><tr><td>F<sub>e</sub></td><td>Electric force (N)</td></tr><tr><td>g</td><td>Gravity acceleration (m/s²)</td></tr><tr><td>c<sub>p</sub></td><td>Specific heat capacity at constant pressure (J/kg K)</td></tr><tr><td>ρ</td><td>Density (kg/m³)</td></tr><tr><td>ε</td><td>Electric permittivity (F/m)</td></tr><tr><td>ζ</td><td>Electric conductivity (S/m)</td></tr><tr><td>p</td><td>Static pressure (kg/m s²)</td></tr><tr><td>T</td><td>Temperature (K)</td></tr><tr><td>T&#773;</td><td>Maxwell stress tensor</td></tr><tr><td>M<sub>a</sub></td><td>Marangoni number</td></tr><tr><td>v</td><td>Local average fluid velocity</td></tr><tr><td>Π</td><td>Stress tensor</td></tr><tr><td>q<sup>v</sup></td><td>Free electric charge density (C/m³)</td></tr><tr><td>Ω</td><td>Domain of integration of the kernel function</td></tr><tr><td>δ</td><td>Dirac-delta function</td></tr><tr><td>C<sub>w</sub></td><td>Normalisation constant of the kernel function</td></tr><tr><td>W(x, h)</td><td>Kernel function</td></tr><tr><td>h</td><td>Smoothing length of the kernel function</td></tr><tr><td>V<sub>i</sub></td><td>Volume of the particle i</td></tr><tr><td>Φ&#773;<sub>ij</sub></td><td>Inter-particle average of the electric potential</td></tr><tr><td>λ<sub>i</sub></td><td>Thermal conductivity (W/m K)</td></tr><tr><td>L<sub>o</sub></td><td>Characteristic length (m)</td></tr><tr><td>τ</td><td>Bulk relaxation time (s)</td></tr></table>
<table><tr><td>R</td><td>Electrical conductivity of the droplet to bulk ratio</td></tr><tr><td>E<sub>o</sub></td><td>Electric field vector magnitude (V/m) set to be E<sub>o</sub> = 1 unless stated otherwise</td></tr><tr><td>S</td><td>Electrical permittivity of the droplet to bulk ratio</td></tr><tr><td>φ</td><td>Electric potential across the domain</td></tr><tr><td>U</td><td>Characteristic velocity</td></tr><tr><td>+</td><td>High-electric potential</td></tr><tr><td>−</td><td>Low-electric potential</td></tr><tr><td>b</td><td>Bulk conditions</td></tr><tr><td>d</td><td>Droplet conditions</td></tr><tr><td>v</td><td>viscous stress</td></tr><tr><td>c</td><td>capillary stress</td></tr><tr><td>ref</td><td>Reference conditions</td></tr><tr><td>T</td><td>Gradient with respect to temperature</td></tr><tr><td>*</td><td>Dimensionless quantity</td></tr><tr><td>H</td><td>Hot conditions</td></tr><tr><td>C</td><td>Cold conditions</td></tr><tr><td>θ</td><td>Tangential component</td></tr></table>

## References

1. Eow, J.; Ghadiri, M.; Sharif, A.; Williams, T. Electrostatic enhancement of coalescence of water droplets in oil: A review of the current understanding. *Chem. Eng. J.* **2001**, *84*, 173–192. [CrossRef]

2. Zolfaghari, R.; FakhruRazi, A.; Luqman Chuah, A.; Elnashaie, S.S.; Pendashteh, A. Demulsification techniques of water-in-oil and oil-in-water emulsions in petroleum industry. *Sep. Purif. Technol.* **2016**, *170*, 377–407. [CrossRef]

3. Yue, X.; Li, Z.; Zhang, T.; Yang, D.; Qiu, F. Design and fabrication of superwetting fiber-based membranes for oil/water separation applications. *Chem. Eng. J.* **2019**, *364*, 292–309. [CrossRef]

4. Das, R.; Cleary, P. Three-dimensional modelling of coupled flow dynamics, heat transfer and residual stress generation in arc welding processes using the mesh-free SPH method. *J. Comput. Sci.* **2016**, *16*, 200–216. [CrossRef]

5. Brennen, C.E. *Fundamentals of Multiphase Flow*; Cambridge University Press: Cambridge, UK, 2005; ISBN 9780511807169. [CrossRef]

6. Necati, O.M.; Orlande, H.R.B.; Colaco, M.J.; Cotta, R.M. *Finite Difference Methods in Heat Transfer*; Taylor and Francis, Chemical Rubber Company (CRC) Press: Boca Raton, FL, USA, 2017.

7. Lewis, R.W.; Schrefler, B.A. *The Finite Element Method in the Deformation and Consolidation of Porous Media*; Wiley: Hoboken, NJ, USA, 2000.

8. Moukalled, F.; Mangani, L.; Darwish, M. *The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® and Matlab®*; Springer: Berlin/Heidelberg, Germany, 2016; Volume 113, 791p.

9. Sethian, J.A.; Smereka, P. Level Set Methods for Fluid Interfaces. *Annu. Rev. Fluid Mech.* **2003**, *35*, 341–372. doi: 10.1146/annurev.fluid.35.101101.161105. [CrossRef]

10. Subramaniam, S. *Lagrangian-Eulerian methods for multiphase flows*. *Prog. Energy Combust. Sci.* **2013**, *39*, 215–245. [CrossRef]

11. Besagni, G.; Inzoli, F.; Ziegenhein, T. Two-Phase Bubble Columns: A Comprehensive Review. *ChemEngineering* **2018**, *2*, 13. [CrossRef]

12. Ishii, M. *Thermo-Fluid Dynamic Theory of Two-Phase Flow*; Eyrolles: Paris, France, 1975.

13. Li, S.; Liu, W.K. Meshfree and particle methods and their applications. *Appl. Mech. Rev.* **2002**, *55*, 1–34. [CrossRef]

14. P Fries, T.; G Matthies, H. *Classification and Overview of Meshfree Methods*; Technical University Braunschweig: Brunswick, Germany, 2003.

15. Shadloo, M.; Oger, G.; Touzé, D.L. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. *Comput. Fluids* **2016**, *136*, 11–34. [CrossRef]

16. Nazaripoor, H. Electrohydrodynamic and Thermocapillary Instability of Thin Liquid Films. Ph.D. Thesis, University of Alberta, Edmonton, AB, Canada, 2016.

17. Corbett, A.; Kumar, S. Combined thermal and electrohydrodynamic patterning of thin liquid films. *J. Eng. Math.* **2015**, *94*, 81–96. [CrossRef]

18. Lau, C.Y.; Russel, W.B. Fundamental Limitations on Ordered Electrohydrodynamic Patterning. *Macromolecules* **2011**, *44*, 7746–7751. [CrossRef]

19. Hopp-Hirschler, M.; Shadloo, M.S.; Nieken, U. A Smoothed Particle Hydrodynamics approach for thermo-capillary flows. *Comput. Fluids* **2018**, *176*, 1–19. [CrossRef]

20. Brackbill, J.; Kothe, D.; Zemach, C. A continuum method for modeling surface tension. *J. Comput. Phys.* **1992**, *100*, 335–354. [CrossRef]

21. Saville, D.A. Electrohydrodynamics: The Taylor-Melcher Leaky Dielectric Model. *Annu. Rev. Fluid Mech.* **1997**, *29*, 27–64. [CrossRef]
<table><tr><td>R</td><td>Electrical conductivity of the droplet to bulk ratio</td></tr><tr><td>E<sub>o</sub></td><td>Electric field vector magnitude (V/m) set to be E<sub>o</sub> = 1 unless stated otherwise</td></tr><tr><td>S</td><td>Electrical permittivity of the droplet to bulk ratio</td></tr><tr><td>φ</td><td>Electric potential across the domain</td></tr><tr><td>U</td><td>Characteristic velocity</td></tr><tr><td>+</td><td>High-electric potential</td></tr><tr><td>-</td><td>Low-electric potential</td></tr><tr><td>b</td><td>Bulk conditions</td></tr><tr><td>d</td><td>Droplet conditions</td></tr><tr><td>v</td><td>viscous stress</td></tr><tr><td>c</td><td>capillary stress</td></tr><tr><td>ref</td><td>Reference conditions</td></tr><tr><td>T</td><td>Gradient with respect to temperature</td></tr><tr><td>*</td><td>Dimensionless quantity</td></tr><tr><td>H</td><td>Hot conditions</td></tr><tr><td>C</td><td>Cold conditions</td></tr><tr><td>θ</td><td>Tangential component</td></tr></table>

## References

1. Eow, J.; Ghadiri, M.; Sharif, A.; Williams, T. Electrostatic enhancement of coalescence of water droplets in oil: A review of the current understanding. *Chem. Eng. J.* **2001**, *84*, 173–192. [CrossRef]

2. Zolfaghari, R.; FakhruRazi, A.; Luqman Chuah, A.; Elnashaie, S.S.; Pendashteh, A. Demulsification techniques of water-in-oil and oil-in-water emulsions in petroleum industry. *Sep. Purif. Technol.* **2016**, *170*, 377–407. [CrossRef]

3. Yue, X.; Li, Z.; Zhang, T.; Yang, D.; Qiu, F. Design and fabrication of superwetting fiber-based membranes for oil/water separation applications. *Chem. Eng. J.* **2019**, *364*, 292–309. [CrossRef]

4. Das, R.; Cleary, P. Three-dimensional modelling of coupled flow dynamics, heat transfer and residual stress generation in arc welding processes using the mesh-free SPH method. *J. Comput. Sci.* **2016**, *16*, 200–216. [CrossRef]

5. Brennen, C.E. *Fundamentals of Multiphase Flow*; Cambridge University Press: Cambridge, UK, 2005; ISBN 9780511807169. [CrossRef]

6. Necati, O.M.; Orlande, H.R.B.; Colaco, M.J.; Cotta, R.M. *Finite Difference Methods in Heat Transfer*; Taylor and Francis, Chemical Rubber Company (CRC) Press: Boca Raton, FL, USA, 2017.

7. Lewis, R.W.; Schrefler, B.A. *The Finite Element Method in the Deformation and Consolidation of Porous Media*; Wiley: Hoboken, NJ, USA, 2000.

8. Moukalled, F.; Mangani, L.; Darwish, M. *The Finite Volume Method in Computational Fluid Dynamics: An Advanced Introduction with OpenFOAM® and Matlab®*; Springer: Berlin/Heidelberg, Germany, 2016; Volume 113, 791p.

9. Sethian, J.A.; Smereka, P. Level Set Methods for Fluid Interfaces. *Annu. Rev. Fluid Mech.* **2003**, *35*, 341–372. doi: 10.1146/annurev.fluid.35.101101.161105. [CrossRef]

10. Subramaniam, S. *Lagrangian-Eulerian methods for multiphase flows*. *Prog. Energy Combust. Sci.* **2013**, *39*, 215–245. [CrossRef]

11. Besagni, G.; Inzoli, F.; Ziegenhein, T. Two-Phase Bubble Columns: A Comprehensive Review. *ChemEngineering* **2018**, *2*, 13. [CrossRef]

12. Ishii, M. *Thermo-Fluid Dynamic Theory of Two-Phase Flow*; Eyrolles: Paris, France, 1975.

13. Li, S.; Liu, W.K. Meshfree and particle methods and their applications. *Appl. Mech. Rev.* **2002**, *55*, 1–34. [CrossRef]

14. P Fries, T.; G Matthies, H. *Classification and Overview of Meshfree Methods*; Technical University Braunschweig: Brunswick, Germany, 2003.

15. Shadloo, M.; Oger, G.; Touzé, D.L. Smoothed particle hydrodynamics method for fluid flows, towards industrial applications: Motivations, current state, and challenges. *Comput. Fluids* **2016**, *136*, 11–34. [CrossRef]

16. Nazaripoor, H. Electrohydrodynamic and Thermocapillary Instability of Thin Liquid Films. Ph.D. Thesis, University of Alberta, Edmonton, AB, Canada, 2016.

17. Corbett, A.; Kumar, S. Combined thermal and electrohydrodynamic patterning of thin liquid films. *J. Eng. Math.* **2015**, *94*, 81–96. [CrossRef]

18. Lau, C.Y.; Russel, W.B. Fundamental Limitations on Ordered Electrohydrodynamic Patterning. *Macromolecules* **2011**, *44*, 7746–7751. [CrossRef]

19. Hopp-Hirschler, M.; Shadloo, M.S.; Nieken, U. A Smoothed Particle Hydrodynamics approach for thermo-capillary flows. *Comput. Fluids* **2018**, *176*, 1–19. [CrossRef]

20. Brackbill, J.; Kothe, D.; Zemach, C. A continuum method for modeling surface tension. *J. Comput. Phys.* **1992**, *100*, 335–354. [CrossRef]

21. Saville, D.A. Electrohydrodynamics: The Taylor-Melcher Leaky Dielectric Model. *Annu. Rev. Fluid Mech.* **1997**, *29*, 27–64. [CrossRef]
22. Taylor, G. Electrically Driven Jets. *Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.* **1969**, *313*, 453–475.

23. Shadloo, M.S. Improved Multiphase Smoothed Particle Hydrodynamics. Ph.D. Thesis, Sabanci University, Tuzla, Bosnia and Herzegovina, 2013.

24. Cummins, S.J.; Rudman, M. An SPH Projection Method. *J. Comput. Phys.* **1999**, *152*, 584–607. [CrossRef]

25. Szewc, K.; Pozorski, J.; Minier, J.P. Analysis of the incompressibility constraint in the smoothed particle hydrodynamics method. *Int. J. Numer. Methods Eng.* **2012**, *92*, 343–369. [CrossRef]

26. Hirschler, M.; Kunz, P.; Huber, M.; Hahn, F.; Nieken, U. Open boundary conditions for ISPH and their application to micro-flow. *J. Comput. Phys.* **2016**, *307*, 614–633. [CrossRef]

27. Gingold, R.A.; Monaghan, J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. *Mon. Not. R. Astron. Soc.* **1977**, *181*, 375–389. [CrossRef]

28. Lucy, L.B. A numerical approach to the testing of the fission hypothesis. *Astron. J.* **1977**, *82*, 1013–1024. [CrossRef]

29. Ye, T.; Pan, D.; Huang, C.; Liu, M. Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications. *Phys. Fluids* **2019**, *31*, 011301. [CrossRef]

30. Monaghan, J.J. Smoothed Particle Hydrodynamics. *Annu. Rev. Astron. Astrophys.* **1992**, *30*, 543–574. [CrossRef]

31. Wendland, H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. *Adv. Comput. Math.* **1995**, *4*, 389–396. [CrossRef]

32. Shadloo, M.S.; Zainali, A.; Yildiz, M.; Suleman, A. A robust weakly compressible SPH method and its comparison with an incompressible SPH. *Int. J. Numer. Methods Eng.* **2012**, *89*, 939–956. [CrossRef]

33. Almasi, F.; Shadloo, M.; Hadjadj, A.; Ozbulut, M.; Tofighi, N.; Yildiz, M. Numerical simulations of multi-phase electro-hydrodynamics flows using a simple incompressible smoothed particle hydrodynamics method. *Comput. Math. Appl.* **2019**, *81*, 772–785. [CrossRef]

34. Monaghan, J.J. Smoothed particle hydrodynamics. *Rep. Prog. Phys.* **2005**, *68*, 1703–1759. [CrossRef]

35. Hu, X.; Adams, N. A multi-phase SPH method for macroscopic and mesoscopic flows. *J. Comput. Phys.* **2006**, *213*, 844–861. [CrossRef]

36. Kim, J.; Moin, P. Application of a fractional-step method to incompressible Navier-Stokes equations. *J. Comput. Phys.* **1985**, *59*, 308–323. [CrossRef]

37. Morris, J.P.; Fox, P.J.; Zhu, Y. Modeling low Reynolds number incompressible flows using SPH. *J. Comput. Phys.* **1997**, *136*, 214–226. [CrossRef]

38. Hu, X.; Adams, N. An incompressible multi-phase SPH method. *J. Comput. Phys.* **2007**, *227*, 264–278. [CrossRef]

39. Hosseini, S.M.; Feng, J.J. Pressure boundary conditions for computing incompressible flows with SPH. *J. Comput. Phys.* **2011**, *230*, 7473–7487. [CrossRef]

40. Tong, M.; Browne, D.J. An incompressible multi-phase smoothed particle hydrodynamics (SPH) method for modelling thermocapillary flow. *Int. J. Heat Mass Transf.* **2014**, *73*, 284–292. [CrossRef]

41. Ma, C.; Bothe, D. Direct numerical simulation of thermocapillary flow based on the Volume of Fluid method. *Int. J. Multiph. Flow* **2011**, *37*, 1045–1058. [CrossRef]

42. Taylor, G.I.; McEwan, A.D.; de Jong, L.N.J. Studies in electrohydrodynamics. I. The circulation produced in a drop by an electric field. *Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.* **1966**, *291*, 159–166. [CrossRef]

43. Shadloo, M.S.; Rahmat, A.; Yildiz, M. A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid. *Comput. Mech.* **2013**, *52*, 693–707. [CrossRef]

44. Feng, J.Q. Electrohydrodynamic behaviour of a drop subjected to a steady uniform electric field at finite electric Reynolds number. *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.* **1999**, *455*, 2245–2269. [CrossRef]
22. Taylor, G. Electrically Driven Jets. *Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.* **1969**, *313*, 453–475.

23. Shadloo, M.S. Improved Multiphase Smoothed Particle Hydrodynamics. Ph.D. Thesis, Sabanci University, Tuzla, Bosnia and Herzegovina, 2013.

24. Cummins, S.J.; Rudman, M. An SPH Projection Method. *J. Comput. Phys.* **1999**, *152*, 584–607. [CrossRef]

25. Szewc, K.; Pozorski, J.; Minier, J.P. Analysis of the incompressibility constraint in the smoothed particle hydrodynamics method. *Int. J. Numer. Methods Eng.* **2012**, *92*, 343–369. [CrossRef]

26. Hirschler, M.; Kunz, P.; Huber, M.; Hahn, F.; Nieken, U. Open boundary conditions for ISPH and their application to micro-flow. *J. Comput. Phys.* **2016**, *307*, 614–633. [CrossRef]

27. Gingold, R.A.; Monaghan, J.J. Smoothed particle hydrodynamics: theory and application to non-spherical stars. *Mon. Not. R. Astron. Soc.* **1977**, *181*, 375–389. [CrossRef]

28. Lucy, L.B. A numerical approach to the testing of the fission hypothesis. *Astron. J.* **1977**, *82*, 1013–1024. [CrossRef]

29. Ye, T.; Pan, D.; Huang, C.; Liu, M. Smoothed particle hydrodynamics (SPH) for complex fluid flows: Recent developments in methodology and applications. *Phys. Fluids* **2019**, *31*, 011301. [CrossRef]

30. Monaghan, J.J. Smoothed Particle Hydrodynamics. *Annu. Rev. Astron. Astrophys.* **1992**, *30*, 543–574. [CrossRef]

31. Wendland, H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. *Adv. Comput. Math.* **1995**, *4*, 389–396. [CrossRef]

32. Shadloo, M.S.; Zainali, A.; Yildiz, M.; Suleman, A. A robust weakly compressible SPH method and its comparison with an incompressible SPH. *Int. J. Numer. Methods Eng.* **2012**, *89*, 939–956. [CrossRef]

33. Almasi, F.; Shadloo, M.; Hadjadj, A.; Ozbulut, M.; Tofighi, N.; Yildiz, M. Numerical simulations of multi-phase electro-hydrodynamics flows using a simple incompressible smoothed particle hydrodynamics method. *Comput. Math. Appl.* **2019**, *81*, 772–785. [CrossRef]

34. Monaghan, J.J. Smoothed particle hydrodynamics. *Rep. Prog. Phys.* **2005**, *68*, 1703–1759. [CrossRef]

35. Hu, X.; Adams, N. A multi-phase SPH method for macroscopic and mesoscopic flows. *J. Comput. Phys.* **2006**, *213*, 844–861. [CrossRef]

36. Kim, J.; Moin, P. Application of a fractional-step method to incompressible Navier-Stokes equations. *J. Comput. Phys.* **1985**, *59*, 308–323. [CrossRef]

37. Morris, J.P.; Fox, P.J.; Zhu, Y. Modeling low Reynolds number incompressible flows using SPH. *J. Comput. Phys.* **1997**, *136*, 214–226. [CrossRef]

38. Hu, X.; Adams, N. An incompressible multi-phase SPH method. *J. Comput. Phys.* **2007**, *227*, 264–278. [CrossRef]

39. Hosseini, S.M.; Feng, J.J. Pressure boundary conditions for computing incompressible flows with SPH. *J. Comput. Phys.* **2011**, *230*, 7473–7487. [CrossRef]

40. Tong, M.; Browne, D.J. An incompressible multi-phase smoothed particle hydrodynamics (SPH) method for modelling thermocapillary flow. *Int. J. Heat Mass Transf.* **2014**, *73*, 284–292. [CrossRef]

41. Ma, C.; Bothe, D. Direct numerical simulation of thermocapillary flow based on the Volume of Fluid method. *Int. J. Multiph. Flow* **2011**, *37*, 1045–1058. [CrossRef]

42. Taylor, G.I.; McEwan, A.D.; de Jong, L.N.J. Studies in electrohydrodynamics. I. The circulation produced in a drop by an electric field. *Proc. R. Soc. Lond. Ser. A Math. Phys. Sci.* **1966**, *291*, 159–166. [CrossRef]

43. Shadloo, M.S.; Rahmat, A.; Yildiz, M. A smoothed particle hydrodynamics study on the electrohydrodynamic deformation of a droplet suspended in a neutrally buoyant Newtonian fluid. *Comput. Mech.* **2013**, *52*, 693–707. [CrossRef]

44. Feng, J.Q. Electrohydrodynamic behaviour of a drop subjected to a steady uniform electric field at finite electric Reynolds number. *Proc. R. Soc. Lond. Ser. A Math. Phys. Eng. Sci.* **1999**, *455*, 2245–2269. [CrossRef]