An Eulerian–Lagrangian mixed discrete least squares meshfree method for incompressible multiphase flow problems

Applied Mathematical Modelling 76 (2019) 193–224

Contents lists available at ScienceDirect

Applied Mathematical Modelling journal homepage: www.elsevier.com/locate/apm

An Eulerian–Lagrangian mixed discrete least squares meshfree method for incompressible multiphase flow problems S. Faraji Gargari a, M. Kolahdoozan a,∗, M.H. Afshar b, S. Dabiri c a

School of Civil Engineering, Amirkabir University of Technology (Tehran Polytechnic), 424 Hafez Ave., Tehran, Iran School of Civil Engineering, Iran University of Science and Technology, Narmak, Tehran, Iran c Department of Agricultural and Biological Engineering and School of Mechanical Engineering, Purdue University, 585 Purdue Mall, West Lafayette, IN 47907, USA b

a r t i c l e

i n f o

Article history: Received 17 July 2018 Revised 25 April 2019 Accepted 5 June 2019 Available online 12 June 2019 Keyword: Meshfree methods Mixed discrete least squares meshfree (MDLSM) method Discrete least squares meshfree (DLSM) method Multiphase problems Moving least squares (MLS)

∗

a b s t r a c t Mixed discrete least squares meshfree (MDLSM) method has been developed as a truly meshfree method and successfully used to solve single-phase flow problems. In the MDLSM, a residual functional is minimized in terms of the nodal unknown parameters leading to a set of positive-definite system of algebraic equations. The functional is defined using a least square summation of the residual of the governing partial differential equations and its boundary conditions at all nodal points discretizing the computational domain. Unlike the discrete least squares meshfree (DLSM) which uses an irreducible form of the governing equations, the MDLSM uses a mixed form of the original governing equations allowing for direct calculation of the gradients leading to more accurate computational results. In this study, an Eulerian–Lagrangian MDLSM method is proposed to solve incompressible multiphase flow problems. In the Eulerian step, the MDLSM method is used to solve the governing phase averaged Navier–Stokes equations discretized at fixed nodal points to get the velocity and pressure fields. A Lagrangian based approach is then used to track different flow phases indexed by a set of marker points. The velocities of marker points are calculated by interpolating the velocity of fixed nodal points using a kernel approximation, which are then used to move the marker points as Lagrangian particles to track phases. To avoid unphysical clustering and dispersing of the marker points, as a common drawback of Lagrangian point tracking methods, a new approach is proposed to smooth the distribution of marker points. The hybrid Eulerian and Lagrangian characteristics of the approach used here provides clear advantages for the proposed method. Since the nodal points are static on the Eulerian step, the time-consuming moving least squares (MLS) approximation is implemented only once making the proposed method more efficient than corresponding fully Lagrangian methods. Furthermore, phases can be simply tracked using the Lagrangian phase tracking procedure. Efficiency of the proposed MDLSM multiphase method is evaluated using several benchmark problems and the results are presented and discussed. The results verify the efficiency and accuracy of the proposed method for solving multiphase flow problems. © 2019 Elsevier Inc. All rights reserved.

Corresponding author. E-mail addresses: [email protected] (M.H. Afshar), [email protected] (S. Dabiri).

https://doi.org/10.1016/j.apm.2019.06.002 0307-904X/© 2019 Elsevier Inc. All rights reserved.

194

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

1. Introduction In recent decades, there has been a growing interest in meshfree methods. The meshfree methods have been developed and used to solve partial differential equations (PDEs) as an alternative to the mesh-based methods such as finite element method (FEM) and finite volume method (FVM). Meshfree methods use a set of scattered nodal points instead of elements/cells to discretize the computational domain. The mesh generation procedure, which can be costly and complicated for the complex geometries, is therefore replaced with a relatively simple node generation process. Unlike the mesh-based methods, which must use a predefined mesh to provide a connection between the nodes, no predefined information is needed for the scattered nodes [1]. Hence, the meshfree methods are more flexible than mesh-based methods in arranging the nodes. The higher flexibility of relocating the scattered nodes allows for adaptively changing the position of nodes to represent the changes in the geometry often encountered in the problems with dynamic boundaries or interfaces. Furthermore, it is reported that the accuracy of some meshfree methods is higher than that of the corresponding mesh-based methods [2]. Two major approaches are normally used in meshfree methods to approximate the primary unknown variables: kernel and series representation approximation. The smoothed particle hydrodynamic (SPH) and moving particle semi-implicit (MPS), are the well-known meshfree methods using the kernel approximation. The SPH was originally used to simulate astrophysical problems [3,4], but was later developed for solving the Newtonian fluid dynamics problems [5,6]. One of the commonly used forms of the SPH method, called weakly compressible smoothed particle hydrodynamic (WCSPH), includes compressibility effects through an equation of state to compute the pressure field [7–9]. Since WCSPH is an explicit method, it is simple to code and can efficiently handle fine nodal distributions. However, for fluid flow problems characterized by high Reynolds numbers, the method leads to large fluid density variations and unphysical void regions in the domain [10,11]. It is reported that the reliability of the WCSPH is much lower than the truly incompressible SPH (ISPH) using the projection method [12–14]. The MPS method is based on the Lagrangian movements of nodes [15,16]. In the classic form of the MPS, the spatial derivatives are approximated by a local weighted averaging method. Even though the MPS and SPH have some similarities, the tensile instability in the MPS can lead to more computational issues compared to the corresponding SPH [17] mainly due to the fact that the classic form of the MPS cannot guarantee the momentum and energy conservation. To overcome the mentioned drawbacks, a modified version of MPS is proposed [18]. The stability is also enhanced using an accurate approximation for Laplacian operator [19]. The error-compensating term is also considered in the source term of the Poisson pressure equation to stabilize the method [17,20–22]. The weakly compressible form of the MPS has been also presented [23] to decrease the computational cost compared to fully incompressible form utilizing the projection method. Although the kernel approximation is more efficient than the series representation method, some difficulties may arise when applying the standard form of the kernel approximations for the meshfree methods due to the fact that the standard form of the kernel approximation cannot satisfy the second order consistency [1,24]. Since the second order consistency is essential for solving quadratic PDEs by the strong form of meshfree methods such as SPH and MPS, more recently some techniques have been suggested to retain the second order consistency for the kernel approximation [25–27] at on increased computational cost. The Element Free Galerkin (EFG) and Meshfree Local Petrov Galerkin (MLPG) are the well-known meshfree methods based on the series representation approximation. The EFG has been originally suggested by Belytschko et al. for modeling solid mechanics problems [28–30] and was later extended to solve problems in fluid dynamics [31–33] and heat transfer [34]. Unlike the SPH and MPS, the EFG is a weak form method. In the weak forms, the order of governing PDEs is reduced by one order leading to the requirement of lower order consistency of the trial functions. These methods, however, require a background mesh for numerical integration. Hence, most researchers do not consider EFG as a truly meshfree method. MLPG, which uses a numerical integration on a local nodal sub-domain, has been proposed to avoid the requirement of background mesh in the weak form of meshfree methods. The method was first applied to the linear Poisson’s equations by Atluri and Zhu [35] and later was extended for solving the fluid mechanic problems [36,37]. For solving incompressible Navier–Stokes equations, the method has been modified to overcome the so-called Babuska–Brezzi stability conditions [38]. The upwind scheme is also applied in the MLPG to solve the convection-dominated problems [38,39]. The stability of the method is enhanced for convection–diffusion problems with large Peclet numbers [40]. Recently, the Discrete Least Squares Meshfree (DLSM) method is suggested using the strong form of the governing differential equations [41]. The method is considered as a truly meshfree method unlike the meshfree methods based on the weak form such as EFG. The method minimizes the residual function defined as the summation of the residuals of the governing PDEs and its boundary conditions. Since the least squares approach is used, the coefficient matrix is always symmetric and positive-definite. The method can therefore be overcome to the Ladyzenskaja–Babuska–Brezzi (LBB) condition. The method uses the series representation approximation and, hence, the consistency and accuracy provided by the method is higher than the kernel based meshfree methods such as SPH and MPS. The method is developed for solving the equilibrium [42–44] and propagation problems [45–48]. Recently, the mixed formulation of the method called Mixed DLSM (MDLSM) is proposed [49–52]. The Second order governing PDE is replaced by first order equations in the MDLSM and, therefore, costly procedure of computing the second derivatives of shape functions is avoided in the MDLSM method at the expense of increased size of the coefficient matrix. Furthermore, since the unknown variables and their gradients are simultaneously computed with the same order of accuracy, the gradients are calculated more accurately in the MDLSM method compared to the DLSM which uses a post-processing to obtain the gradients [49–51]. The MDLSM has been efficiently ap-

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

195

