Non‑local modelling of multiphase flow wetting and thermo‑capillary flow using peridynamic differential operator

Interfaces in multiphase flows are affected by surface tension, and when temperature gradients occur in the flow domain, tangential surface tensions along the interface also arise. As the behaviour of fluids contacting on a solid surface is also governed by surface tension, the description of the wetting phenomenon is challenging. Peridynamic differential operator (PDDO) can express partial differentials of any order by integral equations. Therefore, the governing equations for multiphase fluid motion, such as the Navier–Stokes equations and energy equations, can be reformulated in terms of integral equations. In this study, a novel non-local method is developed for modelling the multiphase fluid flow motion using the PDDO, and the thermal effect on surface tension force is considered. To describe the surface tension forces in the normal and tangential directions, the non-local form of the continuum surface force (CFS) model is presented. Besides, to overcome the inaccuracy of the unit normal vectors at the three-phase flow intersection region, an additional treatment for this region is presented. Finally, several benchmark multiphase fluid flow cases, such as square droplet deformation, surface wetting, and droplet migration in thermo-capillary flow are presented and validated. The results demonstrate that the developed non-local model can accurately capture the surface tension effect in multiphase fluid flow motion.


Introduction
The tension at the interface of two different states of matter is known as surface tension.The surface tension maintains the balance of two-phase immiscible fluids at the interface.Surface tension effects widely appear in multiphase flows, such as liquid droplet deformation, wetting effects of water droplets on solid planes, and thermo-capillary effects, etc.In addition, when the characteristic length scale of the system is sufficiently small, the effect of surface tension on the flow field is more prominent than the inertial effect.Due to the presence of temperature or solute concentration gradients, shear forces are exerted on the fluid surface by the tangential surface tension, which drives Marangoni flow.Therefore, it is of great practical significance to accurately simulate multiphase flow at complex interfaces.
Over the years, an enormous amount of research has been carried out to model the multiphase fluid flow motion.Based on the Navier-Stokes equation, there are two common computational fluid dynamics methods for modelling the multiphase fluid flow motion from the nanoscale to the macroscale.The first are Euler methods based on grids, such as volume of fluid (VOF) method [1].Cano-Lozano et al. [2] performed a numerical study on rising bubbles in still liquids using the VOF method to track the interface between two fluids.Hoang et al. [3] performed numerical simulations of the contact angle and wetted surface properties using the fluid volume interface tracking method and the continuum surface force method.Ma et al. [4] developed a numerical method for directly simulating the thermal Marangoni effect at the interface in two-phase incompressible fluids, and quantitatively comparing the numerical results of liquid droplet thermal capillary migration with experimental and theoretical results.Another class of numerical methods that can be used for multiphase simulations are meshless methods [5].The meshless method is a particle method, such as the smoothed particle hydrodynamics method [6] and peridynamics method [7].Morris [6] devised a technique based on smooth particle hydrodynamics for simulating two-phase flow with surface tension.This method addresses problems involving fluids of similar density and viscosity.Adami et al. [8] extended the method to higher density and viscosity ratios, using a density-weighted colour gradient formulation to reflect the asymmetric distribution of surface tension.Describing wetting phenomena, in addition to including surface tension effects at the interfaces between fluids, the interaction of fluids with solid substrates also requires the implementation of appropriate boundary conditions.Breinlinger et al. [9] extended the surface tension model using additional boundary conditions to explicitly include interactions with solid walls.Moreover, if the temperature or concentration gradient vector is tangent to the interface between the two fluids, an additional force known as the Marangoni force develops.Hopp-Hirschler et al. [10] proposed a smoothed particle hydrodynamics model of surface tension gradientdriven thermo-capillary flow based on a continuum of surface force methods, including Marangoni forces.
As an alternative approach, peridynamics is a new formulation of non-local continuum mechanics [11][12][13][14][15][16][17][18].Peridynamics has a length-scale parameter called horizon.The horizon refers to the domain of influence, which defines the extent of interaction between material points.The size of the horizon varies according to the nature of the problem [19] and influences the computational time significantly [20,21].The non-local theory of the continuum establishes a link between classical continuum mechanics and molecular dynamics.These properties make peridynamics a suitable candidate for multiscale analysis of materials.The peridynamic theory has been widely adopted in mechanical analysis [22], thermal analysis [23,24], moisture analysis [25], and thermo-mechanical coupled analysis [26][27][28][29].According to peridynamic theory, the PDDO was proposed to use integral equations representing partial differentials [30,31].Moreover, various applications of PDDO can be found in the literature [32][33][34].In addition, the governing equations are presented in integral form and do not contain any spatial derivatives.Therefore, they are always applicable regardless of discontinuity in the variable field.Recent studies have explored the application of the theory of peridynamic differential operator to fluid mechanics.Gao et al. [35] developed a non-local Lagrangian model for Newtonian fluids with low Reynolds number laminar flow, and subsequently extended the model to flow of multiphase fluids with low density ratios at low Reynolds numbers [7].Using PDDOs, Nguyen et al. [36] modelled a truly incompressible fluid based on Euler's method, and the pressure difference is no longer calculated by a weakly compressible fluid model.
So far, the development of PDDO in fluid dynamics applications has been very limited.Although the work of Gao et al. [35] proposes a PDDO for simulating multiphase fluid flow, the model does not additionally consider the case of multiphase flow wetting phenomena.Furthermore, to the best of the authors' knowledge, there are currently no peridynamic models available for simulating interface for thermo-fluid problems.Therefore, in this study, a new method for modelling multiphase flow simulations using PDDO is proposed.Thermal effects at the interface of multiphase flows and the problem of solid-liquid contact at three-phase junctions are investigated.
This paper is organised as follows.Section 2 introduces the PDDO and its properties.Section 3 provides the governing equations of motion for multiphase fluid flow.Section 4 provides the non-local model using the PDDO.In addition, the treatment for multiphase fluid contact on a solid surface is explained in this section.Section 5 presents the techniques used for numerical implementation.Section 6 presents a series of numerical case studies using the developed nonlocal multiphase fluid flow models.Finally, a conclusion is given in Sect.7.

