MPS mesh-free particle method for multiphase flows

Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

Contents lists available at SciVerse ScienceDirect

Comput. Methods Appl. Mech. Engrg. journal homepage: www.elsevier.com/locate/cma

MPS mesh-free particle method for multiphase flows Ahmad Shakibaeinia, Yee-Chung Jin ⇑ Faculty of Engineering and Applied Science, University of Regina, Regina, Canada

a r t i c l e

i n f o

Article history: Received 19 May 2011 Received in revised form 21 February 2012 Accepted 19 March 2012 Available online 28 March 2012 Keywords: Mesh-free particle method MPS Multiphase flow Multi-viscosity flow Multi density flow

a b s t r a c t By treating the multiphase system as a multi-density multi-viscosity fluid, a straightforward model has been proposed in this paper based on the Moving Particle Semi-implicit (MPS) mesh-free particle method for incompressible multiphase flow. The weakly-compressible MPS (WC-MPS) formulation (developed by the authors) is used to solve a single set of equations for all of the phases. In the model, the multiphase forces are introduced in a straightforward way. Dealing with multi-viscosity systems, different methods for defining the viscosity, by which particles of different phases interact, is examined. To evaluate the accuracy of each of these methods, a stratified multi-viscosity Poiseuille flow test case is used and model results are compared with the analytical solution. The results show that selection of this interaction viscosity has an important role in the accuracy of the results. The model is then validated and applied to two basic hydrodynamic instability cases (Rayleigh–Taylor and Kelvin–Helmholtz instabilities). The results are then compared to a mesh-based model (with volume of fluid interface tracing method) to examine the ability of the model to deal with the multi-density systems. Comparisons show the model has reasonable accuracy. The results of this work offer the potential of modeling of multiphase incompressible immiscible systems in an extensive range of conditions using MPS. Ó 2012 Elsevier B.V. All rights reserved.

1. Introduction From the point of view of the computational frame, the numerical methods for handling multiphase systems can be generally classified into three main groups. The first group is the fully Eulerian approaches such as Volume-Of-Fluid (VOF) method [1] which is characterized by a mesh that is either stationary or adaptive to the interfaces. The second group is the hybrid Eulerian–Lagrangian approaches such as Particle-In-Cell (PIC) method [2] in which a Lagrangian frame is used for interface tracing where the field variables are computed on a Eulerian frame. The last group is the fully Lagrangian approaches such as Smooth Particle Hydrodynamics (SPH) [3–5], Dissipative Particle Dynamics (DPD) [6], and Moving Particle Semi-implicit (MPS) [7,8] methods. Eulerian methods with a stationary mesh have problems with maintaining a sharp interface between two phases while those with a deforming mesh have problems with adaptation of the mesh to follow the highly deformed and fragmented interfaces [9]. Hybrid Eulerian–Lagrangian methods have difficulties with coupling of the Eulerian and Lagrangian frames, as well as numerical inaccuracy that can occur due to the coupling process. The mesh-free particle (Lagrangian) methods belong to the group of fully Lagrangian approaches and are known to be capable of dealing with any interface deformations and fragmentations. ⇑ Corresponding author. E-mail address: [email protected] (Y.-C. Jin). 0045-7825/$ - see front matter Ó 2012 Elsevier B.V. All rights reserved. http://dx.doi.org/10.1016/j.cma.2012.03.013

One of the most useful characteristics of the particle methods in multiphase systems is that the behavior of phases arises from interaction of the particles alone and is irrelevant to the position of interface. The Smoothed-Particle Hydrodynamics (SPH) method was originally developed for astrophysical applications and was later expanded for applications in solid and fluid mechanics. It works by dividing the fluid into a set of discrete elements, referred as particles and is based on integral representation of quantities and spatial derivatives. A number of SPH models have been developed for multiphase systems (e.g. [10–14]). Moving Particle Semi-implicit method (MPS), which is based on Taylor series expansion, was originally developed for fluid mechanics applications by Koshizuka et al. [7] and Koshizuka and Oka [8]. The MPS method is similar to the SPH method. The fundamental difference of MPS and SPH is that the MPS method of spatial derivative approximation is solely based on a local weighted averaging process without incorporating the gradient of a kernel function. Since its development, MPS has been applied to a wide range of fluid mechanics applications, including modeling of dam breaks [7], breaking waves [15], jet breakup [16], hydraulic jump formation [17], and flow over the sills and in trenches [18]. The sub particlescale (SPS) turbulence model was proposed by Gotoh et al. [19] for simulation of turbulent flow using MPS. Khayyer and Gotoh [20] and Shakibainia and Jin [17] showed that the original MPS formula does not guarantee the conservation of linear and angular momentum and proposed some modification of the MPS pressure gradient calculation as well as the deployed kernel function. The artificial

14

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

Fig. 1. Particle interaction conceptual model.

pressure fluctuation in the MPS was addressed by Khayyer and Gotoh [20] and Shakibaeinia and Jin [17]. A multiphase MPS method was introduced by Gotoh and Fredsøe [21] for solid–liquid and by Ikari et al. [22] for gas–liquid two phase flow by solving two sets of equations for different phases and adding the interaction force to the momentum equation. A hybrid particle-mesh method was proposed by Liu et al. [9] for incompressible viscous multiphase flow, in which one phase is represented by the moving particles and the other one is defined by a stationary mesh. Shakibaeinia and Jin [17] proposed the weakly compressible MPS method (WCMPS) for modeling of the incompressible fluids and showed that this method not only improves the MPS artificial fluctuations but also slightly increases the model efficiency compare to standard MPS (fully compressible MPS). This paper then aims to develop a new model based on WC-MPS for multiphase incompressible systems. Here, a straightforward model in which the multiphase system is treated as a multi-viscosity multi-density fluid is employed, and only a single set of equations is solved for the whole of the flow field. In the present method, the multiphase forces are automatically calculated by retaining all terms of the governing equations, such as the viscosity gradients. A fractional steps algorithm is used for time splitting. The pressure is calculated using the WC-MPS method. Particle motion is modified using a pair-wise collision algorithm. For calculation of the shear stress between two different-phase particles, an interaction viscosity must be defined. In this study, different methods of defining of the interaction viscosity are evaluated by applying the model to the stratified single-density, multi-viscosity Poiseuille flow test case and comparing the results with the analytical solution. The final model is validated using two of the basic hydrodynamic instabilities, the Rayleigh–Taylor (multi-viscosity, multi-density) and Kelvin–Helmholtz (inviscid, multi-density) instability cases. The results are then compared to a mesh-based VOF model (Fluent6.3 software is used for this purpose).