plied for simulating solid mechanics [49,53] and advection-diffusion problems [51]. More recently, the MDLSM method is developed and successfully used to solve the single phase Navier–stokes equations [54]. While the literature is rather rich on the application of DLSM or MDLSM to fluid mechanics problems, none of these methods has yet been used to solve the multiphase flow problems. The present study, therefore, is aimed at developing the MDLSM method to solve the multiphase flow problems. Two main approaches namely Eulerian and Lagrangian can be used for multiphase simulations and in particular capturing phase interfaces. Since in the Lagrangian computation, meshes/nodes are attached to the phases and moved in the path of the phases, different phases and their interfaces can be easily determined. While some mesh-based methods are successfully used the Lagrangian computational approach to solve the multiphase flow problems [55], Lagrangian simulation can lead to the computational difficulties due to mesh distortion. Hence, mesh-based methods mainly use Eulerian approach coupled with an interface-capturing algorithm such as volume of fluid (VOF) [56] or level set (LS) [57] method. Implementation of Lagrangian approaches is more convenient in the meshfree methods compared to the mesh-based ones. Hence, most meshfree methods often use the Lagrangian approach to solve the multiphase problems [58]. Tiwari and Kuhnert proposed a phase averaged finite point method (FPM) to solve the two-phase flow problems [59]. Hu and Adams presented a multiphase method based on WCSPH [60] and ISPH [61]. The SPH was also used to simulate different type of multiphase problems such as particle-fluid fluidization [62], dynamic and contact angle of droplet [63], complex interfaces and high density ratio [64,65], sediment scouring and flushing [66], and droplet coalescence and electro-coalescence [67]. Shakibaeinia and Chung developed a weakly compressible MPS model for multiphase flows [68]. Khayyer and Gotoh used MPS for modeling the multiphase problems with high density ratio [21]. While fully Lagrangian computations have been successfully used to solve multiphase flow problems, some methods combined the Eulerian and Lagrangian computations in order to benefit from the advantages of both. Recently, an Eulerian–Lagrangian incompressible SPH was presented using the Eulerian computation in regions away from the free surface and utilizing the Lagrangian approach near the free surface [69,70]. A hybrid Eulerian FVM and Lagrangian MPS method was also suggested to solve the two-phase flow problems [71]. In the current study, an Eulerian–Lagrangian MDLSM method is proposed to solve incompressible multiphase flow problems. In the proposed method, two types of independent points namely nodal and marker points are used. In the Eulerian part, the governing PDEs are solved using the MDLSM method on fixed nodal points to compute the velocity and pressure fields. The velocities of marker points, which are indexed the phases, are then interpolated with respect to the velocities of the nodal points using the affordable kernel approximation, in place of the costly MLS. Using the interpolated velocities, the marker points move as the Lagrangian particles to track the different phases. A special method is also used for smoothing the distribution of marker points to avoid unphysical clustering and dispersing of points, which is a drawback of the Lagrangian meshfree methods. The nodal points are fixed on the Eulerian part, thus, the costly MLS approximation and its derivatives are calculated only once, which remarkably reduces the computational cost compared to the corresponding fully Lagrangian method. Furthermore, the different phases are easily tracked and the properties of the phases can be simply captured using the information carried by the marker points. 2. Moving least squares (MLS) approximation In the current study, the MLS approximation is employed for constructing the shape functions [72]. The method has also been applied in the EFG and MLPG for approximating the unknown variables [28,35,73]. In the MLS, the unknown function, φ , is approximated by a series polynomial function as follows:

φ (X ) = P T (X )a(X );

(1)

where a(X) contains the coefficient vector, which is defined as a function of X denoting the coordinate of the nodal points. P is a complete polynomial basis function of order m, which must be chosen with respect to the required order of consistency and accuracy. For example, P is defined for a two-dimensional problem as:



P T (X ) = 1, x1 , x2 , x1 2 , x1 x2 , x2 2 , . . . , x1 m , x1 m−1 x2 , . . . , x1 x2 m−1 , x2 m



1×k

; X = ( x1 , x2 )

(2)

A weighted least squares functional, J, is defined as:

J=

ns 

  

w j (X − X j ) P T X j a(X ) − φˆ j

2

;

(3)

j=1

where ns is the number of nodes located inside the support domain and φˆ j describes the unknown nodal parameter at jth node. The following cubic spline function is applied as a weight function at jth node:

w j (d ) =

⎧2 2 3 ⎨ 3 − 4d + 4d 4

⎩3

d≤

− 4d + 4d2 − 43 d3

1 2

1 2

≤ d ≤ 1 ; d = X − X j /dw j

(4)

d≥1

0

Here, dwj is the radius of the support domain for jth node located at Xj . Minimizing J with respect to a(X) yields the following relation:

a(X ) = A−1 (X )B(X )φˆ ; φˆ T =



φˆ 1

φˆ 2

...

φˆ j

...

φˆ ns



(5)

196

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Based on the relations described by Eq. (5) and Eq. (1), the unknown nodal value φ is approximated in terms of the unknown nodal parameter φˆ as:

φ (X ) = N (X )φˆ ;

(6)

where the shape function N(X) is

N (X ) = P T (X )A−1 (X )B(X );

(7)

A(X), the moment matrix, and B(X) are

A (X ) =

ns 



  

 

w j X − X j P X j PT X j ,

(8)

j=1





 





 



B(X ) = w1 X − X 1 P X 1 , w2 X − X 2 P X 2 , . . . , wns (X − X ns )P (X ns ) .

(9)

The first order derivative in terms of ith component of the coordination (xi ) can be computed as follows:

∂ N ∂ P T −1 ∂ A−1 ∂B = A B + PT B + P T A−1 . ∂ xi ∂ xi ∂ xi ∂ xi

(10)

The number of nodes located inside the support domain must be at least equal to the number of terms in the basis function (k) to satisfy the essential condition for the invertibility of moment matrix A(X). To satisfy this condition, the radius of the support domain at node j is defined by dwj = λdjk and λ is larger than unity, where djk is the distance of the node j from the kth nearest nodal point. The usual form of the MLS approximation employed in the current study, does not satisfy the Kronecker delta property. Therefore, direct imposition of the Dirichlete boundary conditions is not possible and a special approach must be used (see the discussion in Section 3.3.2). 3. Governing equations and numerical scheme 3.1. Governing equations for multiphase flow problems Based on the one-fluid model, the mass and momentum conservation equations for immiscible and incompressible multiphase system are described as:

∂ ui = 0, on  ∂ xi

(11)

1 ∂P ∂ ui ∂ ui 1 ∂ τi j 1 + uj =− + + gi + Fs i ; on  ∂t ∂xj ρ˜ ∂ xi ρ˜ ∂ x j ρ˜

(12)

where, ui is the averaged Cartesian components of the velocity field, P is the averaged pressure, gi is the source term, ρ˜ is the phase averaged density, and Fs is the surface tension force. For Newtonian incompressible fluids, the gradient of the averaged stress tensor τ ij can be expressed by:







∂ τi j ∂ ∂ μ˜ ∂ ui ∂ μ˜ ∂ u j ∂ ui ∂ u j ∂ ∂ ui ∂ ∂uj = μ˜ + = + +μ ˜ +μ ˜ ; ∂xj ∂xj ∂ x j ∂ xi ∂ x j ∂ x j ∂ x j ∂ xi ∂xj ∂xj ∂ xi ∂ x j

(13)

where μ ˜ is the phase averaged viscosity. Eq. (13) can be simplified using the continuity equation, Eq. (11), as:



∂ τi j ∂ μ˜ ∂ ui ∂ μ˜ ∂ u j ∂ ∂ ui = + +μ ˜ . ∂xj ∂ x j ∂ x j ∂ x j ∂ xi ∂xj ∂xj

(14)

Finally, the momentum equation, Eq. (12), can be rewritten as [59]:

1 ∂P ∂ ui ∂ ui 1 ∂μ ˜ ∂ ui ∂ μ˜ ∂ u j 1 ∂ 2 ui + uj =− + + +μ ˜ + gi + Fs i . ∂t ∂xj ρ˜ ∂ xi ρ˜ ∂ x j ∂ x j ∂ x j ∂ xi ρ˜ ∂ x2j

(15)

The Dirichlet boundary conditions for velocity are defined as:

ui ηi = u¯ η ui ti = u¯ t ;

on

(16)

where ηi and ti describe the unit normal and tangent vectors to the boundary, respectively. u¯ η and u¯ t are the normal and tangential components of the known velocity vector to boundaries, respectively.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

197

3.2. Computing the surface tension force The following relation is used to calculate the surface tension force.

Fs i = Fs = σ κ ns δs ;

(17)

where ns is the unit vector normal to the surface, σ is the surface tension coefficient, κ denotes the curvature of the interface, and δ s is the surface delta function. The color function, c, is often used to calculate the interface properties, which is taken c = 1 on the target phase and c = 0 on the other phases. The color function is smoothed by:

n s

c˜ = k=1 ns

wk c k

k=1

;

wk

(18)

where ck and wk are the color and kernel functions at the kth node, respectively. The cubic spline function, described by Eq. (4), is used as a kernel function. The unit normal vector and curvature of surface can be computed by:

ns =

∇ c˜ , |∇ c˜|

(19)

κ = −∇ .ns ;

(20)

while δ s is calculated by using the following relation:

δs = |∇ c˜|.

(21)

First derivatives of MLS shape functions are used to compute the gradient operator (∇ ). 3.3. Proposed MDLSM method for solving the governing PDEs of multiphase problem In this section, the solution procedure for the governing PDEs, represented by Eqs. (11) and (15), is described. The solution procedure consists of a discretization in time followed by a spatial discretization in space using the MDLSM method. The shape functions and their derivatives are required to be calculated once since a fixed nodal distribution is used. Hence, time-consuming MLS approximation is avoided in each time step unlike fully Lagrangian approaches. In what follows, the projection method, used for time discretization, is described which is followed by the MDLSM for spatial discretization. 3.3.1. Projection method In the projection method, each time step is divided into two pseudo-time steps. The split form of the momentum equation in time is then written as:

uni +1 − u∗i + u∗i − uni

t

+ uj

∂ ui 1 ∂P 1 ∂ τi j 1 =− + + gi + Fs i ; ∂xj ρ˜ ∂ xi ρ˜ ∂ x j ρ˜

(22)

where the superscript containing n denotes the corresponding time step. In the semi-incremental fractional step prediction/correction approach, the intermediate velocity is explicitly predicted in the first step by neglecting the implicit part of the pressure leading to the following decoupled equation for the velocity components as:



u∗i

=

uni

+ t

−unj

∂ 2 uni ∂ uni 1 ∂μ ˜ ∂ uni ∂ μ˜ ∂ unj 1 γ ∂ Pn + + +μ ˜ + gi + Fs i − ; ∂ x j ρ˜ ∂ x j ∂ x j ∂ x j ∂ xi ρ˜ ρ˜ ∂ xi ∂ x2j

(23)

While an explicit approach to calculate the velocity reduces the computational cost significantly, it imposes very strict stability conditions. A discussion on the stability of the method is, therefore, presented in Appendix A for a typical simplified non-linear equation. The effect of the implicit part of the pressure is considered in the second step as:



uni +1 = u∗i + t −

1 ρ˜

∂ Pn+1 γ ∂ Pn + . ∂ xi ρ˜ ∂ xi

(24)

Applying the divergence operator to Eq. (24) and considering the continuity equation, the pressure Poisson equation is derived as:





∂ 1 ∂ Pn+1 1 ∂ u∗i ∂ 1 ∂ Pn = +γ ; ∂ xi ρ˜ ∂ xi

t ∂ xi ∂ xi ρ˜ ∂ xi

(25)

where γ , named fractional parameter, is a user-specified constant taken between 0 and 1. The method is equivalent to the well-known non-incremental and incremental fractional step methods, assuming γ = 0 and γ = 1, respectively [46]. An