Peridynamic differential operator
PDDO was recently proposed [22] to represent any order of partial differentials.It inherits similar ideas to peridynamics theory [11,12], which uses integral equations instead of differential equations.A function f ( ) in the two-dimensional case and up to second-order partial differentials can be computed using PDDO as in which is the initial relative position vector between material point and its family material points ′ within its horizon size H x .The relative position vector can be denoted as in which The terms 1 and 2 represent the projections of the vector with respect to x 1 and x 2 axes, and 1 , 2 are the unit (1) g 10 1 ( ) g 01  1 ( ) g 20  2 ( ) g 02 2 ( ) g 11  2 ( ) vectors.The volume of material point ′ is donated as V ′ .The peridynamic function in first-order 1 ( ) and peridy- namic function in second-order 2 ( ) can be represented as and where g 10 1 ( ),g 01 1 ( ), g 11 2 ( ), g 02 2 ( ) and g 20 2 ( ) are the peridy- namic functions, and they possess the orthogonal property as [22] in which the superscript on relative position vector components represents the power of 1 and 2 .This can also be written in the compact form as in which g P 1 P 2 N ( ) is the peridynamic function up to second order, 0 < P 1 + P 2 ≤ 2, N ≤ 2 , n i = 1, 2 , and the term P i (i = 1, 2) is the differential order with respect to variables x i (i = 1, 2) and n i P i is the Kronecker delta, so that if n i = P i then n i P i = 1 , otherwise n i P i = 0.
The peridynamic function for each material point can be constructed as [22] in which a P 1 P 2 q 1 q 2 is the coefficient matrix, (| |) is a weight function which specifies the level of non-local interaction of material points interacting with its family members.The influence of the weight function reflects the decrease of the interaction level with increase in distance between material points.A preferred weight function used in this study can be denoted as [22] (4) The unknown coefficients . in the peridynamic function in Eq. ( 8) can be expressed as As a result, using Eqs.( 8)- (10), the peridynamic function g up to second order can be represented as (10) a 01 a 02 Alternatively, these can be written in a compact form as Substituting Eqs. ( 12) and ( 9) in orthogonality property relationship given in Eq. ( 6) leads to Solving Eq. ( 13) leads to the explicit form of the coefficient matrix a P 1 P 2 q 1 q 2 , and then explicit form of peridynamic function g P 1 P 2 N can be obtained using Eq. ( 12).

The governing equations of motion for multiphase fluid flow
Fluid dynamics in multiphase fluid flow is governed by continuity equation, Navier-Stokes equation, and energy equation.

Mass conservation
The mass conservation in multiphase fluid flow motion can be described by the continuity equation as in which is density, is velocity, and t is time.

Momentum conservation
The Navier-Stokes equation in Lagrangian description has the form of [37] where represents the body force, and s is the surface tension force.The divergence of stress ∇ ⋅ can be represented as [35] (12) in which p is the hydrostatic pressure, and is the second- order unit tensor.The shear-rate tensor can be defined as [35] in which is the dynamic viscosity, and ̇ is the shear strain rate.The divergence of the shear-rate tensor can be represented as [35] (16) If the dynamic viscosity is assumed as a constant after substituting Eqs. ( 16)-( 18) into Eq.( 15), the Navier-Stokes equation can be rewritten as [35] In addition, the shear strain rate tensor ̇ can be expressed as [35] (18) Fig. 1 The transition band and unit normal vectors between two fluids If the fluid is incompressible, Eq. ( 20) can be further simplified by continuity equation.As a result, the Navier-Stokes equation for incompressible fluid flow becomes [10] The continuum surface force method [38] is adopted for modelling the surface tension force in multiphase fluid flow as the pressure jump occurred at the phase interface.The surface tension force is applied as a volumetric force and is distributed along a transition band (Fig. 1) along the interface.A weight function is acted on the transition band area to convert surface tension s into force per unit volume s [6].This can be represented as [38] where lg is the weight function for surface tension that represents the magnitude distribution of the surface tension force at the transition band, which has a peak at the interface and decays with the distance away from the interface.The weight function, lg , is further described in Sect.4.2.3.
The volumetric surface tension force s comprises the contributions from normal and tangential directions [6] which can be expressed as The normal component of surface tension force s,n represents the surface tension force due to the local curvature at the transition band area.This can be defined as in which is the temperature-dependent surface tension coefficient in N∕m , is the interface curvature, and ̂ lg is the unit normal vector at the interface between two different fluids (Fig. 1).
The temperature-dependent surface tension coefficient can be written as [10] in which d (T)  dT is the surface tension temperature coefficient [39], T is the current temperature, and 0 is the surface ten- sion coefficient at reference temperature T 0 .As the surface tension coefficient is a function of temperature, the surface tension force in tangential direction can occur due to temperature gradient and lead to Marangoni convection [10].Therefore, the term, s,t , on the right-hand side of Eq. ( 23) represents the Marangoni force and acts tangentially  to the interface, which drives the fluid from low surface tension region to high tension regions.The tangential component of the surface tension force is given as [6], where ∇ S is the surface gradient, and ∇ S can be represented as [40] in which ∇ S T is the surface temperature gradient and can be expressed as [40] As a result, the Marangoni force can be written as

Equation of state
Assuming the fluid is barotropic, an additional equation is required to uncouple the mass and momentum equations [41].In this study, the incompressible fluid flow motion is constrained by a weakly compressible equation of state, whose density is only a function of pressure.A typical equation of state is given as [42] in which 0 is the initial density, c 0 is the numerical speed of sound, and is the adiabatic exponent.p 0 is the background pressure which prevents a negative pressure field and provides tension stabilities [43].As the density changes and it is updated using continuity equation given in Eq. ( 14).The pressure field is calculated by the change between the updated density, and its initial density, 0 in the equation of state (Eq.( 30)).In two-phase fluid flow motion, the equation of state for each type of fluid flow can be expressed as [44] and where the subscripts l and g denote the denser and lighter fluids, respectively.
The adiabatic exponent defines the degree of incompressibility and pressure of fluid response to density perturbations.As density perturbations increase, a high adiabatic exponent can cause progressively large error in the pressure field.For laminar flow with low Reynolds numbers, the adiabatic exponent is taken as one ( l = g = 1) to keep the error in density and pressure proportional [42].
In weakly compressible approach, the density variation Δ in each fluid domain need to be [39]   This criterion is checked at the end of the simulation in each case.
The numerical speed of sound, c 0 , in Eq. ( 30) needs to be chosen large enough to limit the density change threshold up to 1% [39].On the other hand, the numerical stability is dependent on the time step size.The numerical speed of sound should not be too large to make the time step excessively small [42].
In this study, as the fluid domain is composed of multiphase flows with different density ratios, the numerical speed of the sound is estimated by the highest pressure change Δp in the denser fluid as [45]   For gravity-based flow, the maximum pressure variation in Eq. ( 34) is estimated by [42] where g is the force of gravity, and H is the reference depth.
For surface tension-driven flows, the pressure changes in Eq. ( 34) is approximated using Young-Laplace equation [9].The work carried out by the pressure on an interfacial area can be represented as in which dV and dA are infinitely small volume and area at the interface, respectively.
For a two-dimensional circular droplet, the pressure change at the interfacial area can be computed as where r is the radius and R is the characteristic radius of droplet curvature.The surface tension coefficient, , is cal- culated from Eq. (25).
On the other hand, the numerical speed of sound in lighter fluid is calculated as [44] (33) Δ  ≪ 1. in which c 0,l is obtained from Eq. (34).Note that by comparing the numerical speed of sound in denser in Eq. (34) and lighter fluids in Eq. (38), it can be found that the numerical speed of sound c 0,g in lighter fluid is higher than the numerical speed of sound c 0,l in denser fluid when the density ratio 0,l ∕ 0,g is significant.