Fig. 2. Third order polynomial spiky kernel function.

q

Du ¼ rp þ rðlr  uÞ þ f; Dt

where u is the flow velocity, q is the density, p is the pressure , l is the dynamic viscosity, and f represents the body forces (per unit volume). Needless to say that, in Lagrangian systems, there is no convective acceleration term in the left hand side of the momentum equation. The motion of each particle is simply calculated by Dr/ Dt = u, (r being the position vector). Spatial derivatives in the governing equations need to be approximated on the particle domain. For calculation of spatial derivatives as well as interpolation of physical quantities, the MPS particle interaction model is used. 3. Numerical method 3.1. MPS appoximations As a mesh-free method MPS interaction is built on a set of disordered points in a continuum without a grid or mesh. The particle i, interacts with other particles in its vicinity, j, covered with a kernel (weight) function W(rij, re), where rij = |rj  ri| is the distance between particle i and j (either in the same phase or different phases) and re is the radius of the interaction area around each particle (Fig. 1). The kernel function is considered to be a smoothing function of physical quantities around particles. This article employs a third order polynomial spiky kernel function (Fig. 2) proposed by Shakibaeinia and Jin [17] and written as:

( Wðr; r e Þ ¼

2. Governing equations The multiphase system is treated as a multi-density multi-viscosity fluid, hence a single set of equations can be used for all of the phases. The governing equation is expressed by the continuity and momentum equations. In the Lagrangian system, they can be written as: Continuity:

Dq þ qðr  uÞ ¼ 0: Dt Momentum:

ð1Þ

ð2Þ

ð1  r=re Þ3

0 6 r=r e < 1;

0

r=r e P 1:

ð3Þ

A particle interacts with a finite number of adjacent particles located within a distance of re. A dimensionless parameter, namely, particle number density, n, defined in Eq. (4) represents the density of particles around a specific particle. Real fluid density, q, is defined by Eq. (5).

hnii ¼

P Wðr ij ; re Þ;

ð4Þ

j–i

P j–i mj Wðr ij ; r e Þ hqii ¼ R ; Wðr; r e Þdv V

ð5Þ

15

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

where mj is mass of particle j. Operator h i is kernel approximation (kernel smoothing or interpolating) operator. The denominator in Eq. (5) is the integral of kernel function in the whole region (except the region occupied by the central particle i). In MPS the smoothed value of the physical property of f on particle i is given by weight averaging of that property over the neighboring particles as:

P 1 P j–i ðfj Wðr ij ; r e ÞÞ hf ii ¼ P ðfj Wðrij ; r e ÞÞ: ¼ hnii j–i j–i ðWðr ij ; r e ÞÞ

ð6Þ

The gradient vector of a physical quantity such as p on a particle i can be achieved by weighted averaging of gradient vectors between particle i and its neighboring particles j as:

  d P pj  pi hrpii ¼ 0 eij Wðr ij ; re Þ ; n j–i r ij

rij eij ¼ : rij

ð7Þ

In case of pressure gradient, Koshizuka et al. [15] replaced the pi in ^i defined by: the above equation with p

^i ¼ minðpi ; pj Þ; p j2J

J ¼ fj : Wðrij ; r e Þ–0g:

ð8Þ

This means that the minimum value is selected among the neighboring particles with the distance of re from the target particle. According to Koshizuka et al. [15], this replacement improves the stability of the model by ensuring the inter-particle repulsive force. The pressure gradient term is then written as:

hrpii ¼

  ^i d P pj  p eij Wðr ij ; re Þ : n0 j–i r ij

ð9Þ

pnþ1 ¼ i

qi c20 ððhn ii =n0 Þc  1Þ c

ð13Þ

in which, the typical value of c = 7 is used; c0 is sound speed in reference density. Since, using the real sound speed for the fluid would result in extremely small time steps, an artificial smaller sound speed is usually used. Morris et al. [25] proposed an estimation of the numerical sound speed by considering the balance between the pressure, viscous force and body force as:

c20 ¼ max

! V 2b tV b FL ; ; ; d dL d

d¼

Dq

q0

ð14Þ

;

where F is the magnitude of the body force, Vb is the bulk velocity (i.e., average velocity magnitude), d is density variation, and L is characteristic length scale. To keep density variation of fluid less than one percent of the reference density, the Mach number must be smaller than 0.1. In other words, sound speed must be higher than ten times the maximum fluid velocity [26]. Proportion of particle number density is calculated in the prediction step, n⁄, and initial average particle number density, n0, is used in the equation of state to compute pressures in the new time step. In this article, the term MPS is used for the WC-MPS method. Note that, although the WC-MPS is not a fully incompressible model, the use of a divergence-free viscous term in the momentum equations still can be reasonable. It is due to this fact that, the order of magnitude of the term including the divergence of the velocity vector is very small comparing to the divergence-free tem. 3.3. Multiphase solver

The divergence formula of velocity vector u at the particle i is given by:

  d P uj  ui hr  uii ¼ 0  eij Wðr ij ; r e Þ : n j–i r ij

ð10Þ

In the MPS, Laplacian formula is given by weighted averaging of distribution of a physical quantity from a particle to its neighboring particles [15]. With a similar concept, the vector Laplacian on each particle is also calculated as:

! 2d P ðuj  ui Þ hr uii ¼ 0 Wðr ij ; re Þ : n j–i r 2ij 2

ð11Þ

In the MPS Laplacian formula, r 2ij is replaced by a representative parameter k, using the weighted averaging of r 2ij as:

2d P ððuj  ui ÞWðr ij ; r e ÞÞ; kn0 j–i P P k ¼ hr 2ij ii ¼ ðr 2ij Wðrij ; r e ÞÞ= ð Wðrij ; r e ÞÞ:

hr2 uii ¼

j–i

The method of this study, assume the multiphase system to be a multi-density multi-viscosity system, hence a single set of equation can be used for all of the phases (similar to the approach used in DPD by Visser et al. [27] and in SPH by Hu and Adams [33]). The density (or mass) and viscosity are inherent properties of each particle and are initially assigned to each of particles based on the phase which that particle belongs to. These inherent properties remain unchanged during the simulations (the exemption is for the non-Newtonian fluids where the viscosity of fluid can change). The model needs to take care of multiphase forces arising from the density and viscosity discontinuities. The density difference between phases explicitly affects the momentum terms. The density appears in the momentum equation and also is used for calculation of the pressure as follows:

pnþ1 ¼ qi ai ; i ð12Þ

j–i

Transferring of a physical property from a particle to its neighboring particles using a Gaussian function results in a variant increase of k. The applied k to the MPS Laplacian formula make the increase of variance equal to that of the analytical solution.

ai ¼

c20

c

ððhn ii =n0 Þc  1Þ:

ð15Þ

For the calculation of the pressure of a certain particle, the density of that particle is directly used. By replacing the value of pressure from (15) into the MPS gradient formula, the pressure gradient term

ρ

Phase I i

Interface

Phase II

3.2. WC-MPS pressure calculation Here in this paper the weakly compressible MPS [17] is used for pressure calculation. Fluid is assumed to be weakly compressible, making it possible to use an explicit equation (equation of state) to determine pressure in each time step instead of solving an implicit relation (Poisson equation). By keeping the compressibility value very small, the fluid is treated as nearly incompressible fluid. The Tait’s equation of state, described by Batchelor [23] and Monaghan [24] and modified by Shakibaeinia and Jin [17] for the MPS is used.

r Transition region

Transition region

Fig. 3. Sketch smooth density used for pressure calculation.

16

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

will be directly related to the ratio of density of neighboring particles j to that of the target particle i (qj/qi) as:

1

qi

hrpii ¼

  d P ðqj =qi Þaj  ai e Wðr ; r Þ : e ij ij n0 j–i rij

ð16Þ

For calculation of velocity of a specific particle, its density appears in the equations directly. Therefore the density difference is automatically taken care of. Due to the direct relation between the density and pressure, the density discontinuity near the interface can cause a discontinuity in the pressure field. To avoid this problem, one can use the smooth value hqii of density instead of the real density for the pressure calculation.

pnþ1 ¼ hqii ai : i

ð17Þ

Far from the interface the smooth value of the density is equal to each of fluids’ density, while within a transition region around the interface the smooth value of the density varies between the fluids’ density (see Fig. 3). When the viscosity of phases in the multiphase system is different, a viscosity discontinuity happens near the interfaces. The existence of a different phase particle, within the support area of a specific particle leads to a non-zero viscosity gradient across the interface and creates a multiphase force due to the viscosity difference. The viscous term can be written as:

rs ¼ rðlr  uÞ ¼ rl r  u þ lr2 u:

ð18Þ

In the area around interfaces the viscosity gradient is not equal to zero. Using the MPS gradient and divergence formula, the term containing the viscosity gradient can be approximated by:     2  P lj  li d P uj  ui hrlðr  uÞii ¼ 0  eij Wðrij ;r e Þ eij Wðrij ;re Þ : n r ij rij j–i j–i

ð19Þ Limiting the value of compressibility, the order of magnitude this term rl (u) is much smaller than that of the second term (lr2u). Although the value of this term may not be significant compared to the second part of the viscous term, its appearance near the interfaces can have some small effects on the momentum exchange in this area.; therefore this term is not neglected here. For the second part of the viscous term, the MPS approximation can be written as:

hlr2 uii ¼

Fig. 4. Sketch of the pair-wise collision.

2d P ðl ðuj  ui ÞWðr ij ; r e ÞÞ; kn0 j–i ij

ð20Þ

respectively. The above equation returns the arithmetic mean for h = 1 and harmonic mean for h = 1. Application of the interaction viscosity in calculation of the viscous term automatically adds the multiphase forces (arising from the viscosity difference) to the momentum equation. 3.4. Pair-wise particle collision One of the problems with particle methods such as MPS is interpenetration which causes particle clustering. For the MPS model in this study, a collision algorithm is used to modify the velocity field and avoid their penetration. The algorithm was originally used by Koshizuka et al. [15] for the single-phase systems. The idea of collision of fluid particles in MPS comes from the conservation of linear momentum in collision of two granular gas or solid particles in the normal direction [28]. The pair-wise particle collision is characterized by the particle mass and distance and mechanism of dissipation of mechanical energy during collisions. In MPS the collision algorithm is used for particles whose distance is less than h times of the average particle distance (h called dimensionless collision distance). The relative velocity of the particles, uin j = uj  ui, can be decomposed into the normal part, uij , and the t tangential part, uij , (see Fig. 4). Here only the normal part is considered and assume that the tangential velocity is not affected by the collision. Collision of particle i and j with density of qi and qj is then given by:

8   > ^ i ¼ ui  q1 ð1 þ eÞ qqþi qqj unij :u ^ j ¼ uj þ q1 ð1 þ eÞ qqþi qqj unij j

where lij is an interaction viscosity between particles i and j. When particle i is interacting with a same phase particle j, it is obvious that the real viscosity of that fluid must be used (lij = li = lj). However, when particle i with a viscosity li is interacting with a different phase particle j with viscosity of li, none of the fluid viscosities can be used as the interacting viscosity lij because this could underestimate or overestimate the friction between particles i and j. The unknown value of lij lays somewhere between li and lj and is difficult to determine. Modeling of multi-viscosity fluids using DPD, Visser et al. [27] showed that for the cases where the position of the interface is unknown, the mean value of friction factors of interacting fluids provides a good accuracy in the velocity profile. To evaluate the effects of different values for the interaction viscosity, in this study a two-phase viscosity stratified flow (i.e. multiphase Poiseuille flow) test case is used. We examine the extreme values lij = li and lij = lj and compare results with the mean values calculated by: 1

lij ¼ ððlhi þ lhj Þ=2Þh :

ð21Þ

With h approaching the positive and negative infinities, the value of the interaction viscosity approaches to max{li, lj } and min{li, lj},

i

where unij ¼ ðuij  eij Þeij ;

ð22Þ

j

^ is the velocity vector after the collision, and e is collision where u ratio. This relation is only applied to the particles that are approaching each other, or in other words, the value of (uij  eij) must be positive. According to the dissipative nature of particle collisions the kinetic energy of the relative motion decreases and, hence,

d2

μ2

Fluid 2 U0

p0 d1

μ1

u(y)

y x

p0 + Δp

Fluid 1

l Fig. 5. Two-phase Poiseuille flow between two parallel walls.

17

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

Fig. 6. Comparison the longitudinal velocity profile computed by the MPS with different values of the interaction viscosities lij.

Fig. 7. MPS solution of the multiphase Poiseuille flow for M = 1/4.

Fig. 8. Comparison the longitudinal velocity profile computed by the MPS (discrete points) and the analytical solution (solid line).

0 6 e 6 1. For e = 1 the total kinetic energy is conserved, that is, the collision occurs elastically. If e = 0, the complete energy of the relative motion is lost, which leads the particles to be bound together after a collision. In this study e = 0.5 is used. If the particles have the same density (qi = qj) the equation can be simplified :

(

^ i ¼ ui  12 ð1 þ e Þunij ; u ^ j ¼ uj þ 12 ð1 þ e Þunij : u

ð23Þ

In general, for the particle i the velocity after collision is given by:

^ i ¼ ui  u

1 P

qi

j2J

!

qi qj n ð1 þ eÞ u ; J ¼ fj : r ij < hDlg: qi þ qj ij

ð24Þ

In the momentum exchanges between particles, the total linear and angular momentum are preserved and therefore, this pairwise collision algorithm does not disturb the satisfaction of the momentum equation. Eq. (24) is comparable to what is used in

18

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

time step the velocity field is modified using the pair-wise collision algorithm. The position of a particle in the new time step is calculated based on these velocities. 4. Numerical examples To illustrate the ability and accuracy of the proposed model to deal with different multiphase systems, several standard test cases are used. The first case, which is a two-layer multiphase Poiseuille flow, is used to evaluate the effect of application of different values as the interaction viscosity and to examine the ability of the model to simulate the viscosity stratification. Furthermore, two types of fundamental hydrodynamics instability cases are simulated in this study. Rayleigh–Taylor and Kelvin–Helmholtz instabilities under different conditions are modeled to examine the ability of the model to handle the density stratification and to predict deformations and fragmentations in these interfacial instabilities. 4.1. Multiphase Poiseuille flow

Fig. 9. Configuration of the Rayleigh–Taylor instability test.

the XSPH method [29] to avoid particle interpenetration and to move particles in a velocity closer to the average velocity of the neighboring particles.

Poiseuille flow is a sensitive test case to examine the accuracy of the model. This test has been used to examine the accuracy of the calculated velocity profile, especially near the interface, using different methods of estimation of the interaction viscosity. The horizontal stratified two-phase flow between two stationary infinite parallel plates is considered (Fig. 5). The lower layer is considered as the first fluid and the top layer as the second one. The gravity and surface tension is neglected. Both of the fluid movements are subject to pressure difference Dp between two sections with distance of l. The interface has a distance of d1 from the bottom plate, and d2 from the top plate. The velocity profile predicted by the model is compared to the analytical solution. The analytical velocity field is given by [30]:

8 < uðyÞ ¼ U 0 ð1 þ

2 y Md 2 Mðd þd Þ d2

: uðyÞ ¼ U ð1 þ Md 2 0 2 

3.5. Solution algorithm

d

The governing equations can be solved using a fractional step method splitting each time step into two pseudo steps of prediction and correction. For particle i, the velocity in the new time step is given as the summation of the velocity calculated in the prediction step and the correction step.

uikþ1 ¼ ui þ u0i :

ð25Þ

Body forces as well the viscous terms are explicitly calculated in the prediction step. The predicted velocity is given by:

Dui ¼

Dt

qi

ðf i þ rli ðr  ui Þ þ lij r2 ui Þ:

ð26Þ

Particles are simply moved according to their predicted velocity as follows:

ri ¼ ri þ ui Dt:

ð27Þ

Based on the predicted particle position, the particle number density n⁄ is calculated. The difference of n⁄ and the initial value of the particle number density n0 provides a method for calculation of the pressures. The velocity is corrected using the pressure gradient term as:

u0 ¼

Dt

qi

rpi :

ð28Þ

The summation of the predicted velocities and their correction gives the velocities in the new time step. Before preparing for the next



Mþd ð y Þ2 Þ 2 Mðd þd Þ d2

Mþd ð y Þ2 Þ 2 d þd d2

d1 6 y 6 0; 0 < y 6 d2

ð29Þ

v ðyÞ ¼ 0: U0 is interfacial velocity given by: 

U0 ¼ 

1 dp d1 d2 d þ 1 : 2 dx l2 M þ d

ð30Þ

Relevant dimensionless parameters include the ratio of the dynamic viscosities M, the ratio of interfaces distances to the walls d⁄ and dimensionless velocities u⁄ and v⁄ are:

M¼

ui ¼ uki þ Dui ;

y þd d2



l1 d1 ðu; v Þ  ; d ¼ ; ðu ; v  Þ ¼ ; U0 l2 d2

ð31Þ

where u and v are the corresponding dimensional quantities. The test domain is defined by 0.50 6 x 6 0.50, 0.50 6 y 6 0.50 and is subject to a pressure gradient Dp/l = 0.50. The interface is located at y = 0 (d1 = d2 = 0.50). The lower and the top fluids both have the density of q1 = q2 = 2 in the absence of gravitational acceleration. The lower fluid dynamic viscosity is set to l1 = 0.05 while that of the top one is determine by l2 = l1/M (M = 1/8, 1/4 and 1/2). This case is totally numerical and all quantities are dimensionless. The periodic boundary condition is assigned to the inflow and the outflow boundaries, and a no-slip boundary condition pertains for the lower and upper walls. The initial particle distance Dl = 1/80 resulting a total of 6400 fluid particles. The model results (after achieving a relatively steady state) are compared to analytical solution. The MPS velocities is extracted from discrete particles located on x = 0 ± Dl/2. Fig. 6 illustrates the comparison of the analytical solution and the normalized velocity profile from the MPS model with a viscosity ratio of M = 1/4 using

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

19

Fig. 10. Evolution of R–T instability for At = 1/2, viscosity = 0.010 (tRT = 1.1, 2.2, 3.3, 4.4).

four different methods of estimation of the interaction viscosity lij (li, lj, arithmetic mean and harmonic mean of li and lj) in the Eq. (20). As it can be seen in this figure, using lij = li, the velocity profile near the interface not only overestimates the analytical value, but also has been linearized. Conversely, lij = lj underestimates the velocity profile. By applying the arithmetic mean, the both the shape and value of the velocity approaches to those of the analytical solution, but the velocity has still been underestimated. The best results are achieved by using the harmonic mean of viscosities as the interaction viscosity. Fig. 7 also shows the particle representation of the phases and the normalized velocity, resulted from MPS with a harmonic mean of viscosities near the interface. The transverse component of the velocity is insignificant; therefore, the interface shape is kept almost as flat as the initial condition. Since the longitudinal velocity is almost uniform along with the flow direction, the cross section from which the velocity profile is extracted does not affect the results. Fig. 8 compares the steady normalized velocity profile for three different dynamic viscosity ratios, M = 1/2, 1/4 and 1/8. These results are achieved by using the harmonic mean of the fluids viscosities as the

interaction viscosity. It can be seen that the MPS results show a reasonable agreement with the analytical solution though some underestimations are evident near the velocity peak point. For the lower m values, this discrepancy is more apparent. Here the reliability of these results can be evaluated by L2 norm error, defined as:

eL2 ¼

N P i¼1

ðDUÞ2i

N P i¼1

ðUÞ2i

1=2 ð32Þ

in which DU is the deviation of numerical velocity from the analytical velocity, U, and N is the number particles from which the velocity profiles are extracted. The L2 norm error calculated from the above formula are 0.007, 0.017 and 0.036 for M = 1/2, 1/4 and 1/8, respectively. The value of the errors seems to be not very significant; therefore, the model has reasonable accuracy. 4.2. Rayleigh–Taylor instability The Rayleigh–Taylor (R–T) instability provides an excellent opening into the mathematical study of stability because of the

20

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

exceptionally simple nature of the base state. Consider two completely plane-parallel layers of immiscible and unequal-densities fluid, the heavier on top of the lighter one. Both fluids are subject to the acceleration g, normal to their planar interface in the direction of heavy to light fluid. The equilibrium between the two layers becomes unstable due to small perturbations or disturbances. An unstable disturbance will grow exponentially in time (linear regime) and lead to the release of potential energy as the heavier material moves down under the gravitational field and the lighter material is displaced upwards. The lighter fluid forms a bubble that eventually breaks up into a two-phase mixing zone (a mushroom shape) that grows algebraically in time (non-linear regime). The problem is characterized by Atwood number, At, and viscosity ratio, m, defined by:

At ¼

qh  ql l ; M ¼ h; qh þ ql ll

ð33Þ

where h and l are denoted for the heavy and the light fluid, respectively. The perturbation has the wavelength of k (the wave number of k = 2p/k) and amplitude of g0. According to the linear growth regime the perturbation grows exponentially and the maximum interfacial displacement g is given by:

gðtÞ ¼ g0 expðxtÞ;

ð34Þ

where x is the growth rate. Neglecting the surface tension, the R–T instability growth rate for viscous fluids is simply given by the positive root of the following [9]:

x2 þ 2tk2 x  Atgk ¼ 0:

ð35Þ 1/2

Fig. 11. A snapshot of R–T instability results (At = 0.50, viscosity = 0.010 and tRT = 4.4).

Scaled velocity and characteristic time, are defined by (kg) and (k/g)1/2, respectively. Other relevant dimensionless parameters include dimensionless time t RT and the maximum interfacial displacement g⁄:

Fig. 12. Comparison of the maximum interfacial displacement in MPS, VOF and theoretical (from the linear regime).

Fig. 13. Comparison of MPS results with the VOF interface (solid line).

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

21

Fig. 14. Evolution of R–T instability for (a) At = 1/2, viscosity = 0.008; (b) At = 1/3, viscosity = 0.010.

pffiffiffiffiffiffiffiffi t RT ¼ t= k=g ;

g ¼ g=k:

ð36Þ

In this study, a rectangular domain, defined by 0.25 6 x 6 0.25 and 0.50 6 y 6 0.50 is used (Fig. 9). At |x| = 0.25 and |y| = 0.50, the noslip boundary condition is used. A constant gravitational acceleration g = 10 is added. A single mode perturbation wave with a wavelength of k = 0.5 and amplitude of 0.025 is added to the system. Fluids have the same kinematic viscosities and different densities. The density of the lighter fluid is set to be one while that of the heavier fluid is calculated based on the selected Atwood number. This is a totally numerical test case and all quantities are dimensionless. Here, the Atwood numbers 1/2 and 1/3 with two different viscosities of m = 0.010 and m = 0.008, are selected. This test problem is used to demonstrate the model’s ability to handle the density stratification and to evolve a linear perturbation into a nonlinear hydrodynamic turbulence. For the MPS model the initial particle distance Dl = 1/200 results in a total of 20,000 fluid particles. The results are compared to a mesh-based finite volume model in which VOF is used for modeling of the interfaces. A mesh resolution of 100 by 200 (cell size of 1/200) for VOF model is comparable to the resolution which is used in the particle method. Fig. 10 illustrates the results for the growth of a single perturbation driven by the R–T instability and its pressure and the y-velocity

for an Atwood number of At = 1/2 and kinematic viscosity of m = 0.010. This figure demonstrates the ability of the model to show the common features of the R–T instability. As the perturbation grows the lighter fluid moves upward and forms a bubble. After approximately t RT = 1.1 the lighter fluid start to form a mushroom-shape feature. The mushroom size grows vertically and horizontally in response to the gravitational force and the vortex roll, respectively (see Fig. 11). Eventually the interface breaks up into a two-phase mixing zone that grows in time. The model is able to properly reflect the R–T features for different conditions. The velocity and pressure has some fluctuations. Unphysical pressure fluctuation, which is intensified near the interface (Fig. 10), is one of the problems with MPS. Due to the density difference, the frequency of the pressure wave is different in light and heavy fluid; therefore, at the interface, some extra noise is visible due to collision of the different mode pressure pulses. Fig. 11 shows an enlarged snapshot of the particles position as well as the velocity vectors. It illustrates a relatively smooth interface represented by the particles although there is some phase jump (i.e., small scale interpenetration of one phase into the other) near the interface due to the unphysical fluctuations as well as the lack of the surface tension. The vortex rolls and the momentum exchange along and perpendicular to the interfaces are also shown in this figure.

22

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

Fig. 15. Results of the K–H instability test for a density ratio of 2:1.

The maximum interfacial displacement achieved by the MPS, VOF and the analytical solution (from the linear regime) are shown in Fig. 12. At early times (before t RT = 1.1) the MPS results are quite compatible with the VOF and the analytical solution which validates the effect of the viscosity and density difference on the growth of the instability. Late time behavior of the R–T growth does not follow the linear regime therefore, the figure show a lower value for the displacement of the interfaces achieved by the numerical methods. The MPS and VOF have close interfacial displacement, although, it seems that as the time proceeded the MPS shows a faster growth. Fig. 13 compares the R–T features that resulted from the MPS model with the VOF interface. The shape of the interface is compatible in these models, although some discrepancies are observed. The reason of these discrepancies can be the fundamental difference in the methods and algorithms. To examine the ability of the model under different conditions the model is applied to two other cases with a lower viscosity and a lower Atwood number, respectively. In Fig. 14(a) the kinematic viscosity is reduced to m = 0.008 while, in Fig. 14 (b) the lower Atwood number of At = 1/3 is used. The lower viscosity increases the width of the mushroom feature and also intensifies the vortices as it can be seen in downward-moving part of the interface. As expected a lower Atwood number decreases the size and growth rate of the bubble and mushroom features.

the horizontal direction has been examined. The domain is a periodic box defined by 0.5 6 x 6 0.5, 0.5 6 y 6 0.5. The fluid has a density of q1 = 2 for |y| 6 0.25 and q2 = 1 for |y| > 0.25 in absence of gravitational acceleration. The initial u velocity is:



u ¼ þU 0 jyj 6 0:25; u ¼ U 0 jyj > 0:25:

ð37Þ

An initial velocity perturbation is imposed to the system as:

v ¼ g0

       

2 px 2px þ p exp ð10ðy  025ÞÞ2  sin exp ð10ðy þ 0:25ÞÞ2 ; sin k k

ð38Þ where g0 is the perturbation amplitude. This small perturbation wave, initiate an instability that start to grow. Initially, this growth is linear, but as the time proceeds, the wave turns into a non-linear form. In this case, the linear growth rate (x) of the K–H instability is given by [31]:

x¼

2p ðq1 q2 Þ1=2 ð2U 0 Þ : k q1 þ q2

ð39Þ

And the characteristic growth time is given by:

sKH ¼

2p ðq1 þ q2 Þk : ¼ w ðq1 q2 Þ1=2 ð2U 0 Þ

4.3. Kelvin–Helmholtz instability

The dimensionless time is defined as:

The Kelvin–Helmholtz (K–H) instability occurs in the presence of sufficient velocity difference across the interface between two fluids or velocity shear within a continuous fluid. In this study, instability of a two-phase incompressible inviscid fluid flowing in

tKH ¼

t

sKH

:

ð40Þ

ð41Þ

For the case selected here, k = 0.25, g0 = 0.05 and U0 = 0.5. For density ratio of 2:1, the linear growth rate x = 11.85 and the character-

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

23

Fig. 16. Results of the K–H instability test for a density ratio of 5:1.

istic growth time sKH = 0.53, while, for the density ratio of 5:1, growth rate x = 9.38 and sKH = 0.67. This test problem is used to demonstrate the model’s ability to evolve a linear perturbation into nonlinear hydrodynamic turbulence. For the MPS model the initial particle distance Dl = 1/250 resulted in a total of 62,500 fluid particles. Ten times of the horizontal velocity (U0) is selected as the numerical value of the sound speed, the higher values increase the frequency of unphysical fluctuations which may result in a numerical instability in the results. The K–H instability is a particularly challenging test for the MPS because the unphysical pressure fluctuations can cause an unphysical and growing perturbation. Furthermore, the velocity due to particle unphysical fluctuation can approach the numerical value of the sound speed which can take away the imposed velocity perturbation. The results are compared to a mesh-based finite volume model in which VOF method is used for modeling interfaces. A mesh resolution of 250 by 250 (cell size of 1/250) is comparable to the resolution is used in the particle method. Fig. 15 shows the evolution of the K–H instability test using a density ratio of 2:1. The initial velocity perturbation creates a waveform interfaces with growing amplitude. The shear velocity and the non-zero curvature of the interface lead to formation of a ripple shape interfaces at approximately t KH = 0.5. The flow of one fluid around the other fluid creates a centrifugal force that amplifies the ripples and finally creates some vortex rolls at t KH = 1.0. Fig. 16 shows the evolution of the system for the density ration of 5:1. This figure illustrates the ability of the model to evolve K–H instability-like features in the higher density ratios where the multiphase standard SPH has some difficulties [12,13].

To make the mesh-free results comparable to mesh-based results, the noise introduced by using discrete particles must be removed. This can be done by mapping the discrete values on a uniform mesh with the same resolution as the mesh-based model. Fig. 17 compares the MPS density variation (mapped on a mesh) with the mesh-based VOF results for the density ratio of 2:1. The K–H instability-like features are similar, and the vortex rolls are comparable to the mesh-based results. However, the MPS wave length seems to be greater than the VOF. Even by mapping the MPS results on a mesh, the fluid interfaces are still not as smooth as the mesh-based interfaces calculated using VOF. Part of this discrepancy is due to the fact that the MPS artificial pressure fluctuation causes some noise on the velocity field [18]. As it has been proved in Fig. 18 the MPS interfacial fluctuations is irrelevant to the multiphase model because such interfacial fluctuations also can be seen for the case with density ratio of 1:1. These interfacial fluctuations are most likely due to some unphysical pressure fluctuations arising from original MPS formulation. Another reason for the discrepancy is that, the VOF method removes the natural interfacial fluctuations smaller than the mesh size in the procedure of interpolation of the volume of fraction near the interface. Fig. 19 magnifies a single vortex roll in Fig. 17 at t KH = 1.0 that resulted from the model and compares it with the mesh-based (VOF) results. At the top, the figure shows the comparison of the configuration of phases represented by MPS particles and VOF mesh. The position, shape and the size of the vortex roll is similar in these models. As can be seen, the sharpness of the interfaces is preserved by the MPS model. Since VOF has difficulty with maintenance of a sharp interface [9], it is reasonable to see a transition area for the VOF interface. There are some fluctuations at the MPS interface

24

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

while the VOF interface remains relatively smooth. The fluctuation in the MPS velocity is also visible in the configuration of the velocity vectors and u and v velocity magnitude (Fig. 19). In terms of the velocity magnitude the MPS and VOF results appear similar.

5. Summary and conclusion In this study a simple and straightforward method for modeling of multiphase incompressible immiscible systems, based on the

WC-MPS formulation, was applied. The method is characterized by treating the multiphase system as a multi-viscosity and multidensity system and solving a single set of momentum equations for all of phases. The multiphase forces due to the viscosity and density differences are automatically considered. The equations were then discretized and solved using the MPS method. Different methods of defining the interaction viscosity were examined. The results of the multiphase Poiseuille flow test case showed that the harmonic mean of fluid viscosities for interaction of two unlike fluid particles, gives a velocity profile compatible with the analyt-

Fig. 17. Comparison of MPS and mesh-based (VOF) results for the K–H instability test with a density ratio of 2:1.

Fig. 18. A snapshot of MPS volume fraction and pressure field for a density ratios of 2:1 and 1:1 at tKH = 1.0.

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

25

Fig. 19. Comparison of MPS and mesh-based (VOF) results of a single wave test for a density ratio of 2:1 at t KH = 1.0.

ical solution. Several standard test cases, with different conditions, were used to examine the ability of the model to deal with the multiphase systems. The results showed that the model can successfully predict the basic hydrodynamics instabilities and their related features in the high and low density rations. The behavior of the phases is directly related to the interaction between particles

and each of particles can have a certain density and viscosity therefore the number of phase that can be modeled is as many as the number of particles are used. As the extension of this work, the surface tension model can be added to the method to make it capable of dealing with the interfaces in which the surface tension is important. This can be done by implementation of the continuum

26

A. Shakibaeinia, Y.-C. Jin / Comput. Methods Appl. Mech. Engrg. 229–232 (2012) 13–26

surface force (CFS) [32] or the new method propose by Adami et al. [33] for the SPH applications. Acknowledgment This research was supported by the Natural Science and Engineering Research Council of Canada. References [1] C.W. Hirt, B.D. Nichols, Volume of fluid (VOF) method for the dynamics of free boundaries, J. Comput. Phys. 3 (9) (1981) 201. [2] F.H. Harlow, The particle-in-cell computing method for fluid dynamics, Methods Comput. Phys. 3 (1964) 319–343. [3] R.A. Gingold, J.J. Monaghan, Smoothed particle hydrodynamics: theory and application to non-spherical stars, Mon. Not. R. Soc. 181 (1977) 375–389. [4] L.B. Lucy, Numerical approach to testing the fission hypothesis, Astron. J. 82 (1977) 1013–1024. [5] J.J. Monaghan, An introduction to SPH, Comput. Phys. Commun. 48 (1988) 89– 96. [6] P.J. Hoogerbrugge, J.M.V.A. Koelman, Simulating microscopic hydrodynamics phenomena with dissipative particle dynamics, Europhys. Lett. 19 (1992) 155– 160. [7] S. Koshizuka, H. Tamako, Y. Oka, A particle method for incompressible viscous flow with fluid fragmentation, Comput. Fluid Dynam. J. 4 (1995) 29–46. [8] S. Koshizuka, Y. Oka, Moving-particle semi-implicit method for fragmentation of incompressible fluid, Nucl. Sci. Engrg. 123 (1996) 421–434. [9] J. Liu, S. Koshizuka, Y. Oka, A hybrid particle-mesh method for viscous, incompressible, multiphase flows, J. Comput. Phys. 202 (2005) 65–93. [10] J.J. Monaghan, A. Kocharyan, SPH simulation of multi-phase flow, Comput. Phys. Commun. 87 (1995) 225–235. [11] X.Y. Hu, N.A. Adams, A multi-phase SPH method for macroscopic and mesoscopic flows, J. Comput. Phys. 213 (2006) 844. [12] O. Agertz, B. Moore, J. Stadel, D. Potter, F. Miniati, J. Read, L. Mayer, A. Gawryszczak, A. Kravtsov, Å. Nordlund, F. Pearce, V. Quilis, D. Rudd, V. Springel, J. Stone, E. Tasker, R. Teyssier, J. Wadsley, R. Walder, Fundamental differences between SPH and grid methods, Mon. Not. Roy. Astron. Soc. 380 (2007) 963–978. [13] D.J. Price, Modelling discontinuities and Kelvin–Helmholtz instabilities in SPH, J. Comput. Phys. 227 (2008) 10040–10057. [14] 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 (22) (2009) 8380–8393.

[15] S. Koshizuka, A. Nobe, Y. Oka, Numerical analysis of breaking waves using the moving particle semi-implicit method, Int. J. Numer. Methods Fluids 26 (1998) 751. [16] K. Shibata, S. Koshizuka, Y. Oka, Numerical analysis of jet breakup behavior using particle method, J. Nucl. Sci. Technol. 41 (2004) 715–722. [17] A. Shakibaeinia, Y.C. Jin, A weakly compressible MPS method for simulation open-boundary free-surface flow, Int. J. Numer. Methods Fluids 63 (2010). 1208–123. [18] A. Shakibaeinia, Y.C. Jin, MPS-based mesh-free particle method for modeling open-channel flows, J. Hydraul. Engrg. 137 (2011) 1375–1385. [19] H. Gotoh, T. Shibahara, T. Sakai, Sub-particle-scale turbulence model for the MPS method-Lagrangian flow model for hydraulic engineering, Comput. Fluid Dynam. J. 9 (2001) 339–347. [20] A. Khayyer, H. Gotoh, Modified moving particle semi-implicit methods for the prediction of 2D wave impact pressure, Coast. Engrg. 56 (2009) 419–440. [21] H. Gotoh, J. Fredsøe, Lagrangian two-phase flow model of the settling behavior of fine sediment dumped into water, in: Proceedings of ICCE (2000) Sydney, Australia, pp. 3906–3919. [22] H. Ikari, H. Gotoh, T. Sakai, Simulation of wave breaking by the particle method with liquid–gas two-phase-flow model, Ann. J. Coast. Engrg. (JSCE) 111–115 (in Japanese). [23] G.K. Batchelor, An Introduction to Fluid Dynamics, Cambridge University Press, Cambridge, UK, 1967. [24] J.J. Monaghan, Simulating free surface flows with SPH, J. Comput. Phys. 110 (1994) 399–406. [25] J.P. Morris, P.J. Fox, Y. Zhu, Modeling low Reynolds number incompressible flows using SPH, J. Comput. Phys. 136 (1997) 214–226. [26] R.A. Dalrymple, B.D. Rogers, Numerical modeling of water waves with the SPH method, Coast. Engrg. 53 (2–3) (2006) 141–147. [27] D.C. Visser, H.C.J. Hoefsloot, P.D. Iedema, Modelling multi-viscosity systems with dissipative particle dynamics, J. Comput. Phys. 214 (2) (2006) 491–504. [28] T. Pöschel, T. Schwager, Computational Granular Dynamics: Models and Algorithms, Springer, Berlin, 2005. [29] J.J. Monaghan, On the problem of penetration in particle methods, J. Comput. Phys. 82 (1989) 1–15. [30] Q. Cao, K. Sarkar, A.K. Prasad, Direct numerical simulations of two-layer viscosity-stratified flow, Int. J. Multiphase Flow 30 (2004) 1485–1508. [31] S. Chandrasekhar, Hydrodynamic and Hydromagnetic Stability, Oxford University Press, 1961. [32] J.U. Brackbill, D.B. Kothe, C. Zemach, A continuum method for modeling surface tension, J. Comput. Phys. 100 (2) (1992) 335–354. [33] S. Adami, X.Y. Hu, N.A. Adams, A new surface-tension formulation for multiphase SPH using a reproducing divergence approximation, J. Comput. Phys. 229 (13) (2010) 5011–5021.