198

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

appropriate value for γ can provide the possibility of using larger time step sizes compared to the corresponding nonincremental one. A more detailed discussion on the appropriate value of γ is available in the literature [46]. The projection method is normally implemented in three steps. The intermediate velocity is calculated using Eq. (23) at the first step. The pressure field is then computed in the second step by solving the pressure Poisson equation, represented by Eq. (25), and finally the velocity field at the next time step is obtained using Eq. (24). The Dirichlet boundary conditions for the velocity, described by Eq. (16), can be directly imposed in Eqs. (23) and (24). However, the effect of such boundary conditions must be also considered in the pressure Poisson equation. Based on Eq. (12), the momentum equations along the normal and tangent to the boundary can be written as:

∂ P 1 ∂ τi j 1 + + gi + Fs i ηi , ∂ xi ρ˜ ∂ x j ρ˜

1 ∂P Du i Du¯ t 1 ∂ τi j 1 ti = = − + + gi + Fs i ti ; Dt Dt ρ˜ ∂ xi ρ˜ ∂ x j ρ˜

Du¯ Du i ηi = η = Dt Dt



1 − ρ˜

(26)

where D(.)i /Dt = ∂ (.)i /∂ t + uj (∂ (.)i /∂ xj ). The condition of slip boundary conditions at static walls where u¯ η = 0 and Du¯ η /Dt = 0, can be represented by:

∂P η = ∂ xi i



∂ τi j + ρ˜ gi + Fs i ηi . ∂xj

(27)

In a similar way, the condition of no-slip boundary conditions at static walls, where u¯ η = 0, Du¯ η /Dt = 0 and u¯ t = 0, Du¯ t /Dt = 0, can be defined by:

∂P t = ∂ xi i



∂ τi j + ρ˜ gi + Fs i ti . ∂xj

(28)

The pressure boundary condition is simply imposed as:

P = P¯ ;

(29)

where P¯ is the specified pressure at the boundary. 3.3.2. Proposed MDLSM method for pressure Poisson equation The pressure Poisson equation, described by Eq. (25), is discretized by the proposed MDLSM method. In the proposed mixed formulation, the gradients of the pressure are considered as new unknown variables defined as:

∂P = qi . ∂ xi

(30)

Substituting Eq. (30) into the pressure Poisson equation, Eq. (25), yields the following system of differential equations:

⎧ n+1 ∂P ⎪ ⎨ − qni +1 = 0 ∂ xi ;   ⎪ ⎩ ∂∂x 1 qni +1 = Sn i ρ˜

with

1 S =

t n

(31)



∂ u∗i ∂ 1 ∂ Pn +γ ; ∂ xi ∂ xi ρ˜ ∂ xi

(32)

This system of equations can be written in a compact form using the matrix notation as:

Ei

∂φ + Gφ = S; ∂ xi

(33)

where the unknown nodal value φ and source term S are defined in 2-D problems as:

φ=



P



S= 0 with

⎡

1 E 1 = ⎣0 0

q1 0

q2

n+1

,

(34)



Sn ;

0 0

∂ (1/ρ˜ ) ∂ x1

(35)

⎤

0 0⎦ 0

3 × 3

⎡

0 , E2 = ⎣1 0

0 0 0

0 0

∂ (1/ρ˜ ) ∂ x2

⎤



⎦

G= 3 × 3

0 0 0

−1 0 0

0 −1 0

 ; 3 × 3

(36)

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

199

Fig. 1. The normal and tangent vectors on the boundary of an arbitrary two-dimensional problem.

Similarly, the Dirichlet and Neumann boundary conditions of Eqs. (27)–(29) can be re-written in a matrix form as a Dirichlet type conditions defined as:

with

Dφ = φ¯

(37)

  φ¯ T = P¯, Q¯ η , Q¯t ,

(38)

and

 D=

1 0 0

0

0



η1

η2 ;

t1

t2

(39)

where η = [η1 ,η2 ] and t = [t1 ,t2 ] denote the normal and tangent vectors to the boundary, respectively.while Q¯ η and Q¯t are





∂ τ1 j ∂ τ2 j + ρ˜ g1 + Fs 1 η1 + + ρ˜ g2 + Fs 2 η2 , ∂xj ∂xj



∂ τ1 j ∂ τ2 j ¯ + ρ˜ g1 + Fs 1 t1 + + ρ˜ g2 + Fs 2 t2 ; Qt = ∂xj ∂xj

Q¯ η =

(40)

The first row of system of Eq. (37) corresponds to the imposition of the specified pressure condition P¯ ; the second equation corresponds to the enforcement of pressure gradient normal to the boundary, and finally the third equation represents the specification of the tangential pressure gradients. Fig. 1 shows the tangent and normal vectors on the boundary of an arbitrary two-dimensional problem. Discretizing the computational domain with a set of nt nodal points arbitrarily chosen on the problem domain and its boundaries, the MLS approximation is used for function approximation. Using MLS shape function to independently approximate each component of the unknown vector φ, leads to the following MLS approximation:

φ (X ) = Nmatrix (X )φˆ ;

(41)

where Nmatrix is the shape function matrix defined as:





Nmatrix (X ) = N 1 (X ), N 2 (X ), . . . , N k , . . . , N nt (X ) ; with

 N (X ) = k

N k (X ) 0 0

0 N (X ) 0 k

(42)



0 ; 0 N k (X )

(43)

where Nk (X) denotes the shape function of kth node at X computed using Eq. (7). The nodal parameter vector, φˆ , is stated as:

with

 n+1 φˆ = φˆ 1 , φˆ 2 , . . . , φˆ k , . . . , φˆ nt ;

(44)

  φˆ k = Pk , qk1 , qk2 .

(45)

200

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Pk denotes the pressure at kth node, qk1 and qk2 are the corresponding pressure gradients in x1 and x2 directions at point k, respectively. The residual of governing PDE of Eq. (33) and the residual of corresponding boundary conditions of Eq. (37) at a typical point X are now defined as:

R ( X ) = L



R ( X ) = D

 φ (X ) − S (X ) = L(Nmatrix (X ) )φˆ − S (X );

(46)



 φ (X ) − φ¯ (X ) = D(Nmatrix (X ) )φˆ − φ¯ (X );

(47)

where the differential operator, L, is defined as:

L ( . ) = Ei

∂ (. ) + G (. ). ∂ xi

(48)

The least-squares residual functional is now defined as:



1 R= 2

nt  

 k 2

R X

+α

k=1

nb 

  k 2



;

R X

(49)

k=1

where Xk is the coordinate at node k, and nb represents the number of nodes on the boundaries. Since the MLS shape functions used here does not satisfy the Kronecker delta property, the boundary conditions cannot be directly imposed. A penalty method is, therefore, used here to enforce the boundary conditions which are all of Dirichlet type. The penalty parameter, α , should be large enough to impose the desired boundary conditions [1,24]. Larger value of the penalty parameter leads to exact satisfaction of the boundary conditions. This will, however, increase the condition number of the coefficient matrix leading to the ill-posed system of equations which may in turn reduce the accuracy of the solution procedure [74]. A brief discussion on the effect of the penalty parameter on the condition number of coefficient matrix is presented in Appendix B while more detailed discussion is available in the literature [51,75]. Minimizing Eq. (49) with respect to the unknown nodal parameters φˆ leads to: nt 

   ∂ R X k  ∂ φˆ

k=1

 k 

R X

 

nb    ∂ R X k +α R X k = 0. ∂ φˆ

(50)

k=1

Substituting Eqs. (41), (46) and (47) into Eq. (50) yields the following symmetric positive-definite system of equations:

K φˆ = F ;

(51)

where the coefficient matrix, K, and right-hand-side vector, F, are calculated as:

K=

nt  



 T  

L Nmatrix X k

 

L Nmatrix X k

+α

k=1

F=

nt  

nb  



D Nmatrix X k

T  

 

D Nmatrix X k

;

(52)

k=1



 T  k 

L Nmatrix X k

k=1

S X

+α

nb  