Energy equation
The local form of total energy in a fluid system can be represented as [46] in which e is the internal energy per unit mass, 1  2 | | 2 is the kinetic energy per unit mass, and S is the source term.The first term on right-hand side, ∇ ⋅ q , is the net rate of heat addition due to conduction.The third term on the right-hand side, ∇ ⋅ (p ), represents the rate of doing work against pres- sure.The term ∇ ⋅ ( ⋅ ) represents the rate of doing work against viscous force, and the term (g ⋅ ) represents the rate of doing work against the body force.
Using product rule within the divergence operator ( ∇⋅ ), the rate of doing work against pressure and viscous force can be rewritten as and The mechanical energy equation can be derived from the momentum equation by multiplying velocity with momentum equation which leads to [46] Using Eq. ( 40)-Eq.( 41), mechanical energy equation in Eq. ( 42) can be rewritten as [46] In this study, the fluid is assumed to be incompressible for which the speed of fluid flow is lower than the compressible flow.Therefore, the mechanical energy can be subtracted (38) from the total energy equation and leads to the internal energy equation as [17] Defining internal energy as [17] in which C p is the specific heat capacity and substituting Eq. ( 45) into Eq.( 44) leads to [46] The heat flux q based on Fourier's Law can be represented as [46] where k is the thermal conductivity.In the case of sudden expansion or compression phenomenon, the term p∇ ⋅ rep- resents energy for the cooling or heating a fluid internally [47].Since the focus of this study is on multiphase flow and there are no significant sudden volume changes in the fluid domain, this term is omitted from the energy equation.On the other hand, the term ∶ (∇ ⊗ ) representing the motion energy is irreversibly exchanged into thermal energy, and it is considerable if the speed of the fluid is relatively high [47].As the current study focuses on the multiphase flow motion at a low Reynolds number, this term is also not considered in the energy equation.The thermal conductivity k is assumed to be a constant number as a result, the internal energy can be rewritten as

Velocity divergence in local form can be written as
The partial derivative terms in velocity divergence can be replaced by the first-order peridynamic function as Therefore, the velocity divergence in non-local form can be constructed as As a result, the non-local form of mass conservation can be represented as

Non-local form of terms in Navier-Stokes equation
As discussed in Sect.3, the Navier-Stokes equation incorporates terms for pressure gradient, viscosity, surface tension, and body forces.The non-local form of each term is expressed in this session.

Pressure gradient
The pressure gradient term ∇ ⋅ (−p ) in local form can be written as Correspondingly, using the first-order peridynamic function, its non-local form can be expressed as

Viscous force
Local form of the velocity gradient can be written as The velocity gradient matrix can be transformed into its non-local form using peridynamic function which can be expressed as Local form of Laplacian operator is defined by the divergence of the gradient as Hence, the non-local form Laplacian operator can be constructed as Subsequently, the non-local form of viscous force in a compact form can be written as (57)

