Journal of Applied Sciences ISSN 1872-5654 science alertJournal of Applied Sciences 8 (21); 3817-3826, 2008 ISSN 1812-5654 © 2008 Asian Network for Scientific Information Modeling Two-Phase Flows Using SPH Method A, Valizadeh, 'M, Shafieefar, *1.J. Monaghan and 'S.A.A. Salehi Neyshaboori “Department of Civil Engineering, Tarbiat Modares University, Tehran, Iran *School of Mathematical Sciences, Monash University, Clayton, Australia Abstract: The basic equations governing the free surface flow of a fluid are considered and approximated by using a semi-compressible Smoothed Particle Hydrodynamics (SPH) method which is a grid-less Lagrangian approach. Tt is found that the standard SPH method, which has been successfully applied to simulate the transient free surface flows, could be used for simulation of multiphase flows with high density ratios, by doing a correction on mass of interface particles. The SPH modeling is shown to provide a promising tool to predict the two liquid characteristics of free surface flows. To show the validity of the present method mimerical examples on bubble rising in a heavier fluid and luck-exchange flow are considered and a good agreement is observed, Key words: Multi-phase flows, free surface flows, smoothed particle hydrodynamics INTRODUCTION ‘The Smoothed Particle Hydrodynamics (SPH) method is a powerful particle approach which was originally developed for the study of astrophysies (Monaghan and Gingold, 1983, Monaghan, 1988) and later employed to study the wave overtopping over coastal structures (Dalrymple and Rogers, 2006, Gomez-Gesteira and Dalrymple, 2004), In the simulations of incompressible flows using the SPH model, the incompressiblity is achieved by employing an equation of state and the fluid is assumed to be slightly compressible, which is called the weakly compressible SPH method. In SPE method, materials are approximated by particles that are free to move around rather than by fixed grids or meshes. The particles are basically moving interpolation points that carry with them physical properties, such as the mass, velocity, viscosity, density of the fluid and any other properties that are relevant. The inter-particle forces are calculated by smoothing the information from nearby particles in a way that ensures that the resultant particle motion is consistent with the motion of a corresponding real fluid, as determined by the Navier- Stokes equations. Despite its fully Lagrangian property, when the standard formulation of SPH is applied to multi-phase flows only small density differences are permitted between the considered phases because it is implicitly assumed that the density gradient is much smaller than that of the smoothing Kernel (Monaghan, 1994; Monaghan and Kocharayan, 1995; Monaghan et al, 1999), As a remedy, Ritchie and Thomas (2001) suggest a summation of the particle-averaged pressure, not density, to handle large density gradients. However, their method does not satisfy mass conservation. Colagrossi_ and Landini (2003) modified the approximation form of spatial derivatives to diminish the effects of large density difference across the interface. Nevertheless, since the density summation is replaced by @ non-conservative density evolution equation mass conservation is not satisfied either. Although the conservation errors are decreased somewhat by a special density re-initialization approach, they may accumulate and affect the flow behavior considerably in long time computations, Hu and Adams (2006) employed a particle smoothing function in Which neighboring particles only contribute to the specific volume but not to the density. The resulting algorithm resolves a density discontinuity at a phase interface naturally and satisfies mass conservation since a density summation equation was employed, this method ‘was suffering from clumping and clustering, so Hu and Adams (2007) developed an incompressible SPH model which solves the Poisson equation for pressure by a projection formulation. This model is a good model for two phase flows. The scope of this study is finding a quick and simple way for correction the standard semi- compressible SPH model, for high density difference ‘multiphase flows, The essential step is that by changing the mass of particles which have interaction with particles from a different phase, the undesirable force between two particles is controlled, this simple idea improves the results very well and comparison model's results with previous numerical, analytical and experimental results supports this claim, Corresponding Author: Mehdi Shafieefar, Department of Civil Engineering, Tarbiat Modares University, P.O, Box 14115/143, Tehran, Iran Tel: +98-21-8288 3318 Fax; +98-21-8800 6544 3817J. Applied Sei., 8 (21): 3817-3826, 2008 In the present study, first the general concepts of the standard SPH modeling for incompressible flows are given. On this ground, the present SPH-implementation is, introduced and described in details and derived graphs for mass coefficient in a case that light fluid is surrounded by a heavy fluid (eg, bubble rising in water) are presented. Using these graphs the lock-exchange flow for different fluids is simulated and compared with experimental works and other numerical models MATERIALS AND METHODS, SPHis a particle based method for modeling coupled fluid flows, solid structure deformation and heat transfer. ‘The particles represent blobs of discretised fluid or solids, that move around in response to the fluid or solid stresses, produced by the interaction with other particles Importantly, SPH does not use any fixed grids or meshes, to track the fluid and calculate the fluid velocities Fonmally, SPH is a Lagrangian continuum method for solving systems of partial differential equations. The fluid (or solid) is discretised and the properties of each of these fluidsolid elements are attributed to their centers, whieh are then interpreted as a particle SPH uses an interpolation kernel to smooth the values of any information held by the particles to give smooth continuous interpolated fields. Governing equations: The almost incompressible medium. water is not made totally incompressible but more compressible in SPH. The full mass conservation equation has to be used: Dp a Pee pve a) pb vp yust, ab) (lc) The Eq, la-e give a time derivative of the density, velocity and the position of fluid particles. p, u, v, x, t represent the density, velocity, viscosity, position and time, respectively and f, is a body force such as gravity With a relation between the pressure and the density it is no longer necessary to solve the implicit Poisson equation for the pressure, Numerical model: In SPH, the fhuid field is represented as a collection of N neighbor particles, interacting each other through evolution equations of the general form: nf -t (2a) 118) 70 K-07, (a Du, (P+P. PH sm | PL, + 2b) mPa TMM PS (20) Bt where. p, u. P,, m, x, are the density, velocity, pressure, ‘mass and position of particle i fis a body force, I, is the viscous term and W, = Wr, h) is the smoothing function 1, = tis the position vector from particle j to particle i and h ‘is the smoothing length. By neglecting the almospheric pressure, the state equation between pressure and density is as follow (Batchelor, 1974): 38 where, P, and p, ate the reference pressure and density of fluid, y = 7, Py = o'py/y and c is the used speed of sound and instead of using the real sound speed which needs using a very small time step for the stability reasons, an artificial sound speed of ¢ « 10 ig is used in SPH. With this value for speed of sound, the maximum relative density differences ate small in order of ~ 1% (Monaghan, 1994), Now the calculation time stays reasonable and water is only slightly compressible in SPH, In multi-phase flows P,P, earel Y may be different for each phase, The quartic spline kemel which is used here is in this form 25-94-s0s-9)'+1005-9!' gzqcos Wi 2s-qt-sas-qi ossqsls Hees @s-a)4 13sq525 ° 2554 @ where, q = ry/h,, 1 is the absolute distance between particles i and j and. The spatial derivative of kemel funetion for two dimension problems is calculated as: () The choice of the smoothing length h is determining the number of interactions for each particle. When h is too small, there are not enough particles nearby to interact with, giving a low accuracy. When hi is too big, local 3818J. Applied Sei., 8 (21): 3817-3826, 2008 properties are smeared out too much, this gives a low accuracy again and the calculation becomes slow. Here an adaptive kemel funetion is used by defining a variable smooth length, using the following equation (Monaghan, 2005), (6) where, d is the number of dimensions, According to this, formula the smooth length is increasing by decreasing in density and vice-versa. So, in places such as interface between two phases or near the boundaries, kemel function adopts itself to the conditions. With a smoothing length around 1.4 times the initial particle spacing, Ap, the circle of influence has a radius of 3.5 and 28 Ap for quartic and cubic spline functions, respectively. With this smoothing length every particle has 36 and 20 neighbors in two dimensions for mentioned functions, respectively. That isa good optimum for both accuraey and calculation speed I, is the viscous term and is defined as an artificial viscous stress term, Different methods for defining this, term have been produced by scientists (Monaghan, 1994; Cleary, 1998; Gomez-Gesteira and Dalrymple, 2004). One of the prevalent definitions which is used in this study is, (Monaghan, 2005) (a) (7b) (0) B= (0;+P)/? is the average demsity. Uy . this is the velocity difference between two particles. @=(@ +a)/2 is the constant for defining the viscosity difference between two particles and a 1s chosen normally between limits 1 > c > 0, B is a coefficient which in all simulations is equal to 2. This artificial viscosity introduces. both shear and bulk viscosity, but with neglectable changes in the density it is almost entirely shear viscosity. This approach to model viscosity guarantees that momentum, is transferred from the particle with the highest velocity to the particle with lowest velocity Colagrosi and Landrini (2003) showed that the averaging of the densities helps ensure that the free surfaces are smooth and physically acceptable and Dalrymple and Rogers (2006) performed this filtering every 40 time steps by reinitializing the density of each particle according to: sm, oe Em, yi @) In this study, the reinitialisation of density is done every 40 time steps for each phase separately, it prevents of density deviations and keeps the density of particles near to reference density. The XSPH correction for movement of particle is included according to Eq, 9, it means an average between the particle's velocity and the average field velocity from all neighbor particles is used for moving a partie in Eq, 2c. This smoothing on velocity prevents particles penetration. Also for consistency reasons, the smoothed velocity has to be used in continuity equation, It should be noted that the new smoothed velocity is used just in Eq, 2a and c and in the momentum equation the non- smoothed velocity is used (Monaghan, 1989), o where, ¢ is a constant (0 <€ <1). With € =1 the particles are moved with the field velocity, with © = 0 the particles are moved with their own velocity, normally € = 0.5 ‘Time integration: By using the SPH approximation, three differential equations should be integrated in time: Pe _ ax (10a) 20 Fax) = FON) D (ux (0b) (0X9) (10e) ‘These differential equations are integrated in time by using the predictor-corrector time integration which predicts the velocity, density and particle position at the intermediate time steps, as follows: Predictor step: yh gM, og? f-U2 x02) a) oP oo Beaty aby (le) 3819J. Applied Sei., 8 (21): 3817-3826, 2008 where, A’, A" and A"? denote the value of variable A at the beginning, intermediate and at the end of (n+l )th time step, respectively. Corrector step: re) (2b) hap? , gn ney (12¢) With these predictors the time integration is second order acourate. The time step, &t, is calculated according to the Courant condition modified for the existence of viscosity: a3) &t, isthe initial time step, h is the smoothing length, 0 is, a coefficient between zero and one, normally @ = 0.5 and v, is defined in Eq, 7b. Boundaries: Boundaries in SPH can be modeled by virtual particles that characterize the system limits, Basically, three different types of particles can be distinguished: (@) Boundary particles with repulsive forces: ‘The closed boundary is modeled with a single row of boundary particles with mass, fixed density equal to fluid particles density and no pressure, Their positions are fixed or extemal imposed during the calculation (Fig. 1a). A repulsive boundary force is, used to stop approaching fluid particles in analogy with the forces among molecules (Monaghan, 1994, 2008) (b) Quasi fluid boundary particles: One or more rows of fixed fluid particles are placed at the boundary. No estimates of boundary forees or dieetions are needed; the boundaty particles just like fluid particles verify equations of continuity and state, except that their position is fixed or externally imposed in time ‘They build up pressure just like normal fluid particles, this is how they prevent particle penetration. Fluid particles near the boundary have interaction with a double staggered tow of quasi-fluid particles. When a slip boundary is needed these quasi fluid particles are used in the calculation of the viscous stress terms, for a free-slip boundary they are not used (Crespo et al, 2007). w Moving bound particles Fd particles \ ere OPtSeocccce CReccccoccce Cocccccce Eotesoscee eoce, eteccccecet| Me Cer erent Secencccce evenscceseccssogsesessoes Fixed boundry particles Gost panicles Fig. 1: Modeling of closed boundaries in SPH method, (a) boundary particles with repulsive force and (b) ghost particles, dark particles are fluid and gray particles are ghost particles (c) Ghost particles: For every particle within a distance shorter than the kernel support radius (e.g,, 2 h, for cubic spline function) from the boundary, a virtual (ghost) particle is mirrored outside the boundary every time step (Randles and Libersky, 1996). Both particles have same mass, density and pressure, but velocities normal to boundary are opposite. The velocities parallel to boundary can be the same in sleep boundaries and opposite for non-slip ones (Fig, 1b). In this research the closed boundaries are simulated by ghost particles; the main advantage of this method is that boundary is compatible with different phases and it can adopt itself with fluid particles. Because we are working with particles with different characteristics, especially with different mass and densities, the boundary is adaptive with particles. By using fixed particles for closed boundary, either quasi-fMuid particles or repulsive force particles, there will be difficult to simulate an adaptive boundary which can be compatible with different 3820J. Applied Sei., 8 (21): 3817-3826, 2008 phases, Also, using particles with repulsive forces impose an undesirable kinetic energy in system which needs some initial time steps for damping this energy. For problems which start with an instability such as bubble rising or dam breaking, initial damping of the kinetic energy is difficult but by using ghost particles the amount of initial kinetic energy is small and there is no need to initial damping Inthis study, the bottom and wall friction are defined by selecting the velocity of ghost particles in direction parallel to boundary in this form Ug = Cutty aa) where, u, and uy are the parallel to boundary velocity components of ghost and fluid particles, respectively and C, is a coefficient which defines the amount of friction and slip or non-slip boundary (-1 om) Time ee) Fig. 7 all cases, Exp: Experimental results of Grobelbauer e SPH: Current SPH model, p,~ 1000 kg m~ anda Time (see) Comparison of the numerical and experimental results, (a) light front and (b) dense front. For lock-exchange, in sval, (1993), FE: Finite Element model of Btienne ef al. (2005), (PPPs os ® # Gribeibauer er al. (1993) () © Geivetbauer eal (1993) b 1B Etienne ea, (2008) asf & Flieaeer al (2008) A SPH 4 SPH — 3 “ery 2 - E04. je 04 00 00. 80 05 To 80 05 v0 i Fig. & Comparison of the numerical and experimental res f= 1000 kg m™ and «= (p.-P)/ps the results of Btienne ef al, (2005), For the dense front, « nearly constant shift in time between the calculated and measured front arrivals is observed (Fig. 7b) similar to Etienne ef al (2005) numerical results. Possible explanations for this time shift according to reasoning of ftienne et al. (2005) are either a large numerical inconsistency in time around t =O or a time lag in the measurements. Thus, we are left to suppose that there is either a uniform time lag in the measured arrival time of the front, which may be due to a detection problem, or that the time shift is due to the opening of the gate. In Fig. 8a and b, the experimental and mumerical Froude numbers ¥; - u, ai’ and ¥y~ UyieH of light and dense fronts for different density’ ratios are compared Since the Froude mumber only accounts for the established velocity of the front, the shift in front arrival between mumerical and experimental results has no effect It is shown from Fig, 8a that, for the light front, the ‘SPH model is in close agreement with the experiments of Grdbelbauer ef al. (1993) and also with numerical results, wults for Froude numbers, (a) Light front and (b) dense front, of Etienne e¢ al. (2005), nnmerical and experimental results concerning Fr, deviate from the straight line p° i Figure 8b shows that for the dense front, the SPH model as well as the model of Etienne ef al. (2005), for Reynolds numbers corresponding, to the experiments of Grobelbauer et al. (1993), fit closely the experiments. The data points can be closely fitted by a power law of the form Fi18via-a-p5) In Fig, 9 the shape of lock-exchange flow and nondimensional vorticity for B = 0.11, at two stages in the flow development is shown. The vorticity of each particle is defined as (Monaghan, 1992): i as) Fn The general shape of vorticity is approximately similar but is not in excellent agreement. This may be improved by doing on vorticity formula in SPH and by increasing the resolution of SPH model some correction 3824J. Applied Sei., 8 (21): 3817-3826, 2008 Fig. 9 Comparison of vorticity atthe interface of two fluids, top to bottom, vorticity from Fitienne ef al. (2005), vorticity from SPH and the density profile, respectively CONCLUSION ‘A multi-phase free surface flow implementation of ‘Smoothed+-Particle Hydrodynamics (SPH) is presented and a new correction on standard SPH is provided. In this model without defining the interface between two phases cor using complicated correction for standard SPH method and just by changing the mass of particles which have interaction with particles from other phases according to derived grap from bobble rising in a heavier fluid and an engineering, judgment, itis possible to model the violent two-phase free surface fluids with density ratios (py/p) up to 1000, Numerical examples are investigated and compared with analytic solutions, previous results and experiments, The results suggest that this method can be faithfully applied to multi-phase flows. Since its construction is based on the standard SPH method the involved corrections are simple to implement and suitable for straightforward extension to three dimensions. ACKNOWLEDGMENT ‘The authors gratefully acknowledge Iranian Water Resources Management Co, support through grant ENVI-84120, REFERENCES, Batchelor, GK. 1974. An Introduction to Fluid Mechanics, 4th Edn., Cambridge University Press, Cambridge, UK., ISBN: 0521663962, pp: 635. Cleary, P.W., 1998, Modelling confined multi-material heat and mass flows using SPH. Applied Math. Model, 22: 981-993, Colegrossi, A. and M, Landrini, 2003, Numerical simulation of interfacial flows by smoothed particle hydrodynamics, J. Comput, Phys,, 191: 448-475. Crespo, A.J.C., M. Gomez-Gesteria and RA. Dalrymple, 2007. Dynamic boundary particles in SPH. 2nd Intemational Workshop on SPHERIC-Smoothed Particle Hydrodynamics. May 23rd-25th, Madrid, Spain, pp: 152-155. Dalrymple, RA. and B.D. Rogers, 2006. Numerical modeling of water waves with the SPH method, Coastal Eng., 53: 141-147 Davies, RM. and G. Taylor, 1950. The mechanies of large bubbles rising through extended liquids and through liquids in tubes, Proc. R, Soc, London Series A, Math. Phys. Sci,, 200; 375-390 Etienne, J, EJ. Hopfinger andP. Saramito, 2005. Numerical simulations of high density ratio lock-exchange flows. Phys. Fluids, 17: 036601-036601-12, Gomez-Gesteira, M. and R.A. Dalrymple, 2004, Using a 3D SPH method for wave impact on a tall structure J. Waterway Port Coast, Ocean Eng,, 130: 63-69. Grobelbauer, HP., LK. Fannelop and R.E, Britter, 1993, ‘The propagation of intrusion fronts of high density ratios. J. Fluid Mech., 250: 669-687 Hu, XY. and N.A. Adams, 2006. A multi-phase SPH method for macroscopic and mesoscopic flows, J. Comput. Phys., 213: 844-861 3825J. Applied Sei., 8 (21): 3817-3826, 2008 Hu, X.Y. andN.A, Adams, 2007, An incompressible multi- phase SPH method. J. Comput, Phys, 227: 264-278, Liu, Z. and Y, Zheng, 2006. PIV study of bubble rising behavior. Powder Technol., 168: 10-20. Monaghan, J.J and R.A, Gingold, 1983. Shock simulation by the particle method SPH. J. Comput. Phys., 52: 374.389. Monaghan, J.., 1988. An introduction to SPH, Comput Phys. Commun. 48: 89-96, Monaghan, J.J,, 1989. On the problem of penetration in particle methods, J. Comp, Phys., 82: 1-15 Monaghan, J.J., 1992, Smoothed patticle hydrodynamics ‘Ann, Rey. Astronom, Astrophys., 30: 543-574 Monaghan, J.J., 1994. Simulating free surface flows with SPH, J, Comput, Phys., 10: 399-406, ‘Monaghan, J.J. and A. Kocharyan, 1995. SPH simulation of multi-phase flow. Comput. Phys, Commun, 87: 225-235, Monaghan, J.J, RAF. Cas, AM, Kos and M. Hallworth, 1999, Gravity currents descending a ramp in a stratified tank. J. Fluid Mech. 379: 39-69, Monaghan, J.J, 2008. Smoothed particle hydrodynamics. Rep. Prog, Phys., 68: 1703-1759. Raniles, P.W. and LD. Libersky, 1996, Smoothed Particle Hydodynamies-some recent improvements and applications. Comput. Methods Applied Mech. Eng., 139: 375-408, Ritchie, BW. and PA smoothed-particle hydrodynamics. Astron, Soe,, 323: 743-758, ‘Thomas, 2001. Multiphase Mon, Not. R. 3826