D Nmatrix (X k

T  k  φ¯ X .

(53)

k=1

The unknown nodal parameters φˆ can now be obtained by solving the algebraic system of Eq. (51). Since the Kronecker delta property is not satisfied using the MLS approximation, the nodal parameters are not equal to the nodal values. The nodal values must, therefore, be retrieved by applying Eq. (41) using the obtained nodal parameters. It is seen that application of MDLSM method to the second order pressure Poisson equation converts to the first order PDE only requires first order derivative of the shape function, which is both more accurate and efficient compared to second order derivatives required in the DLSM method. Moreover, the pressure and its gradients are simultaneously and independently approximated using the same MLS shape function, which naturally leads to more accurate results for the pressure gradients to be used in the momentum equations compared to the DLSM method, which uses a post-processing approach to calculate the pressure gradients [51,54]. While the proposed MDLSM method avoids the costly calculation of the second order MLS shape function derivatives, like any other mixed approach increases the number of unknowns in the pressure Poisson equation from one to three in 2-D and to four in 3-D problems. 3.4. Particle tracking scheme As noted previously, two types of independently distributed points are used in the proposed Eulerian–Lagrangian method to simulate the underlying multiphase problem namely nodal points and marker points. The nodal points are static and used in the Eulerian part of the method while the dynamic marker points are used in the Lagrangian part. In the Eulerian part, the flow properties, including the velocity and pressure fields, are calculated on the static nodal points, which are used to discretize the computational domain. A set of marker points, acting as Lagrangian particles, are then employed for tracking the phases. The proposed particle tracking procedure contains two steps: marker points tracking and marker points smoothing.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

201

3.4.1. Marker point tracking Once the velocity and pressure fields are computed at the nodal points, the velocities at the marker points are computed using the following interpolation of the velocity field at the nodal points:

n s

U _m =

wkU k

k=1  ns

k=1

wk

, U k = uki ;

(54)

where wk is the value of cubic spline function, defined in Eq. (4), at kth nodal point. U _m is the velocity vector at the marker points, and Uk denotes the velocity vector at kth nodal point. ns is the number of nodal points located inside the support domain of the considered marker point. The velocity of the marker points is then used to calculate the location of marker points in the next time step as follows:



X _mn+1 =

U _mn+1 + U _mn

t + X _mn ; 2

(55)

here X _mn and U _mn are the coordinate and velocity vector of marker points at nth time step, respectively. Since the marker points acts as Lagrangian particles, the density and viscosity are constant on each particle path over time for incompressible fluids. Hence, the following relations hold:

dρ _m(X ) = 0, dt

(56)

dμ_m(X ) = 0. dt

(57)

The averaged density, ρ˜ , and viscosity, μ ˜ , at the nodal points, required for the solution of the governing equations (Eqs. (11) and (12)), are now calculated at the nodal points in terms of the density and viscosity of neighboring marker points as follows:

ρ˜ = μ˜ =

n s _ m

ρ

wk _mk k=1 ns_m k , w k=1

n s _ m

wk k=1 ns _m k=1

μ_mk wk

;

(58)

(59)

where, ρ _mk and μ_mk are the density and viscosity at kth marker point, respectively and ns_m denotes the number of marker points in the support domain of the considered node. The nodal and marker points and their support domains are schematically illustrated in Fig. 2. A discussion is presented in Appendix C to investigate the effect of radius of support domain on stability and accuracy. 3.4.2. Special consideration around the boundaries The velocity normal to the wall boundaries is equal to zero so the marker points located on the walls will remain on the wall. However, Eq. (54) may over-estimate the velocity of the marker points next to the walls leading to large displacement calculated by Eq. (55) causing these points to cross the wall boundary which obviously undermines the performance of the particle tracking procedure. Such problem especially may arise when the large time step sizes are used. Special considerations are, therefore, considered to circumvent this problem and keep the marker points inside the computational domain. The special consideration is schematically illustrated in Fig. 3. For the marker points located next to the walls, the velocity component normal to the wall is virtually set to zero, U _mn+1 .η = U _mn .η = 0, in case the calculated coordinate using Eq. (55), X _mn+1 , makes the points to cross the wall boundaries. 3.4.3. Marker point smoothing Since the marker points move with an interpolated velocity obtained from the velocity of nodal points, the marker points may cluster in some regions or disperse in some other leading to the void regions in the computational domain, which could in turn deteriorate the accuracy of the approximation in the crowded or void regions. The same drawback has been also reported and discussed for the methods such as SPH and MPS which are based on the Lagrangian approach [71,76–78]. This problem can be resolved or at least moderated by using a finer discretization at the expense of increased computational cost. In the current study, the drawback is circumvented by using a smoothing procedure to find the marker point distribution. The proposed smoothing procedure is inspired by the refinement strategies used in solid [42,53] and fluid [79,80] mechanics problems. The purpose of the smoothing procedure is avoid clustering and dispersing of marker points. In this method, the marker points are connected to each other using a Delaunay diagram. The set of points that are connected to the point under consideration are defined as the neighboring points. Fig. 4 illustrates the connection between the points in a 2-D problem along with the Delaunay diagram, and the neighboring points. Since the smoothing procedure is not based on the physical nature of the flow problem, the smoothing procedure must not change the position of points located on the interface and boundaries. Hence, the marker points are fixed on both the

202

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 2. Nodal and marker points and their support domains.

Fig. 3. Schematically illustrating the special considerations to avoid penetration of points to the walls.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

203

Fig. 4. Connected marker points using the Delaunay diagram, and neighborhood and interface points.

boundaries and the interfaces. A point is considered as an interface point if its neighboring points contain at least one marker point of a different phase. This can be mathematically defined as:

k ∈ Cinter f ace i f C k neighborhood ∩





i= p

C _mi = ∅,

∀ k ∈ C _m p ;

(60)

where C _m p denotes the set of marker points of the pth phase, Ck neighborhood represents the set of neighboring marker points of the point k, and Cinterface denotes the set of points located on the interface. Fig. 4 clearly shows that only the interface points are connected to points of different phases. Once the interface points are recognized, each marker point is connected to its neighboring points by a set of virtual springs. The spring system is established in each phase separately which is schematically shown in Fig. 5 for a 2D problem. For a spring connecting the marker points 1 and 2, the following system of equations, representing the equilibrium of the end points, can be written as:



1 s12 F _smoothing 2 s12 F _smoothing





=

s12 Z

−s12 Z

−s12 Z

s12 Z



X  _m1 X  _m2



;

(61)

where X  _m1 and X  _m2 are the modified coordinate of marker points in smoothed distribution, respectively. 1 2 s12 F _smoothing and s12 F _smoothing are respectively the local point forces at points 1 and 2 produced by the action of the spring s12, and s12 Z is the local stiffness matrix. Fig. 5 shows the point force for the end point 1. Since a regular distribution of the marker points is desired, the local stiffness for all springs set equal to unity defined as: s (. ) Z

= I;

(62)

where I is the identity matrix. Using the assembly process, the global system of equations for the entire system of virtual springs can be written as:

K _smoothing X  _m = F _smoothing;

(63)

where K _smooting and F _smoothing represent the coefficient matrix and the force vector for the entire spring system, respectively, and X  _m is the smoothed coordinates of the marker points. For equilibrium, the global force vector, F _smoothing,

204

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 5. Virtual spring system for each phase.

is set to zero leading to the following linear system of equations to be solved for the new coordinates of the marker point.

K _smoothing X  _m = 0.

(64)

The system of Eq. (64) is singular before imposing any boundary conditions. The following essential boundary conditions are used to fix the boundary and interface points:

X  _mk = X _mk ;

∀ k ∈ Cinter f ace ∨ k ∈ Cboundary ;

(65)

where Cinterface and Cboundary are the set of marker points located on the interfaces and boundaries, respectively. Since the computational cost for the smoothing procedure is relatively high, the process is only applied when it is needed. The flowchart of the proposed MDLSM method is described in Fig. 6. 4. Numerical examples In this section, the efficiency and accuracy of the proposed multiphase model is tested on several benchmark problems. Unless otherwise stated, the linear basis function is used for constructing the MLS shape functions. The penalty parameter, α , is taken as 108 and γ is set as 0.5. The smoothing procedure for the marker point distribution is applied in every 10 time steps. 4.1. Laplace’s law A circular drop surrounded by a background fluid of different phase is modeled here as a well-known problem for evaluating the surface tension. The computational domain is a unit square. No-slip boundary condition is applied at all boundaries, and the pressure is set to zero at the walls. The pressure jump across the interface can be obtained by the Laplace’s law defined as:

Pin − Pout = σ κ =

σ

R

;

(66)

where Pin and Pout denote the inner and outer pressure values with Pout taken as zero, and R is the radius of the drop. The Density of background fluid is set to 2 (kg/m3 ) while the density of drop is equal to 1 (kg/m3 ). The problem is defined in a horizontal plane so the gravity is set to zero in both directions. Other variables are taken as R = 0.2 (m), σ = 1 (N/m), μ = 10−2 (Pa .s). The problem is solved using different nodal distributions shown in Fig. 7. The pressure field, calculated

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

205

Fig. 6. The flowchart of MDLSM method for multiphase problems.

Fig. 7. Nodal and marker points distributions for the Laplace’s law problem.

by the proposed method at 200th time step, is shown in Fig. 8. In Fig. 9, the averaged pressure inside the drop is presented and compared with the analytical pressure jump calculated using the Laplace’s law. The results clearly show the ability of the proposed method to accurately solve this problem. Furthermore, it is seen that accuracy of results is improved by using the finer nodal distributions as expected. The convergence rate of the method for the problem is shown in Fig. 10. The convergence property of the MDLSM method is also investigated with more details in the literature [49,52] and Appendix D.

206

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 8. Pressure field for the Laplace’s law problem.

4.2. Two-fluid sloshing This example considers the sloshing of two fluids of different density. The heavy fluid is located above the lighter one. The density of heavy and light fluids are ρ = 2 (kg/m3 ) and ρ = 1 (kg/m3 ), respectively. The dynamic viscosity is set to μ = 10−3 (Pa .s) for both fluids and gravity is taken as g = −0.294 (m/s2 ) in the vertical direction. Details of initial conditions for the problem are illustrated in Fig. 11. The no-slip and slip boundary conditions are applied on horizontal and vertical walls, respectively. An oscillatory behavior is expected on the interface until the systems comes to an equilibrium state in which the interface is parallel to the horizontal plane. Three uniform nodal point distributions of 9 × 7, 17 × 13 and 41 × 31 nodes are used to solve the problem using three uniform marker point distributions of 17 × 13, 33 × 25 and 81 × 61, respectively. The relative height of the interface at the reference point “A” is shown and compared with benchmark results [81,82] in Fig. 12. Comparing the results of proposed MDLSM method with the benchmark results verifies the acceptable efficiency of the proposed method. The position of the interface at t = 50 (s), obtained on 41 × 31 nodal distribution, is shown in Fig. 13. The horizontal interface computed by the proposed MDLSM method is seen to be fairly in agreement with those of other numerical predictions available in the literature [81,82], in spite of a tiny gap present between the phases on the right side near the interface. This issue will be briefly addressed later at the end of Section 4.4. 4.3. Rayleigh–Taylor instability The Rayleigh–Taylor instability test problem is widely used for evaluating the efficiency of multi-phase numerical methods [13,59,64,68]. The initial state is illustrated in Fig. 14 where the heavier fluid is located on top of the lighter fluid with a small perturbation in the interface defined as y = d + 0.1cos (2π x). The problem domain is rectangular with the dimension 1d × 2d. The width of domain is taken as d = 1 (m). Density of heavy and light fluids are set as ρ = 3 (kg/m3 ) and ρ = 1 (kg/m3 ), respectively. Kinematic viscosity for both fluids is taken ν = 0.0124 (m2 /s), and the gravity is set as g = −10 (m/s2 ) in the vertical direction. The slip boundary condition is used at the walls. The surface tension force is ignored in this test case. The time step size is t = 0.001 (s). Quadratic basis function is used to construct the MLS shape functions. The problem is solved on the nodal and marker point distributions illustrated in Fig. 14. Fig. 14 also shows the evolution of the interface at different times. Vertical position of middle bottom point of the heavier phase, marked as point B in  Fig. 14, and left/right top point of the lighter phase, point A, are shown in Fig. 15 in terms of the dimensionless time (t ∗ = t/ d/g) and compared with the existing results in the literature [83] using 62,500 particles in total. Points A and B represent the maximum and minimum vertical coordinates in the light and heavy phases, respectively. The problem is also solved imposing an upward disturbance at the interface defined as y = d − ∈ 0 cos (2π x/ξ ). The amplitude of the initial perturbation is taken as ∈ 0 = 0.05 with the wavelength chosen as ξ = d = 1 and the initial condition shown in Fig. 16. The other constant parameters are taken as ν = 0.01 (m2 /s), g = −10 (m/s2 ), ρ heavy phase = 3 (kg/m3 ) and ρ light phase = 1 (kg/m3 ). The variation of the computed dimensionless maximum interfacial displacement, ∈ ∗ = ∈/ξ , versus



dimensionless time (t ∗ = t/ ξ /g) is shown in Fig. 17 and compared to the VOF and weakly compressible MPS results [68]. The results obtained using the proposed multiphase MDLSM method are seen to be in well agreement with that predicted

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

207

Fig. 9. Averaged pressure inside the drop for the Laplace’s law problem.

using VOF and weakly compressible MPS considering the fact that the proposed MDLSM method uses coarser nodal distributions than VOF and weakly compressible MPS methods. 4.4. Heterogeneous flow with three different phases This example is used to evaluate the ability of the proposed method to solve density stratified fluids. Three different phases (A, B and C) are initially located as shown in Fig. 18. The computational domain is a 1 × 1 square. All fluids have the same volume. The same dynamic viscosity is taken for all fluids as μ = 10−1 (Pa. s). The density of the fluids are ρ A = 0.1 (kg/m3 ), ρ B = 1 (kg/m3 ) and ρ C = 2 (kg/m3 ). The initial 21 × 21 nodal point distribution is shown in Fig. 19. The entire computational domain is surrounded with the wall boundaries. The slip boundary condition is used at all walls. All of the pressure boundary conditions are of Neumann type. Hence, a reference pressure value, P = 0, is considered on the top-right corner of the domain to guarantee uniqueness of the solution. The position of phases obtained using the proposed method is illustrated in Fig. 19 at different times. The results at t = 0.4 (s) and t = 1 (s) are compared with the available results computed using particle finite element method with 41 × 41 nodes [55]. The phases almost reach the stable state at t = 5 (s) which is in agreement with what is expected theoretically and reported in the literature [55]. The obtained results indicate proper ability of the proposed MDLSM method to solve this problem. However, a tiny gap between the phases shown in blue and red colors is seen at t = 5 (s) at the left boundary. This issue is briefly discussed in Appendix C. While the gap may disappear for finer point distribution, it may well be due to the fact that the smoothing procedure is carried out separately. This problem can be solved using a special point enrichment process on the interface which will be investigated in future studies. 4.5. Dam break problem The dam break problem is a well-known test case to evaluate the multi-phase methods. The initial condition of the problem is shown in Fig. 20. Starting from the initial conditions, the water column collapses and flows to the downstream as reported by the experimental investigations [84]. This problem is considered as a more challenging test case for multi-phase methods due to the high density ratio between the phases. The width of the water column is considered as L = 0.6 (m). The gravity is set as g = −9.81 (m/s2 ). A slip boundary condition is imposed on the walls and surface tension is ignored. Three

208

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 10. Convergence rate of the proposed MDLSM method for the Laplace’s law problem.

Fig. 11. The initial condition for the two-fluid sloshing problem.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 12. The results of proposed MDLSM method for the sloshing problem obtained using different nodal distributions.

Fig. 13. Position of marker points at t = 50 for the sloshing problem using 41 × 31 nodal points.

209

210

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 14. Initial condition and snapshots of the Rayleigh–Taylor instability problem using the MDLSM method (downward disturbance).

different nodal and marker point distributions, shown in Fig. 21, are used to solve the problem. Snapshots of the phases and velocity vectors obtained on a 31 × 41 nodal point distribution are shown in Fig. 22. Fig. 23 compares  the dimensionless leading edge of water in horizontal direction (Xfront /L) as a function of the dimensionless time (t ∗ = t 2g/L) to the available experimental data [84,85]. For the leading edge position in time, the error norm of MDLSM method in terms of the VOF prediction [85] is presented in Table 1. The results of present MDLSM multiphase method, especially the one obtained on finer nodal distribution, are in good agreement with the reported experimental [84] and numerical [85,86] results.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

211

Fig. 15. The maximum and minimum vertical coordinates in the light and heavy phases for the Rayleigh-Taylor instability problem (downward disturbance).

Fig. 16. Initial condition and snapshots of the Rayleigh–Taylor instability problem using the MDLSM method (upward disturbance).

Fig. 17. Maximum interfacial displacement for the Rayleigh–Taylor instability problem (upward disturbance). Table 1 Error norm of MDLSM method, using different set of nodal distributions, with respect to the VOF method computed using 50 × 38 grid points (Nguyen et al. [85]). Nodal distribution

25 × 19

31 × 24

41 × 31

Error

0.080

0.047

0.013

212

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 18. Initial and final stable state for the heterogeneous flow problem.

Fig. 19. Initial condition and snapshots of heterogeneous flow problem based on the MDLSM method.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

213

Fig. 20. The initial condition for the dam break problem.

Fig. 21. Three different nodal and marker point distributions used for solving the dam break problem.

5. Discussion More recently, an Eulerian MDLSM method is presented to solve single-phase flow problems [54]. It has been shown that using the mixed formulation leads to more accurate solution for the MDLSM method compared to the existing DLSM method. In this paper, an Eulerian–Lagrangian MDLSM method is presented to solve multiphase flow problems with interfacial forces. A novel Lagrangian type particle tracking procedure is also presented that uses a particle smoothing procedure to enhance its performance. The proposed Eulerian–Lagrangian MDLSM multiphase method solves the governing PDEs on a set of fixed nodal point distribution and tracks the phases using a Lagrangian concept. Advantages of the proposed Eulerian–Lagrangian MDLSM method are multifold. First of all, the MLS approximation used for the shape function construction paves the way for the method to satisfy the required order of consistency, unlike the SPH and MPS methods using the well-known form of the kernel approximations. The proposed method, therefore, is not theoretically subjected to the sort of numerical problem reported and discussed in the literature [1,17,22,24,87–90] for the kernel-based meshfree methods. Drawbacks of the kernelbased methods mainly arise due to lack of required order of consistency and in particular essential second order consistency when dealing with the second order PDEs [1,24,88,90,91] which can deteriorate the quality of the final solution. This explains the fact that proposed MDLSM produces solutions with the same accuracy as those obtained by the SPH or MPS methods with coarser nodal distribution as illustrated in Figs. 15, 17 and 23. The superior accuracy and convergence rate of the

214

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. 22. Snapshots of dam break problem based on the MDLSM method, obtained using 41 × 31 nodal point distribution.

Fig. 23. The water front propagation along the bottom wall for the dam break problem.

proposed MDLSM method is also shown and compared with those obtained using the well-known form of the kernel-based meshfree methods in Appendix D for a quadratic PDE problem (Fig. D.1.). Superior performance of the proposed MDLSM method compared to the kernel-based meshfree methods comes at an extra cost. This is mostly due to the fact that the computational cost of the MLS approximation is remarkably higher than the corresponding kernel approximation. Table D.1 of Appendix D compares the computational cost of the MLS and kernel approximations. However, since the flow simulation part of the proposed MDLSM method is of Eulerian nature, time-

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

215

consuming MLS approximation is used only once while a fully Lagrangian version of the method requires the use of MLS approximation at each time step.

6. Conclusion Recent researches reported in the literature shows that the Mixed DLSM method is more accurate than the irreducible DLSM method mainly due to the fact that the mixed formulation used in the former method allows for direct calculation of the pressure gradient that is more accurate than those obtained from pressure using a post-processing procedure as used in the DLSM. The MDLSM method, as a truly meshfree method, was proposed in this paper to simulate multiphase flow problems. The method employs an Eulerian–Lagrangian approach. In the Eulerian part, the governing phase averaged Navier–Stokes equation were discretized using MDLSM method at a set of fixed nodal distribution and solved to obtain the velocity and pressure field while in the Lagrangian part, different phases were tracked. Using a fixed nodal distribution on Eulerian part allows us to implement the time-consuming MLS approximation only once. Hence, the computational cost is reduced compared to corresponding fully Lagrangian approaches using the MLS approximation at each time step. A Lagrangian phase tracking was implemented by indexing each phase with a set of marker points, which move in time according to their velocity interpolated from the velocity field at the nodal points using an affordable kernel approximation. The Phase tracking scheme was coupled with a marker point smoothing procedure to avoid clustering and dispersing of the marker points. The ability of proposed MDLSM multiphase method was tested using a set of benchmark problems indicating the appropriate accuracy of the proposed method to solve a series of classical multiphase flow problems from the literature.

Appendix A In this section, the stability of the proposed MDLSM method is investigated for a typical one-dimensional non-linear problem defined as:

∂U + U .(∇ U ) = ν ∇ 2U; ∂t

(A.1)

Using an explicit Euler method, as used in Eq. (23) for time discretization, leads to:



UIn+1 = UIn + t −UIn .(∇ UIn ) + ν

 ∇ 2UIn ; UIn = U (XI , tn )

(A.2)

For simplification, the equation is linearized as follows:



UIn+1 = UIn + t −Umax .(∇ UIn ) + ν

 ∇ 2UIn ;

(A.3)

where Umax is maximum value of U which is obviously more critical for stability conditions. Using MLS approximation, the semi-discretized equation can be written as:



UIn+1

=

UIn

+ t



−Umax

         ∂ 2 NJ (XI ) ∂ NJ (XI ) n UJ + ν UJn ; ∂x ∂ x2 J

 J

(A.4)

NJ (XI ) defines the value of shape function of J-th node at XI . Where NJ (XI ) defines the value of MLS shape function of the Jth node at XI . The MLS approximation satisfies the following conditions:



NJ (X ) = 1;

J

 ∂ NJ (X )  = NJ (X ) = 0; ∂x J

(A.5)

J

Defining the exact solution, U, as sum of the numerical approximate solution, Unumerical , and the error in the solution, u, (U = Unumerical + u), leads to the requirement that the solution error, u, should satisfy the semi-discredited equation of (A.5) since both the exact and numerical solutions satisfy this equation.



unI +1

=

unI

+ t



−Umax

 J

         ∂ 2 NJ (XI ) ∂ NJ (XI ) n uJ + ν unJ ; ∂x ∂ x2 J

(A.6)

Following the conventions used in von Neumann stability analysis, a typical term of the Fourier series is used to define the error as:

u(XI , tn ) = unI = eat eikXI ; u(XI , tn + t ) = unI +1 = ea(t+ t ) eikXI ;

216

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. A.1. Mapping between reference and physical domains that are discretized using uniform and symmetric arrangement of nodes.

u(XI + x, tn ) = unI+1 = eat eik(XI + x) ; u(XI − x, tn ) = unI−1 = eat eik(XI − x) ;

(A.7)

Substituting Eq. (A.7) into the semi-discretized equation of. (A.6), and assuming a uniform nodal distribution leads to:



unI +1

=

unI

+ t



−Umax

J=s 

 





at ik(XI +J x )

NJ (XI )e e

+

ν

J=−s

J=s 





at ik(XI +J x )

NJ (XI )e e

.

(A.8)

J=−s

where x is the nodal distance. Considering the uniform nodal distributions inside the support domain of the typical point XI (see Fig. A.1), we can write:

N0 = 0;

= −N and N = N ; ∀ J > 0 N−J J −J J

(A.9)

Substituting Eq. (A.9) into Eq. (A.8), dividing both sides of the resulting equation by unI , and using some mathematical simplifications, the following relation is arrived at for the amplification factor, G:



G = 1 + t



−Umax

m =s 

 



2N m (XI )i sin (mk x )

+

ν N 0 +

m=1

m =s 



2N 1 (XI ) cos (mk x )

;

(A.10)

m=1

With the amplification factor defined as G =

unI +1 . unI

For stability, the amplification factor should comply with the condition

of |G| ≤ 1. Consider that ξ N and ξ N respectively denote first and second derivatives of the MLS shape functions on a reference

x ξ between reference (ξ ) and physical (X) domain, domain with uniform nodal spacing of ξ . Using linear mapping X = ξ as shown in Fig. A.1, the following relations can be written:

N (XI ) = ξ N (ξI )

ξ ;

x



N (XI ) = ξ N (ξI )

Now Eq. (A.10) rewrites:

⎛



ξ −Umax

x



m =s 

ξ

x

2

;

(A.11)



⎞

2ξ N m i sin (mk x )

⎜ ⎟ ⎜ ⎟ m=1 ⎜ ⎟.     G = 1 + t ⎜ ⎟ m =s    ⎝+ ν ξ 2 N + ⎠ 2ξ N m cos (mk x ) 0

x

(A.12)

m=1

The stability condition now leads to:

& 



2 && &

ξ

ξ & & G = &1 + t −Umax λi + ν λ2 & ≤ 1.

x 1

x & &

(A.13)

The stability requirement strongly depends on the value of λ1 and λ2 defined as

λ1 =

m =s 

2ξ N m sin (mk x );

m=1

λ2 = N 0 +

m =s  m=1

(A.14) 2ξ N m cos (mk x );

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

217

Fig. A.2. First and second derivatives of MLS shape functions computed using different number of nodes inside the support domain (linear basis function).

Fig. A.3. Calculated |λ1 | and |λ2 | using different number of nodes inside the support domain (linear basis function).

For a specified number of nodes inside the support domain, λ1 and λ2 are constant and, therefore, the stability conditions of (|Umax | ν ) for convection-dominated and (|Umax |  ν ) for diffusion-dominated problems, leads to the requirement 2 x of t ≤ | U CF Lconvection |and t ≤ | vx CF Ldi f f usion |, respectively. The stability requirement, therefore, reduces to relations max

t ∼ o( x) and t ∼ o( x)2 for convection and diffusion dominated problems, respectively.

In what follows, the stability condition is investigated in terms of number of nodes located inside the support domain. Increasing the number of nodes in MLS approximation, decreases the absolute values of ξ N and ξ N at each typical node (see Fig. A.2), which in turn leads to a reduction of |λ1 | and |λ2 |. The stability region is, therefore, expected to expand with increasing number of nodes inside the support domain. To validate the claim, three different number of nodes are assumed inside the reference support domain with ξ = 0.1 and then ξ N and ξ N are computed assuming linear basis function. The results are schematically illustrated in Fig. A.2. |λ1 | and |λ2 | are also calculated in terms of θ (0 ≤ θ = mk x ≤ 2π ) and shown in Fig. A.3. As shown in Fig. A.3, both |λ1 | and |λ2 | decay with increasing the number of nodes inside the support domain. Therefore, it can be concluded that stability of the method can be enhanced by using more nodes in the support 2

ξ

ξ domain, especially for the convection-dominated problems (|Umax

x λ1 | |ν ( x ) λ2 |).

Appendix B A sensitivity analysis with respect to the value of the penalty parameter, α , is presented here for the dam break problem using 25 × 19 nodal distribution shown in Fig. 21. Fig. B.1 shows the penalty parameter, α , in terms of the corresponding condition number of the coefficient matrix. As shown in Fig. B.1, the condition number increases with increasing value of the penalty parameter.

218

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. B.1. Condition number of the coefficient matrix in terms of the penalty parameter.

Appendix C It can be seen from Eq. (24) that velocity, u, is proportional to ρ1˜ ∂∂xP . Since the pressure is continuous, it can be reai sonably assumed that the pressure gradient, ∂∂xP , does not sharply change around the interface [64]. In high-density ratio i

multiphase flows, however, the sharp variation of density across the interface leads to high variation of velocity on and around the interface schematically shown in Fig. C.1. This in turn may lead to a sharp variation in the velocity of the marker points and consequently the marker points distribution across the interface which may grow in time and lead to instability of the method. Different remedies are suggested in the literature to overcome this problem. The most common approach is to smooth the density on and around the interface by spatial weighted averaging [21,65,92,93]. This approach is implemented in the proposed MDLSM method via Eq. (58) where cubic spline kernel function is used as a weight function. The effectiveness of the density smoothing procedure strongly depends on the radius of support domain assigned to the nodal points, dw_node , (see Fig. 2). With increasing size of the support domains, smoothing of the averaged density, ρ˜ , is increased leading to improved stability of the method. To support this claim, the averaged density is calculated here for a simple one-dimensional case for different size of the support domain size, dw_node , and the results are illustrated in Fig. C.2 which shows ρ˜ is smoothed more for larger dw_node . However, a question now arises as to how much the accuracy is affected by increasing the size of the support domain. In what follows, the accuracy of the kernel approximation is investigated in terms of the support domain size. Similar to the Eq. (58), consider the following kernel approximation of a typical function, f(x),

f ( x0 ) =

 x−x  0 =dw ∫xx−x w | dw 0 | f (x )dx −x =−d  x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d w

0

0

w

;

(C.1)

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

219

Fig. C.1. Schematic variation of pressure gradient and velocity across interface.

where f(x0 ) is the approximated function, dw and x0 are radius and center of the support domain, respectively. Using Taylor expansion for f(x):

f ( x0 ) =

 x−x  0 =dw ∫xx−x w | dw 0 | −x =−d w

0



f ( x0 ) + f ( x0 ) ( x − x0 ) +

f (x0 )(x−x0 )2 2

 x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d 0

+ ...+

f (n ) (x0 )(x−x0 )n n!



+ . . . dx

;

(C.2)

w

For the above approximation to be exact, f(x0 ) = f(x0 ), the following conditions should hold:

 x−x  0 =dw ∫xx−x w | dw 0 | f (x0 )dx −x0 =−dw  x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d w

0

 x−x  0 =dw ∫xx−x w | dw 0 | −x =−d 0

= f ( x0 )

w

0 =dw ∫xx−x w −x =−d



w

0

f (n ) (x−x0 )n dx n! |x−x0 | dx dw



)

 x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d