Surface tension force in normal direction
The classical form of normal surface tension s,n represented in terms of surface tension coefficient , unit normal vector ̂ lg , weight function lg and interface curvature is denoted in Eq. (24).To construct the normal surface tension in a non-local form, it is first necessary to construct the fluid interface normal vector, ̂ lg , and curvature, , in a non-local form.
According to the continuous surface force method [38], as the colour function has a unit jump at the interface, it can be used to identify and track the position of the interface.The normal vectors between fluid and gas interface can be represented by the gradient of the colour function, ∇c lg ( ) , as where c lg ( ) is the colour function at material point .The difference of colour function between a pair of material points can be represented by In numerical simulations, if there is a large density difference between two fluids at the interface, a weighted-density approach is used as an alternative method to determine the difference of colour function [8], i.e.
where is the density of material point located at , and ′ is density of material point located at ′ .Following the approach purposed by Morris [6],the unit normal vector between two fluid domains, ̂ lg ( ), as shown in Fig. 1 then can be formulated using the gradient of colour function in Eq. ( 60) as ( 60) if and � areinsamef luiddomain otherwise . (62) if and � areinsamef luiddomain otherwise , In addition, the weight function for surface tension in Eq. ( 24) in continuum surface force method is taken as the magnitude of the gradient of the colour function [6] and its non-local form is provided as [7] the weight function for surface tension is used to convert the surface tension force into volumetric force, and distribute the force along the fluid interface transition band area.
Figure 2 shows a droplet surrounded by a gaseous fluid and resting on a solid surface.The unit normal vectors, ̂ lg ( ) , from Eq. (63) represent the normal direction of the interface between liquid and gaseous fluid, and they can be accurately computed when material points in each fluid domain fully interact with their family material points in fluid domain.However, at the triple line region, where the liquid-gas interface meets the solid-liquid interface in Fig. 2, the material points in the fluid domain close to the solid wall do not have enough family material points to contribute to the integral equation, unit normal vectors between fluids calculated according to Eq. ( 63) can be corrupted.Moreover, as the curvature is calculated from the divergence of unit normal vectors at the interface, these corrupted unit normal vectors result in erroneous curvature calculations.
Therefore, consideration is required when computing unit normal vectors at the triple line region.In this study, the corrupted unit normal vectors at this region are corrected through a normal prescription scheme [9].
The unit normal vectors at the triple line region, ̂ lg,cor ( ) , as shown in Fig. 3 can be prescribed as [9] in which ̂ t is the projection of unit normal vector, ̂ lg ( ) , between the denser fluid and the lighter fluid on the solid-fluid interface, ̂ sf ( ) is the unit normal vector between (63) ) ̂ lg,cor ( ) = ̂ t ( )sin eq − ̂ sf ( )cos eq solid phase and fluid phase, and eq is the equilibrium contact angle.
When a droplet is in contact with a solid surface, the balance between the adhesive and cohesive forces in the droplet forms the equilibrium contact angle eq .Equation (65) prescribes unit normal vectors at the triple line region point in the direction of the interface normal that forms the equilibrium contact angle eq .If a droplet comes into contact with a solid surface and forms an instantaneous contact angle not equal to the equilibrium contact angle eq , the curvature obtained by the divergence of the prescribed unit normal vectors ̂ lg,cor ( ) will drive the interface to move until the droplet forms an equilibrium contact angle with the solid interface [9].
Following a similar approach introduced in continuum surface force method, the unit normal vector between solid phase and fluid phase, ̂ sf ( ), in Eq. ( 65), can be formulated as [39] in which ∇c sf ( ) is the gradient of the colour function between fluid and solid phase.Similarly, its non-local form can be represented as where c sf ( ) is the colour function at material point for distinguishing fluid and solid phase.The difference of colour function between a pair of material points can be defined as Consequently, the projection of unit normal vector between liquid and gaseous fluid on the solid-fluid interface, ̂ t ( ) , in Eq. ( 65) can be computed as [9]   (66) if and � areinsamephase otherwise .where ̂ sf ( ) is provided in Eq. ( 66) and ̂ lg ( ) is provided in Eq. ( 63).
The ordinary computed unit normal vectors ̂ lg ( ) from Eq. ( 63) and prescribed unit normal vectors ̂ lg,cor ( ) at the triple line region from Eq. ( 65) are represented in Fig. 3. Since the calculation of ̂ lg,cor ( ) depends on the equilibrium con- tact angle eq , when current contact angle is greatly different from equilibrium contact angle eq , a sharp transition can be observed between ̂ lg ( ) and ̂ lg,cor ( ) .This sharp transition will cause a discontinuity when computing the curvature.As a result, a smoothed unit normal correction scheme is implemented here to ensure a smooth transition from prescribed unit normal vectors to ordinary computed unit normal vectors, the smoothed interface unit normal vectors can be obtained as [9] where ̂ lg,cor ( ) can be calculated from Eq. (65) and ̂ lg ( ) can be calculated from Eq. (63).The parameter f w,x is a transition function and determines the influence of the prescribed normal vector at the triple line region which depends on the distance to the wall.It is provided as [35] in which d w is the distance between fluid material points and the solid-fluid interface.
The d max in Eq. ( 71) is denoted as the maximum smooth distance from the wall, in this work, it is taken as 2Δx , where Δx is the spacing between material point.The schematic dia- gram of the smoothed interface unit normal vectors calculated from Eq. ( 70) is shown in Fig. 4. Section 6.2 gives a numerical example of the effectiveness of this smoothing correction method. (69) Subsequently, the surface curvature in Eq. ( 24) then can be calculated via the divergence of the smoothed interface unit normal vector, i.e.
As suggested in Morris' work [6], curvature directly computed from Eq. ( 72) can lead to errors at the edges of the transition region, as the smoothed interface unit normal vectors are relatively small and can have erroneous directions when they are away from the interface.Therefore, the surface curvature cannot be approximated accurately.This problem can be addressed by introducing selection criteria to determine if a 'reliable' normal vector can be obtained for divergence computation.A function at each material point is used to distinguish 'reliable' normal vectors that can contribute to the curvature approximation in Eq. (72) [6], i.e. and where  ≪ 1 is a user-defined tolerance [6].As the unit nor- mal vectors below the tolerance are not contributed to the curvature computation, an intermediate curvature estimation needs to be used to sum over reliable normal vectors [6].As a result, the curvature in Eq. ( 71) can be recomputed as On the other side, considering material points at the edge of the phase transition region, whose family material points are within the horizon but outside the transition region, the interface unit normal vectors at these family material points (72) are zeros.As a result, a correction factor is used to consider the truncated material points, and it can be represented as [7] in which ω( ) is the weight function as represented in Eq. ( 9).Based on Eq. ( 70), the smoothed interface unit normal vectors,̂ * lg ( ) obtained for material points at each fluid domain always point from themselves to the fluid interface.Considering a pair of material points which is in one fluid domain and ′ is in the other fluid domain, the direction of their unit normal vectors will opposite with each other.Therefore, a phase normal coefficient ′ is added to reverse the unit normal vector direction if it is opposite from the unit normal vector at material point [44].Hence, the interface curvature in this study is computed as with Finally, the non-local form of the normal surface tension can be expressed as in which is given by Eq. ( 25), * * ( ) is the curvature as described in Eq. ( 77), ̂ * lg ( ) is the smoothed unit normal vector between liquid-gas interface as described in Eq. ( 70) and lg is weight function for surface tension as described in Eq. (64).

Marangoni force
The classical form of the Marangoni force is given in Eq. (29).
To evaluate the non-local form of Marangoni force, the temperature gradient ∇T in non-local form needs to be developed.The temperature gradient ∇T in local form is given as Therefore, the non-local form of the temperature gradient using the peridynamic function can be expressed as As a result, the non-local form of Marangoni force can be represented as (79) s,n = * * ( )̂ * lg ( ) lg ( ) Fig. 5 Location of material point n i with its family members at t = t n , and the updated location at t = t n+1 where ∇T( ) is calculated from Eq. (81), ̂ * lg ( ) is calculated from Eq. (70), and lg is calculated from Eq. (64).

Thermal model
The classical form of the energy equation is given in Eq. (48).The local form of the divergence of the heat flux in the equation can be written as The second-order partial derivative for temperature can be represented by second-order peridynamic function as As a result, the Lagrangian description of internal energy equation in non-local form can be written as

Numerical implementation
The non-local form of governing equations is solved numerically.The fluid domain is discretised into a series of material points, and each material point carries information such as material density, viscosity, pressure, velocity, displacement, and temperature.Since PDDO inherits analogous concepts with peridynamics theory, the long-range force is considered in the simulation domain representing a material point interacting with a series of family material points within a horizon.As shown in Fig. 5, n i represents the current coordinate of the material point i at time t = t n , and material point i interacts with family material points j at n j within a range of .With the updated Lagrangian description, the location of the material points changes at every time step.Therefore, a material point may interact with different family material points at time t = t n+1 and the peridynamics function needs to be reconstructed based on the updated configuration.