w 0 = f (x0 );  x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d 0

(C.3)

w

 |x−x0 |  x−x0 =dw n f (n ) (x0 ) ∫x−x0 =−dw w dw (x − x0 ) dx = = 0;   x−x 0 =dw n! ∫xx−x w | dw 0 | dx −x =−d 0

∀n ≥ 1

(C.4)

w

Eq. (C.3) is automatically satisfied. Assuming n = 1 in Eq. (C.4) yields:



x−x '0 =dw

w x−x0 =−dw

|x − x0 | dw

(x − x0 )dx = 0;

(C.5)

220

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. C.2. Averaged density calculated for different sizes of the support domain.

Fig. C.3. Schematic graph for kernel approximations using radius dw and Dw .

This condition is always met by a symmetric, even, kernel function. Assuming n = 2, Eq. (C.4) leads to:



x−x '0 =dw

w x−x0 =−dw

|x − x0 | dw

(x − x0 )2 dx = 0;

(C.6)

Condition (C.6) is not necessarily met by the kernel approximation. Using Eqs. (C.6) and C.4, the error of the kernel approximation can be defined as follows:

errordw

 |x−x0 |  x−x0 =dw 2 f (x0 ) ∫x−x0 =−dw w dw (x − x0 ) dx = =C   x−x 0 =dw 2 ∫xx−x w | dw 0 | dx −x =−d 0

w



 x−x  0 =dw ∫xx−x w | dw 0 | (x − x0 )2 dx −x =−d  x−x  0 =dw ∫xx−x w | dw 0 | dx −x =−d

0

w

0

 ;

(C.7)

w

where C is a constant. To study the effect of the support domain size on the accuracy, consider kernel approximations defined on the support domains of radius dw and Dw with equal nodal distances, schematically shown in Fig. C.3. Similar to

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

221

Fig. C.4. Error norm of the kernel approximation versus support domain size.

the Eq. (C.7), the error can be defined as:



errorDw = C

 x−x  0 =Dw ∫xx−x w | Dw 0 | (x − x0 )2 dx −x0 =−Dw  x−x  0 =Dw ∫xx−x w | Dw 0 | dx −x0 =−Dw

 ;

(C.8)

using a change of variables defined as (x − x0 ) = (X − x0 )( Dd w ), Eq.(C.8) can be re-written as:

 ( X−x =d 0

w

 |X−x0 | 

X−x0 =−dw w dw ( X−x0 =Dw

errorDw = C

X−x0 =−Dw

hence

er rorDw =

 D 2 w

dw

er rordw ;

( 

w





2 X − x0 2 Ddww Ddww |X−x0 | Dw dX dw dw

) 

w

dX





=C



Dw 2 Dw dw dw

( X−x0 =dw





| 0| (X − x0 )2 dX X−x0 =−dw w dw   ( X−x Dw X−x0 =dw w | dw 0 | dX dw X−x0 =−d X−x

 ;

(C.9)

(C.10)

Eq. (C.10) shows that the error increases by order of two in terms of the support domain size. This claim is numerically investigated for a function defined as f(x) = sin(x), by calculating the slope of the kernel approximation error in terms of the support domain size for three different uniform nodal spacing, x, of 0.01, 0.05 and 0.1. For each nodal spacing size, the kernel approximation error was computed for three different support domain sizes of dw = 2 x, dw = 3 x and dw = 4 x which were then used to calculate the slopes as shown in Fig. C.4. Fig. C.4 clearly supports the conclusion analytically inferred from Eq. (C.10) that the kernel approximation error changes by an order of two in term of the support domain size. Appendix D An additional numerical example is solved here to provide a comparison between the results of the proposed MDLSM method and those of the well-known form of the SPH method presented by Brookshaw [27,94] which computes the second order derivatives explicitly, not requiring to solve any algebraic system of equations which in turns increases the computational cost. A quadratic PDE, ∂∂Tt = ∇ 2 T , is considered in a unit square which is solved on different set of nodal distributions. The Dirichlet type boundary condition is used on all boundaries for which the analytical solution is available. In Fig. D.1, the

222

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

Fig. D.1. Comparing the error norm and convergence rate of MDLSM and well-known SPH methods. Table D.1 CPU time for kernel and MLS approximations (using 2.50 GHz core E5200). Nodal distribution

Usual form of kernel approximation

MLS approximation

6×6 11 × 11 21 × 21

0.002505 0.008846 0.049473

0.009796 0.036130 0.144643

error norm of the MDLSM method at t = ) 0.1 (s) is compared with that of the well-known SPH method reported in the literN 2 1  ature [27]. The error norm is defined as (T analytical (Xi ) − T numerical (Xi )) . With N denoting the number of nodes. The N( i

result clearly shows that the proposed MDLSM method is more accurate than the SPH method. This may well be due to the fact that the MLS approximation used in the MDLSM enjoys higher order consistency compared to the kernel approximation used by SPH method. The CPU times required for shape function construction by the MLS and kernel approximations are presented in Table D.1 on different sets of uniform nodal distributions. Since the MLS approximation requires the solution of a linear algebraic system of equations (see Eqs. (7) and 10), its computational cost is obviously higher than that of the kernel approximation used in the SPH and MPS. References [1] G.-R. Liu, Meshfree Methods: Moving Beyond the Finite Element Method, CRC press, 2009. [2] S. Suleau, P. Bouillard, One-dimensional dispersion analysis for the element-free Galerkin method for the Helmholtz equation, Int. J. Numer. Methods Eng. 47 (20 0 0) 1169–1188. [3] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Astron. Soc. 181 (1977) 375–389. [4] L.B. Lucy, A numerical approach to the testing of the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [5] J.J. Monaghan, Smoothed particle hydrodynamics, Annu. Rev. Astron. Astrophys 30 (1992) 543–574.

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224 [6] [7] [8] [9] [10] [11] [12] [13] [14] [15] [16] [17] [18] [19] [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] [30] [31] [32] [33] [34] [35] [36] [37] [38] [39] [40] [41] [42] [43] [44] [45] [46] [47] [48] [49] [50] [51] [52] [53] [54]

223

G.-R. Liu, M.B. Liu, Smoothed Particle Hydrodynamics: A Meshfree Particle Method, World Scientific, 2003. J.P. Morris, P.J. Fox, Y. Zhu, Modeling low Reynolds number incompressible flows using SPH, J. Comput. Phys. 136 (1997) 214–226. J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994) 399–406. M. Khanpour, A.R. Zarrati, M. Kolahdoozan, A. Shakibaeinia, S. Jafarinik, Numerical modeling of free surface flow in hydraulic structures using smoothed particle hydrodynamics, Appl. Math. Model. 40 (2016) 9821–9834. M.S. Shadloo, A. Zainali, M. Yildiz, A. Suleman, A robust weakly compressible SPH method and its comparison with an incompressible SPH, Int. J. Numer. Methods Eng. 89 (2012) 939–956. R. Issa, E.S. Lee, D. Violeau, D.R. Laurence, Incompressible separated flows simulations with the smoothed particle hydrodynamics gridless method, Int. J. Numer. Methods Fluids 47 (2005) 1101–1106. E.-S. Lee, C. Moulinec, R. Xu, D. Violeau, D. Laurence, P. Stansby, Comparisons of weakly compressible and truly incompressible algorithms for the SPH mesh free particle method, J. Comput. Phys. 227 (2008) 8417–8436. S.J. Cummins, M. Rudman, An SPH projection method, J. Comput. Phys. 152 (1999) 584–607. S. Shao, E.Y. Lo, Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface, Adv. Water Resour. 26 (2003) 787–800. S. Koshizuka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dyn. J. 4 (1995) 29. S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nuclear Sci. Eng. 123 (1996) 421–434. A. Khayyer, H. Gotoh, Enhancement of stability and accuracy of the moving particle semi-implicit method, J. Comput. Phys. 230 (2011) 3093–3118. A. Khayyer, H. Gotoh, Modified moving particle semi-implicit methods for the prediction of 2D wave impact pressure, Coast. Eng. 56 (2009) 419–440. A. Khayyer, H. Gotoh, A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method, Appl. Ocean Res. 32 (2010) 124–131. M. Kondo, S. Koshizuka, Improvement of stability in moving particle semi-implicit method, Int. J. Numer. Methods Fluids 65 (2011) 638–654. A. Khayyer, H. Gotoh, Enhancement of performance and stability of MPS mesh-free particle method for multiphase flows characterized by high density ratios, J. Comput. Phys. 242 (2013) 211–233. A. Khayyer, H. Gotoh, Y. Shimizu, K. Gotoh, On enhancement of energy conservation properties of projection-based particle methods, Eur. J. Mech.-B/Fluids 66 (2017) 20–37. A. Shakibaeinia, Y.C. Jin, A weakly compressible MPS method for modeling of open-boundary free-surface flow, Int. J. Numer. Methods Fluids 63 (2010) 1208–1232. G.-R. Liu, Y.-T. Gu, An Introduction to Meshfree Methods and Their Programming, Springer Science & Business Media, 2005. N. Trask, M. Maxey, K. Kim, M. Perego, M.L. Parks, K. Yang, J. Xu, A scalable consistent second-order SPH solver for unsteady low Reynolds number flows, Comput. Methods Appl. Mech. Eng. 289 (2015) 155–178. H.F. Schwaiger, An implicit corrected SPH formulation for thermal diffusion with linear free surface boundary conditions, Int. J. Numer. Methods Eng. 75 (2008) 647–671. R. Fatehi, M. Manzari, Error estimation in smoothed particle hydrodynamics and a new scheme for second derivatives, Comput. Math. Appl. 61 (2011) 482–498. T. Belytschko, Y.Y. Lu, L. Gu, Element-free Galerkin methods, Int. J. Numer. Methods Eng. 37 (1994) 229–256. T. Belytschko, L. Gu, Y. Lu, Fracture and crack growth by element free Galerkin methods, Model. Simul. Mater. Sci. Eng. 2 (1994) 519. Y. Lu, T. Belytschko, L. Gu, A new implementation of the element free Galerkin method, Comput. Methods Appl. Mech. Eng. 113 (1994) 397–414. A. Huerta, Y. Vidal, P. Villon, Pseudo-divergence-free element free Galerkin method for incompressible fluid flow, Comput. Methods Appl. Mech. Eng. 193 (2004) 1119–1136. A.R. Firoozjaee, F. Farvizi, E. Hendi, Analysis of shallow water problems using element-free galerkin method, Int. J. Civil Eng. 15 (2017) 223–230. T. Zhang, X. Li, A variational multiscale interpolating element-free Galerkin method for convection-diffusion and stokes problems, Eng. Anal. Bound Elem. 82 (2017) 185–193. H. Yang, Y. He, Solving heat transfer problems with phase change via smoothed effective heat capacity and element-free Galerkin methods, Int. Commun. Heat Mass Transf. 37 (2010) 385–392. S.N. Atluri, T. Zhu, A new meshless local Petrov-Galerkin (MLPG) approach in computational mechanics, Comput. Mech. 22 (1998) 117–127. Q. Ma, Meshless local Petrov–Galerkin method for two-dimensional nonlinear water wave problems, J. Comput. Phys. 205 (2005) 611–625. N. Thamareerat, A. Luadsong, N. Aschariyaphotha, The meshless local Petrov–Galerkin method based on moving Kriging interpolation for solving the time fractional Navier–Stokes equations, SpringerPlus 5 (2016) 417. H. Lin, S. Atluri, The meshless local Petrov-Galerkin (MLPG) method for solving incompressible Navier-Stokes equations, CMES 2 (2001) 117–142. H. Lin, S. Atluri, Meshless local Petrov-Galerkin(MLPG) method for convection diffusion problems, CMES 1 (20 0 0) 45–60. Z.-J. Chen, Z.-Y. Li, W.-L. Xie, X.-H. Wu, A two-level variational multiscale meshless local Petrov–Galerkin (VMS-MLPG) method for convection-diffusion problems with large Peclet number, Comput. Fluids 164 (2018) 73–82. H. Arzani, M. Afshar, Solving Poisson’s equations by the discrete least square meshless method, WIT Trans. Model. Simul. 42 (2006) 23–31. M. Afshar, M. Naisipour, J. Amani, Node moving adaptive refinement strategy for planar elasticity problems using discrete least squares meshless method, Finite Element. Anal. Des. 47 (2011) 1315–1325. M. Afshar, J. Amani, M. Naisipour, A node enrichment adaptive refinement in discrete least squares meshless method for solution of elasticity problems, Eng. Anal. Bound. Elem. 36 (2012) 385–393. A.R. Firoozjaee, M.H. Afshar, Discrete least squares meshless method with sampling points for the solution of elliptic partial differential equations, Eng. Anal. Bound. Elem. 33 (2009) 83–92. M. Afshar, M. Lashckarbolok, G. Shobeyri, Collocated discrete least squares meshless (CDLSM) method for the solution of transient and steady-state hyperbolic problems, Int. J. Numer. Methods Fluids 60 (2009) 1055–1078. A.R. Firoozjaee, M.H. Afshar, Steady-state solution of incompressible Navier–Stokes equations using discrete least-squares meshless method, Int. J. Numer. Methods Fluids 67 (2011) 369–382. G. Shobeyri, M. Afshar, Simulating free surface problems using discrete least squares meshless method, Comput. Fluids 39 (2010) 461–470. G. Shobeyri, M. Afshar, Corrected discrete least-squares meshless method for simulating free surface flows, Eng. Anal. Bound. Elem. 36 (2012) 1581–1594. J. Amani, M. Afshar, M. Naisipour, Mixed discrete least squares meshless method for planar elasticity problems using regular and irregular nodal distributions, Eng. Anal. Bound. Elem. 36 (2012) 894–902. S. Faraji, M. Afshar, J. Amani, Mixed discrete least square meshless method for solution of quadratic partial differential equations, Scientia Iranica. Trans. A Civil Eng. 21 (2014) 492. S. Faraji, M. Kolahdoozan, M.H. Afshar, Mixed discrete least squares meshless method for solving the linear and non-linear propagation problems, Scientia Iranica 25 (2018) 565–578. S. Faraji, M. Kolahdoozan, M.H. Afshar, Collocated mixed discrete least squares meshless (CMDLSM) method for solving quadratic partial differential equations, Scientia Iranica 25 (2018) 20 0 0–2011. S.N. Kazeroni, M. Afshar, An adaptive node regeneration technique for the efficient solution of elasticity problems using MDLSM method, Eng. Anal. Bound. Elem. 50 (2015) 198–211. S.F. Gargari, M. Kolahdoozan, M.H. Afshar, Mixed discrete least squares meshfree method for solving the incompressible Navier–Stokes equations, Eng. Anal. Bound. Elem. 88 (2018) 64–79.

224

S. Faraji Gargari, M. Kolahdoozan and M.H. Afshar et al. / Applied Mathematical Modelling 76 (2019) 193–224

[55] S. Idelsohn, M. Mier-Torrecilla, E. Oñate, Multi-fluid flows with the particle finite element method, Comput. Methods Appl. Mech. Eng. 198 (2009) 2750–2767. [56] R. Reddy, R. Banerjee, GPU accelerated VOF based multiphase flow solver and its application to sprays, Comput. Fluids 117 (2015) 287–303. [57] J. Chessa, T. Belytschko, An extended finite element method for two-phase fluids, J. Appl. Mech. 70 (2003) 10–17. [58] Z.-B. Wang, R. Chen, H. Wang, Q. Liao, X. Zhu, S.-Z. Li, An overview of smoothed particle hydrodynamics for simulating multiphase flow, Appl. Math. Model. 40 (2016) 9625–9655. [59] S. Tiwari, J. Kuhnert, Modeling of two-phase flows with surface tension by finite pointset method (FPM), J. Comput. Appl. Math. 203 (2007) 376–386. [60] X.Y. Hu, N.A. Adams, A multi-phase SPH method for macroscopic and mesoscopic flows, J. Comput. Phys. 213 (2006) 844–861. [61] X. Hu, N.A. Adams, An incompressible multi-phase SPH method, J. Comput. Phys. 227 (2007) 264–278. [62] Q. Xiong, L. Deng, W. Wang, W. Ge, SPH method for two-fluid modeling of particle–fluid fluidization, Chem. Eng. Sci. 66 (2011) 1859–1865. [63] F. Yeganehdoust, M. Yaghoubi, H. Emdad, M. Ordoubadi, Numerical study of multiphase droplet dynamics and contact angles by smoothed particle hydrodynamics, Appl. Math. Model. 40 (2016) 8493–8512. [64] Z. Chen, Z. Zong, M. Liu, L. Zou, H. Li, C. Shu, An SPH model for multiphase flows with complex interfaces and large density differences, J. Comput. Phys. 283 (2015) 169–188. [65] M. Rezavand, M. Taeibi-Rahni, W. Rauch, An ISPH scheme for numerical simulation of multiphase flows with complex interfaces and high density ratios, Comput. Math. Appl. 75 (2018) 2658–2677. [66] M. Khanpour, A. Zarrati, M. Kolahdoozan, A. Shakibaeinia, S. Amirshahi, Mesh-free SPH modeling of sediment scouring and flushing, Comput. Fluids 129 (2016) 67–78. [67] A. Rahmat, M. Yildiz, A multiphase ISPH method for simulation of droplet coalescence and electro-coalescence, Int. J. Multiphase Flow 105 (2018) 32–44. [68] A. Shakibaeinia, Y.-C. Jin, MPS mesh-free particle method for multiphase flows, Comput. Methods Appl. Mech. Eng. 229 (2012) 13–26. [69] S. Lind, P. Stansby, High-order Eulerian incompressible smoothed particle hydrodynamics with transition to Lagrangian free-surface motion, J. Comput. Phys. 326 (2016) 290–311. [70] G. Fourtakas, P. Stansby, B. Rogers, S. Lind, An Eulerian–Lagrangian incompressible SPH formulation (ELI-SPH) connected with a sharp interface, Comput. Methods Appl. Mech. Eng. 329 (2018) 532–552. [71] J. Liu, S. Koshizuka, Y. Oka, A hybrid particle-mesh method for viscous, incompressible, multiphase flows, J. Comput. Phys. 202 (2005) 65–93. [72] P. Lancaster, K. Salkauskas, Surfaces generated by moving least squares methods, Math. Comput. 37 (1981) 141–158. [73] X. Li, S. Li, Analysis of the complex moving least squares approximation and the associated element-free Galerkin method, Appl. Math. Model. 47 (2017) 45–62. [74] A. Quarteroni, R. Sacco, F. Saleri, Numerical Mathematics, Springer Science & Business Media, 2010. [75] X. Li, S. Zhang, Y. Wang, H. Chen, Analysis and application of the element-free Galerkin method for nonlinear sine-Gordon and generalized Sinh-Gordon equations, Comput. Math. Appl. 71 (2016) 1655–1678. [76] R. Xu, P. Stansby, D. Laurence, Accuracy and stability in incompressible SPH (ISPH) based on the projection method and a new approach, J. Comput. Phys. 228 (2009) 6703–6725. [77] A. Mokos, B.D. Rogers, P.K. Stansby, A multi-phase particle shifting algorithm for SPH simulations of violent hydrodynamics with a large number of particles, J. Hydraul. Res. 55 (2017) 143–162. [78] P. Sun, A. Zhang, S. Marrone, F. Ming, An accurate and efficient SPH modeling of the water entry of circular cylinders, Appl. Ocean Res. 72 (2018) 60–75. [79] M.H. Afshar, A.R. Firoozjaee, Adaptive simulation of two dimensional hyperbolic problems by collocated discrete least squares meshless method, Comput. Fluids 39 (2010) 2030–2039. [80] M. Afshar, M. Lashckarbolok, Collocated discrete least-squares (CDLS) meshless method: error estimate and adaptive refinement, Int. J. Numer. Methods Fluids 56 (2008) 1909–1928. [81] T.E. Tezduyar, M. Behr, S. Mittal, J. Liou, A new strategy for finite element computations involving moving boundaries and interfaces—the deforming-spatial-domain/space-time procedure: II. Computation of free-surface flows, two-liquid flows, and flows with drifting cylinders, Comput. Methods Appl. Mech. Eng. 94 (1992) 353–371. [82] M. Cruchaga, D. Celentano, T. Tezduyar, A moving Lagrangian interface technique for flow computations over fixed meshes, Comput. Methods Appl. Mech. Eng. 191 (2001) 525–543. [83] M. Luo, C. Koh, M. Gao, W. Bai, A particle method for two-phase flows with large density difference, Int. J. Numer. Methods Eng. 103 (2015) 235–255. [84] J.C. Martin, W.J. Moyce, A. Price, C. Thornhill, Part IV. An experimental study of the collapse of liquid columns on a rigid horizontal plane, Phil. Trans. R. Soc. Lond. A 244 (1952) 312–324. [85] V.-T. Nguyen, D.-T. Vu, W.-G. Park, Y.-R. Jung, Numerical analysis of water impact forces using a dual-time pseudo-compressibility method and volume-of-fluid interface tracking algorithm, Comput. Fluids 103 (2014) 18–33. [86] R. Takahashi, H. Mamori, M. Yamamoto, Development and elaboration of numerical method for simulating gas–liquid–solid three-phase flows based on particle method, Int. J. Comut. Fluid Dyn. 30 (2016) 120–128. [87] A. Khayyer, H. Gotoh, Y. Shimizu, Comparative study on accuracy and conservation properties of two particle regularization schemes and proposal of an optimized particle shifting scheme in ISPH context, J. Comput. Phys. 332 (2017) 236–256. [88] T. Belytschko, Y. Krongauz, D. Organ, M. Fleming, P. Krysl, Meshless methods: an overview and recent developments, Comput. Methods Appl. Mech. Eng. 139 (1996) 3–47. [89] F. Macia, L.M. González, J.L. Cercos-Pita, A. Souto-Iglesias, A boundary integral SPH formulation: consistency and applications to ISPH and WCSPH, Progr. Theoret. Phys. 128 (2012) 439–462. [90] M. Liu, G. Liu, Smoothed particle hydrodynamics (SPH): an overview and recent developments, Arch. Comput. Methods Eng. 17 (2010) 25–76. [91] Z. Zhang, M. Liu, A decoupled finite particle method for modeling incompressible flows with free surfaces, Appl. Math. Model. 60 (2018) 606–633. [92] A. Colagrossi, M. Landrini, Numerical simulation of interfacial flows by smoothed particle hydrodynamics, J. Comput. Phys. 191 (2003) 448–475. [93] N. Grenier, M. Antuono, A. Colagrossi, D. Le Touzé, B. Alessandrini, An Hamiltonian interface SPH formulation for multi-fluid and free surface flows, J. Comput. Phys. 228 (2009) 8380–8393. [94] L. Brookshaw, A method of calculating radiative heat diffusion in particle simulations 6 (1985) 207–210.