Discretised form of governing equations using PDDO
To ensure the mass conservation is maintained at the fluid domain, each material point in the simulation domain is assigned with an initial mass after the simulation domain is discretised.The mass at each material point remains same during simulation, as the density field of the material point is updated by continuity equation, the volume of material point can be updated correspondingly at each time step as where m i is the initial mass of material point i , and n i is the density at the current time step.
The discretised form of the peridynamic function up to second order is given as and in which the relative distance vector between a pair of material points i and j can be represented as where superscript n represents the current time step.
As introduced in Eq. ( 12), peridynamic function up to second order is represented as in which coefficient matrix in Eq. ( 90) is obtained by numerically solving Eq. ( 13) as with and where V n j is the volume of family material point j , N i repre- sents the number of family material points of material point i.
As suggested by Madenci et al. [31], horizon size in the numerical simulation is chosen as maximum order of differentiation plus one.Since, the maximum order of differentiation is two, the horizon size is chosen as 3Δx , in which Δx is spacing between material points at the initial configuration. (90) The discretised form of the continuity equation is given as in which n+1 i is the updated density, Δt is the time step size, and term n i and n j represent the velocity field of material point i and j , respectively.
In addition, the pressure field of the material point i at time t = t n is computed by equation of state as expressed in Eq. ( 30) as (95 in which p n i is the pressure at the current time step and p 0,i is the background pressure.The background pressure in this study is estimated as Depending on the type of the fluid motion, the numerical speed of sound c 0,i of material points in denser and lighter fluid phase is estimated using Eqs.(34) and (38), respectively.
After the pressure field at each material point is computed by Eq. ( 96), the discretised form of the pressure gradient in momentum equation then can be represented as (97) p 0,i = 0.05 × 0,i c 2 0,i i .

The discretised form of the viscous force in momentum equation can be computed as
The discretised form of the normal surface tension force in momentum equation is calculated from Eq. ( 79) as (98) s,n in which is the temperature-dependent surface tension coefficient and is computed from Eq. ( 25).The term c lg n j − c lg n i represents the difference of colour index between a pair of material points, and i is the phase normal coefficient.̂ * n lg i is the smoothed interface unit normal, at time t = t n this can be expressed as where f w,i is the transition function.The ordinary computed unit normal vector, ̂ lgi n ,is calculated as and the prescribed unit normal vector at triple line region, ̂ lg,cori n ,is calculated as with and (101) in which the term c sf n j − c sf n i represents the difference of colour index for distinguishing fluid and solid phase between a pair of material points.eq is a user predefined equilibrium contact angle before simulation.
N is the unit normal index and is calculated as The discretised form of the Marangoni force in momentum equation then can be calculated from Eq. (82) as in which d (T)  dT is the surface tension temperature coefficient as shown in Eq. ( 25).
As a result, the discretised form of momentum equation can be represented as where n+1 i is the updated acceleration at t = t n+1 and i is the body force density.When the density ratio between two fluids is relatively large, the discontinuity of the fluid density and viscosity field will present at the transition band of the Fig. 6 Schematic drawing of interface between two fluids and the boundary material points interface.As PDDO is a non-local approach, the discontinuity will cause stability issues 0. In order to reduce numerical oscillations and prevent material points' penetration during the simulation, the density and viscosity coefficient can be smoothed by harmonic means as 0 and The harmonic means for density and viscosity treatment is validated for multiphase flow fluid motion at low Reynold's number.In the case of the flow motion with high Reynold's number, additional numerical treatments are required [7].In this study, only the flow with low Reynold's number is investigated.
The discretised form of the internal energy equation can be represented as Similarly, to prevent numerical oscillation in energy equation and to have a smooth transition of the heat conductivity coefficient at the interface region, the coefficient k n i in Eq. ( 111) is also smoothed by harmonic means as

Boundary conditions
As shown in Fig. 6, the boundary conditions are implemented through fictitious layers in the numerical simulation, which is widely adopted in peridynamic studies [24].Three variables are considered in boundaries, which are velocity, pressure, and temperature.
Slip or no-slip velocity boundary conditions are used at wall boundary and implemented by material points in fictitious layers.The velocities of the material points in fictitious layers are computed based on the velocities of material points in the fluid domain.Different methods for implementing these boundary conditions are used in numerical simulations, such as mirroring material points at boundary approach as suggested in Ref. [24], so that the variable at a material point in the fictitious layers is mirrored by a material point at the fluid domain.Since the fluid particles are (109) moving during the simulation, instead of mirroring the moving fluid particles, a simplified boundary implementation using weighted average approach is used to keep the material points in the fictitious region at fixed locations [7].
For slip or no-slip velocity boundary conditions, the velocity of material point i in a fictitious layer at the current time step n is calculated as with in which wall is the wall velocity, the subscripts j and k represent the family material points of i in fluid 1 and fluid 2, respectively.The superscript N 1 and N 2 represents the material point i in the fictitious layer has N 1 family material points in fluid 1 and N 2 family material points in fluid 2, respectively.The velocity of family material points of i in fluid 1 and fluid 2 at current time step n are denoted by n j and n k , respectively.The weight function in Eq. ( 113) between a pair of material points is computed as [43] where the term e represents Euler number.
Without considering the gravity force, the pressure of material points in fictitious layers are calculated from the pressure of the material points in fluid domain as [7] in which p n j and p n k are the pressure of family material points of i in fluid 1 and fluid 2 at current time step n.
Dirichlet and Neuman temperature boundary conditions are also implemented by material points in fictitious layers.With the analogous ideas of computing velocity and pressure for material points at fictitious layers, the temperature of material point i at current time step n is computed from its family material points temperature T n j and T n k in fluid 1 and fluid 2 as (113 with in which T wall is the wall temperature.

Time stepping scheme
The momentum equation is integrated explicitly in time using Velocity Verlet scheme as 0, and the displacement field at each material point in fluid domain is updated by As a result, the updated location of the material points can be found as in which the superscript 0 represents initial coordinate of the material points.
To maintain the numerical stability in time integration, the time step size Δt is constrained by the Courant-Frie- drichs-Lewy (CFL) condition [48].The CFL condition in time integration is based on several conditions.The time step size for numerical speed of sound condition is [49] in which c max is the maximum numerical speed of sound among all phases.The estimation of numerical speed of sound in each fluid phase is provided in Eq. ( 34) and Eq. ( 38), and max is the maximum velocity in the simulation domain.
The time step size for viscous condition is constrained as [49] where 0 max is the maximum kinematic viscosity among all phases. (117) (121) The time step size for body force condition is applied as [49] in which g is the gravity acceleration.
The time step size for surface tension condition is implemented as [8] In addition, the time step size in thermal analysis is restricted as [37] In processing the numerical time integration, the time step size is chosen as the minimum of above criteria as

Material points shifting technique
Since distorted material points induce stability issues in processing the numerical integration [50], the position shifting technique is applied at material points at each time step in fluid domain to avoid clustering problems.The application of position shifting technique is introduced in Ref. [7] for PDDO.At each time step, the displacement for each material point in fluid domain n+1 i is corrected by a shifted distance Δ n+1 i * , which is represented as The shifting distance is defined as [7] in which i is the displacement shifting vector.To ensure the shifted distance can sufficiently prevent the instability and do not cause accuracy issues, the constant C is typically taken between 0.01 and 0.1.The shifting magnitude MPST is set as [50] In addition, the displacement shifting vector in Eq. ( 129) is provided as [7] (124) Δt b ≤ 0.25 The summation of distance vectors in Eq. ( 131) describes the anisotropy of the spacing between material points, and 2 is used as a weight function to evaluate the influence from material point j [21].The averaged material points spacing ξ i in Eq. ( 131) is defined as [50]

Moving least square method
In Lagrangian method, the position of material points is tracked and updated at each time step.When position of material points in the fluid domain changes continuously, the number of family material points may decrease.In this case, the calculated density may be smaller than normal.Therefore, the equation of state predicts wrong pressure values, leading to a gradual deterioration of the entire field [51].To avoid mass conservation issues, and oscillations at the density, pressure, and velocity of material points at fluid domain are smoothed using moving least square method [51].
The velocity field in the fluid domain is smoothed as [7] (131 The pressure field in the fluid domain is smoothed as [7] The density field in the fluid domain is smoothed as [7] with It is worth noting that in Eq. ( 135), the smoothed density n i smoothed is obtained from the density of family material points n j as provided in Eq. ( 136) [7].The weight function, MLS , for smoothing variables in above equations at fluid domain is provided as [43] (133) In addition, to reduce the numerical computation time, the velocity, pressure and density field only being corrected over a period of time steps using moving least square method, and this is optional in numerical simulation.For benchmark cases moving least square method is used for dynamic cases (Sects.6.1, 6.3 and 6.5), however this method is not needed for static cases (Sects.6.2, 6.4).( 138) .

Numerical simulations
In this section, five numerical cases are considered using the developed PDDO for modelling the surface tension forces in multiphase fluid flow motion.First, a two-dimensional square droplet deformation is studied to examine the nonlocal form of surface tension force model in normal direction.Second, when fluid flow is in contact with a solid surface, the difference between unit normal vectors at the triple line region before and after using the unit normal vector prescription scheme are compared through a static droplet wetting case.Third, simulation of droplet contact angle development on a solid surface is performed to show the effect of prescribed normal vectors at the triple line region.Afterwards, capillary stresses tangential to the interface are computed under a heat conduction phenomenon to validate the non-local form of the Marangoni force formulation.Finally, a two-dimensional droplet migration in a thermocapillary flow is presented to test the combination of the surface tension forces in the normal direction and the tangential direction.The predicted migration velocity of the circular droplet in the thermo-capillary flow is compared with the volume of fluid method.

Droplets deformation
In the first case, a two-dimensional square droplet deformation is conducted to validate the surface tension force model using PDDO.As shown in Fig. 7a, the square droplet with dimensions of 0.6m × 0.6m is filled with fluid 2 and surrounded by fluid 1.The fluid domain has a box of size L = W = 1m .The density for fluid 1 and fluid 2 are speci- fied as 1 = 10kg∕m 3 and 2 = 1kg∕m 3 , respectively.Fluid 1 has a viscosity coefficient of 1 = 1Pas while fluid 2 has 2 = 0.2Pas .The surface tension coefficient between fluid 1 and fluid 2 is independent with the temperature and chosen as = 1N∕m.
The fluid is initially at rest for which the initial condition can be illustrated as No-slip boundary conditions are applied at four edges of the fluid 1 as In processing the numerical simulation, as shown in Fig. 7b, three layers for fictitious material points are wrapped along the edges of the fluid domain.No-slip boundary conditions are implemented on these fictitious material points using Eqs.( 113)-( 114) with wall = 0m∕s .The pres- sure field for fictitious material points is computed using Eq. ( 116).The deformation of the square droplet is driven by surface tension force, which transforms square droplet in a circular shape at equilibrium state.With incompressibility hypothesis, the radius of the final circular droplet can be estimated as Therefore, the pressure changes Δp between fluid 1 and fluid 2 can be estimated using Young-Laplace equation as As introduced in Sect.3.3, the incompressible fluid motion is constrained by a weakly compressible equation of state, and the numerical speed of sound in equation of state is computed based on the pressure change.As a result, the numerical speed of sound for fluid 1 can be estimated using Eq.(34), and the numerical speed of sound for fluid 2 can be obtained by Eq. (38).In this case, numerical speed As shown in Fig. 7b, the fluid domain is discretised with a uniform spacing of Δx = 0.0125m .The horizon size is selected as = 3Δx .Simulation is processed for a total time of t = 2.56s and the time step size is set as Δt = 8 × 10 −5 s .The displacement field of the material points in fluid domain is obtained by velocity Verlet scheme [52].
In addition, the moving least square method introduced in Sect.5.5 is used to correct the density, pressure, and velocity field at fluid domain at every 20 time steps.To obtain a smooth distribution of material points in the fluid domain, the material points shifting technique is utilised in every time step of the simulation to smooth the displacement field.The constant C is taken as 0.01 in this case.
Figure 8 shows the snapshots of the droplet transformation from square shape to circular shape.As can be observed from Fig. 8f, at the final state, the droplet has an average radius of 0.0338m , which is close to the estimated value given in Eq. ( 142).The time history of average pressure difference between fluid 1 and fluid 2 is presented in Fig. 9a.As can be observed from the figure, the pressure difference reaches to an equilibrium state after t = 1.5s.The pressure change Δp between fluid 1 and fluid 2 pre- dicted by PDDO shows a close agreement with the analytical solution as computed by Eq. ( 143).The final pressure profile of the fluid domain at t = 2.56s is presented in Fig. 9b.The square droplet transforms into a circular shape with a smooth material point distribution, and the pressure difference between droplet and surrounding fluids match with the analytical value from Young-Laplace equation.It can be concluded that the current surface tension model utilising the PDDO can accurately capture the surface tension effect in the normal direction in multiphase flow.

Unit normal vector correction at triple line region
The first case demonstrates the effectiveness of surface tension force model in normal direction.However, this only validates the case when fluid 2 is fully surrounded by fluid 1.When both fluids are in contact with a solid interface, additional processing of the unit normal vectors between the two fluids at the triple line region is required to prevent curvature errors.As shown in Fig. 10, this section uses cases of droplets forming two different contact angles on a solid surface to demonstrate the difference between the unit normal vectors obtained before and after normal prescription scheme at triple line region.
As shown in Fig. 11, two denser fluid droplets with radius r = 0.0125m are surrounded by lighter fluid, and they are both being placed in a rectangular box.The box has a size of 0.1 m in length and 0.05 m in width.The density for denser fluid is 1 = 1000kg∕m 3 and for lighter fluid is 2 = 1.2kg∕m 3 .Two denser fluid droplets contact with the solid surface and form a contact angle of 30° and 90°, respectively.The spacing between material point is set as Δx = 0.0013m.
The unit normal vectors between the denser and lighter fluid domains, calculated according to Eq. ( 63), are shown in Fig. 12.It can be observed that in the case of a 30° contact angle, the unit normal vectors ̂ lg ( ) at the triple line region in the lighter fluid domain is pointed towards the wrong direction.The issues are also reflected in the 90° contact angle case.This is because, at the triple line region, there are not enough family material points to contribute to the integral equation when computing the unit normal for the lighter fluid.Therefore, incorrect unit normal vectors will result in incorrect curvature calculations and affect the surface tension force modelling.
Considering that the same droplet forms contact angles of 30° and 90°, respectively, the prescription normal vectors ̂ lg,cor ( ) at triple line region computed based on Eq. ( 65), and interface unit normal vectors ̂ lg ( ) not at the triple line region obtained based on Eq. ( 63) are presented in Fig. 13.
As shown in Fig. 13, while the incorrect unit normal vectors at the triple line region are resolved in both contact angle cases, a sharper transition is observed between the two types of unit normal vectors.Therefore, calculating the curvature of the interface remains problematic.
To prevent discontinuity of the unit normal vectors not at the triple line region, the unit normal vectors at the triple line region are smoothed according to Eq. (70).The unit normal vectors between denser and lighter fluid domains after smoothing are presented in Fig. 14.As can be observed, the unit normal vectors close to the triple line region comply with the corrected unit normal vector ̂ lg,cor ( ) while far from the triple line region unit normal vector is computed from ̂ lg ( ) .When the instantaneous contact angle formed by a droplet contacting a solid surface is not equal to the equilibrium contact angle eq , the curvature creates a force to move the triplet until the equilibrium contact angle is reached.

Static contact angle development
After demonstrating the influence of treatment for unit normal vectors at the triple line region, in this section, a two-dimensional liquid droplet deformation on a rigid wall is investigated to study the characteristics of droplet wetting in different stages.As shown in Fig. 15a, a rectangular liquid droplet with a dimension of 2.25 × 10 −2 m 2 in length and 1.25 × 10 −2 m 2 in width is placed in a rectangular box with a dimension of 0.1 × 0.05m 2 .The liquid droplet on the wall is surrounded by gas fluids.The density and viscosity  Fig. 16 Evolution of liquid droplet wetting on a solid surface for various contact angles coefficient of the gas phase are specified as 1 = 1.2kg∕m 3 and 1 = 1 × 10 −3 Pas , respectively.The density and vis- cosity coefficient for liquid is set as 2 = 1000kg∕m 3 and 2 = 3.16 × 10 −5 Pas , respectively.To focus on the physi- cal characteristics of droplet wetting, the surface tension coefficient between liquid and gas fluids is chosen as = 0.07N∕m , which is independent of the temperature variation.
The fluid gravitational acceleration is disregarded in this example.The initial condition of the fluid is provided as Figure 15b shows the numerical simulation domain.The domain is discretised with a uniform spacing of Δx = 1.0 × 10 −3 m .The horizon size is chosen as = 3Δx .The material points in the fluid 1 and fluid 2 are presented in dark red and red, respectively.The gas and liquid fluid are wrapped by three layers of fictitious material points while the rigid solid wall boundary is presented in orange colour.The no-slip boundary conditions are implemented as Therefore, the velocities at fictitious material points are computed based on Eq. ( 113), (114) with v wall = 0 .The pres- sure field at fictitious material points is calculated according to Eq. ( 116).
The time step size is set as Δt = 5 × 10 −5 s with a total simulation time of t = 1s .The displacement field of the material points in fluid domain is acquired using velocity Verlet scheme.In this case, the moving least squares technique is applied in the simulation every 20 steps to smooth the velocity, pressure, and density fields.In addition, the material points shifting technology is used in the simulation, and the constant C in Eq. ( 129) is taken as 0.01 to ensure the smooth distribution of material points at each time step.Figure 16 presents snapshots of the two-dimensional liquid droplet formation on the solid wall for three different cases eq = 60 • (hydrophilic wetting), eq = 90 • and eq = 150 • (hydrophobic wetting).The droplet initially stays as a rectangular shape.The unit normal vectors between the gas and liquid phases at the triple line region are corrected and smoothed based on the prescribed equilibrium contact angle using Eq.(70).The unit normal vectors form a smooth curvature at the interface between gas and liquid.The curvature obtained from the corrected unit normal vectors induced a force along the fluid interface to deform the droplet until the prescribed equilibrium contact angle is reached.

Capillary stress tangential to interface
After investigating the surface tension force in normal direction, in this section, the developed non-local Marangoni force formulation is examined by considering a heat conduction test case.As shown in Fig. 17a, two-layered fluids are placed in a square simulation domain with 5.75 mm in length and width.The heat conduction model and  29).In addition, the surface tension coefficient is a temperature-dependent property which is given in Eq. ( 25), in which the surface tension temperature coefficient is chosen to be d (T)∕dT = 0.002N∕K.
The time step size is chosen as Δt = 1 × 10 −5 s , and the total simulation time is t = 0.1s.
At the initial state, the temperature distribution of the fluid domain is set to zero as The boundary conditions are implemented using the fictitious layers as shown in Fig. 17b and defined   To accurately capture the Marangoni force distribution at the interface, the heat conduction model is first examined with a mesh size of Δx = 9.0 × 10 −5 m .The horizon size is chosen as = 3Δx .The temperature distribution pre- dicted using PDDO at t = 0.1s is compared with ANSYS and presented in Fig. 18.A good agreement is observed between the two methods, which validates the heat conduction model.In addition, as can be computed from the temperature field and the analytical solution, a temperature gradient of ∇T = 200K∕m is distributed along the width of the simulation domain.As a result, Marangoni force is developed vertically at the interface.The theoretical Marangoni force distributed along the interface is computed as s,t = ∇ S = 0.4N∕m 2 .The Marangoni force is examined under different spacings between material points, in which the mesh size is chosen as Δx = 1.8 × 10 −4 m, Δx = 9.0 × 10 −5 m and Δx = 4.5 × 10 −5 m .Figure 19 presents Marangoni force vectors distributed along the interface using a mesh size of Δx = 9.0 × 10 −5 m .The profile of the Marangoni force vectors is perpendicular to the interface.As the continuum surface force method is adopted in the model, the volume Marangoni force is smoothed and distributed symmetrically along the transition band of the fluid interface.
Figure 20 shows the Marangoni force magnitude distribution profile in the horizontal direction at the centre of the simulation domain using various mesh sizes.It can be noticed that the maximum Marangoni force decreases as the material point spacing increases.The continuum surface force model handles the local Marangoni forces at the fluid interface by applying them to material points in the transition zone between the two fluids.As the horizon size is taken as = 3Δx , material points adjacent to the interface but beyond a horizon size from the interface do not capture surface tension forces.As a result, only three layers of material points on either side of the interface capture the Marangoni force regardless of the spacing chosen.The magnitude distribution of the Marangoni force at these material points is governed by a weighting function that controls the decay of the Marangoni force with increasing distance from the interface.
To compare the numerically computed Marangoni force with the analytical solution, Fig. 21 presents the integral of the Marangoni force distributed at the three-layer material points along the interface for different spacings between material points.As can be observed, the calculated Marangoni force has a good agreement with the analytical solution despite the spacing between material points being different.Therefore, the presented Marangoni force formulation can accurately capture the Marangoni force due to the temperature gradient.

Two-dimensional droplet migration in thermo-capillary flow
After validating the surface tension formulation in normal direction, the Marangoni force formulation, and heat conduction model, in this case, surface tension in normal and tangential directions are combined to investigate the motion of a droplet in thermo-capillary flow.Thermocapillary flow motion was studied in past decades experimentally and numerically [53][54][55].In this case, the flow motion is simulated using PDDO.As shown in Fig. 22a, a    The Neumann boundary condition is applied on the lateral edges of the square box.As represented in Fig. 22b, the fluid domain is wrapped with three layers of material points for implementing the temperature, pressure, and velocity boundary conditions.No-slip and free slip boundary conditions are implemented on these fictitious material points using Eqs.( 113), (114) and the pressure field at fictitious material points is calculated using Eq.(116).In addition, Dirichlet and Neumann temperature boundary condition for fluid domain are implemented using Eqs.(117), (118) for material points at fictitious region.
The fluid domain is discretised with a uniform spacing of Δx = 9 × 10 −5 m and the horizon size is taken as = 3Δx .The simulation is conducted with a total time of t = 0.12s and the time step size is chosen as Δt = 1 × 10 −5 s .Material points shifting technique is implemented, and the constant C is taken as 0.01.In addition, the density, pressure, and velocity field are smoothed by moving least square method every 20 time steps.
The interface between the circular droplet and surrounding fluid is subjected to the surface tension force in normal direction.The pressure difference between two fluids maintains the equilibrium of circular shape.The pressure field of the fluid domain at t = 0.12 s is shown in Fig. 23.The pressure difference between two fluids interface can be analytically verified by Young-Laplace equation, in which ΔP = 6.944N∕m 2 .As can be seen from the figure, a good agreement with the analytical value is observed for the pressure difference between fluid 1 and fluid 2.
The dimensionless parameters Reynolds number, Marangoni number, and capillary number are used to characterise thermo-capillary migration.To compare the result with the Volume of fluid (VOF) method, the dimensionless parameters are taken as the same as the case presented in the VOF method [4].The dimensionless parameters in this case are set as and in which U r is the characteristic velocity, and it is being com- puted as The migration of the circular droplet at different times is provided in Fig. 24.To visualise the migration of a droplet, the lowest location of the circular droplet at the initial stage is used as a reference location and presented as a dash line in the figure .The ratio between square box and the radius of the droplet is L∕R = 4 .The migration velocity of the circular droplet using VOF method at this ratio is provided in Ref. [4].The comparison between PDDO method and VOF method for the time evolution of droplet migration velocity is presented in Fig. 25.The velocity is non-dimensionalised as U * = U∕U r , and the time is non-dimensionalised as t * = t∕t r with t r = R∕U r .As can be observed from the figure, the migration velocity predicted using the proposed method has a good agreement with the velocity predicted using the VOF method.
The temperature and velocity field of the fluid domain at t = 0.12s are provided in Figs. 26 and 27, respectively.The uneven temperature distribution in Fig. 26

Conclusion
This study presents a new non-local surface tension model in multiphase fluid flow through the PDDO.The model considers surface tension in the normal direction, Marangoni forces, and surface wetting.The governing equations of multiphase flow motion are represented using the PDDO.The non-local form of normal surface tension is verified by simulating the deformation of a square droplet.Subsequently, this work explains the handling of the unit normal vector at the triple line region using the static and dynamic behaviour of droplet wetting on solid surfaces.Furthermore, this work validates the accuracy of the newly developed nonlocal form of the Marangoni force formulation via a heat conduction model.Finally, the normal surface tension and Marangoni forces formulations are simultaneously examined in the multiphase flow by simulating the migration of droplets in thermal capillary flow.A good agreement is observed by analysing the droplet migration speed and comparing the results with existing methods.The current study is limited for flows with low Reynolds number.However, it can be extended for flows with higher Reynolds number in future studies.

Fig. 2
Fig. 2 Two fluids come into contact at a solid surface, and the triple line region at the point of contact

Fig. 3
Fig.3Unit normal vector ̂ lg,cor ( ) at the triple line region and ordinary computed unit normal vector ̂ lg ( )

Fig. 4
Fig. 4 Corrected unit normal vectors distributed along the fluid interface

Fig. 7
Fig. 7 Investigation of square droplet deformation a geometry, b PDDO discretisation

Fig. 9 Fig. 10 Fig. 11
Fig. 9 Pressure profile a time history of pressure difference between the droplet and surrounding fluid, b distribution at t = 2.56s

Fig. 12 Fig. 13
Fig.12 The unit normal vectors between fluid domains before correction a 30° contact angle, b 90° contact angle

Fig. 14 Fig. 15
Fig. 14 Smoothed unit normal vectors at fluid interface a 30° contact angle, b 90° contact angle

Fig. 17
Fig. 17 Investigation of Marangoni force at the interface between two fluids a geometry, b PDDO discretisation as (146)T = 0 at t = 0.

Fig. 18
Fig. 18 Temperature distribution of the simulation domain at t = 0.1s (a), PDDO (b), ANSYS (c) comparison of PDDO and ANSYS results along central y-axis

Fig. 19 Fig. 20 Fig. 21
Fig. 19 Marangoni force vectors distributed along the interface with Δx = 9.0 × 10 −5 m circular droplet with radius R = 0.00144m is initially located at the centre of the simulation domain, and it is filled with fluid 1.The density and viscosity coefficient for fluid 1 are set as 1 = 250kg∕m 3 and 1 = 0.012Pas ., respectively.The droplet is surrounded by fluid 2 in a square box with a dimension of L = W = 0.00576m .The density and viscos- ity coefficient for fluid 2 are specified as 2 = 500kg∕m 3 and 2 = 0.024Pas .In heat conduction model, the specific heat capacity for fluid 1 and fluid 2 is c p1 = 0.5 × 10 −4 J∕kgK and c p2 = 1.0 × 10 −4 J∕kgK , respectively.The heat con- duction coefficient for fluid 1 and fluid 2 is chosen ask 1 = 1.2 × 10 −6 W∕mK a n d k 2 = 2.4 × 10 −6 W∕mK , respectively.As the incompressible fluid motion is constrained by a weakly compressible equation of state, the numerical speed of sound in the equation of state for each fluid domain is set as c 1 = 1.666m∕s and c 2 = 1.178m∕s , and the material constants are 1 = 2 = 1 .The surface tension

2 =
causes a temperature gradient across the fluid domain.As a result, tangential forces are created at the interface between fluid 1 and fluid 2. Combining with the viscosity of the fluid, the droplet is pushed to move along the thermal gradient upwards, and the (152) Re = 2 RU r

Fig. 27
Fig. 27 Velocity profile of the thermo-capillary migration of the droplet at t = 0.12s