Jtam-A4.dvi JOURNAL OF THEORETICAL AND APPLIED MECHANICS 56, 3, pp. 675-686, Warsaw 2018 DOI: 10.15632/jtam-pl.56.3.675 SMOOTHED PARTICLE HYDRODYNAMICS MODELLING OF THE RAYLEIGH-PLATEAU INSTABILITY Michał Olejnik, Kamil Szewc Institute of Fluid-Flow Machinery, Polish Academy of Sciences, Gdańsk, Poland e-mail: michal.olejnik@imp.gda.pl; kamil.szewc@gmail.com The break-up of liquid ligaments and formation of droplets are elementary phenomena in multiphase flowswhich are of high importance in industrial andmedical applications. From the numerical point of view, they require proper interface and surface tension treatment. In the present work, we apply Smoothed Particle Hydrodynamics, a meshless approach, to simulate the break-up of a liquid cylinder inside the gaseous phase, i.e. theRayleigh-Plateau instability. Results obtained in 3D show that even a relatively coarse resolution allows one to predict correctly the size of droplets formed in the process. The detailed analysis of the break-up time in 2D setup implies that a certain level of spatial discretisation needs to be reached to determine this moment precisely. Keywords: meshless methods, SPH, capillary jet break-up, interfacial flows 1. Introduction Multiphase flows have been of interest for both academic and industrial communities for a long period of time. Accurate interface tracking and surface tension modelling are of particular importance due the influence on the flow solution as a whole. Another challenge has been the ability to properly treat high density ratios, frequently appearing in industrial problems. Over theyears,manydifferent approacheshavebeendeveloped to tackle thismatter, includingVolume of Fluid, Level Set andFrontTrackingmethods, for a comprehensive overview see the handbook by Tryggvason (2011). Among existingmethods for fluid flowmodelling, Smoothed Particle Hydrodynamics (SPH) remains a relatively newalternative. It is a particle basedmeshlessmethodofLagrangiannature. SPHwas originally developed in the 1970s for astrophysical problems, and later adapted for fluid mechanics simulations (Monaghan, 2012), where it gainedmore andmore interest over the recent years.Brancheswhere itsmeshless nature is notably favorable are free-surfaceflows (Violeau and Rogers, 2016) and general two-fluid flows with interfacial surfaces (Das and Das, 2010b; Szewc et al., 2013; Olejnik et al., 2016). SPH advantages in the latter case are mainly the easiness in dealing with high density ratios and straightforward treatment of the interface. Recently, applications of SPH have appeared in the area of microfluidics (Wieth et al., 2016). In the present paper, we focus on a particular application of SPH to interfacial flows which is theRayleigh-Plateau (R-P) instability. This fundamental phenomenon leads to atomisation of liquid jets and, likewise, the Kevin-Helmholtz instability (Boeck et al., 2007) is one of the basic mechanisms of the regime change in complexmultiphase flows. Another area where breaking-up liquid ligaments is of high concern is the formation of droplets of a desired size in drug delivery and lab-on-a-chip devices (Abate et al., 2009; Guzowski et al., 2013). In our work, we test feasibility of SPH in these problems. Although attempts in this matter have already beenmade for the case of gravity-driven dripping break-up, see Sirotkin and Yoh (2012), they focused on qualitative results and influence of dimensionless numbers. In the present study, we pay special attention to the resolution needed for accurate prediction of the liquid ligament break-up time. 676 M. Olejnik, K. Szewc 2. SPH for multiphase flows 2.1. Basics of the method The general idea behind SPH lies in interpolation theory. Let us consider any scalar (for simplicity) fieldA. The integral formula A(r)= ∫ Ω A(r′)δ(r−r′)dr′, (2.1) where δ(r) is the Dirac delta function, can be used to express the field value at the point r in space Ω. To obtain SPH approximation, we first replace δ(r) with the weighting kernel func- tionW(r,h) which should be normalised, symmetrical and converging to δ(r) with h→ 0 (Mo- naghan, 1992). Argumenth is the so-called smoothing length and it determines the interpolation range. In our work, we use the quintic kernel proposed byWendland (1995) W(r,h) =C    ( 1− q 2 )4 (2q+1) for q < 2 0 otherwise (2.2) where q = |r|/h and C is the normalisation constant (C = 7/4πh2 in 2D and C = 21/16πh3 in 3D), as guaranteeing good stability of computations (Dehnen and Aly, 2012; Szewc et al., 2012a). The second step consists in discretisation of the space into a set of particles of the volume Ωb = mb/̺b, where mb is mass and ̺b is density of the b-th particle. As a result, the integral from Eq. (2.1) is approximated by a sum, i.e. A(r)≃ ∑ b A(rb)W(r−rb,h)Ωb (2.3) In the shorthand notation, the SPH approximation 〈A〉a of the fieldA at any point a is defined as 〈A〉a = ∑ b AbWab(h)Ωb (2.4) where Ab = A(rb) and Wab(h) = W(ra− rb,h). Thanks to the properties of W(r,h), differen- tiation can be shifted from the field to the kernel function yielding 〈∇A〉a = ∑ b Ab∇Wab(h)Ωb (2.5) Further derivatives can be obtained in a similar way. Using theabovemethod, various kindsof differential equations canbe rewritten into theSPH formalism and solved by calculating interactions between particles, hence its wide application. More detailed information on the derivation of SPH and basics of the method can be found in the handbook by Violeau (2012). 2.2. Governing equations The set of governing equations for viscous flow consists of the Navier-Stokes (momentum) equation du dt =− 1 ̺ ∇p+ µ ̺ ∆u+a (2.6) Smoothed particle hydrodynamics modelling of the Rayleigh-Plateau instability 677 and the continuity equation d̺ dt =−̺∇·u (2.7) where u denotes the velocity vector, ̺ is fluid density, p is pressure and µ is dynamic viscosity. In Eq. (2.6), a stands for the interfacial term, detailed in Section 2.3. Due to the Lagrangian nature of SPH, we also include the advection equation dr dt =u (2.8) Depending on the purpose and assumptions, different SPH formalisms for the fluid flow can be obtained by using Eqs. (2.4) and (2.5). In the present study, we use a formulation proposed by Hu andAdams (2006). To the best of our knowledge, their approach is well suited formodelling multiphase flows with large density ratios (Szewc et al., 2012b). The pressure term in Eq. (2.6) will become 〈∇p ̺ 〉 a = 1 ma ∑ b (pa θ2a + pb θ2 b ) ∇aWab(h) (2.9) whereθ is the inverseof theparticle volume.Theviscous termofEq. (2.6), obtainedbycombining SPH formalism and Finite DifferenceMethod, takes the form 〈µ ̺ ∆u 〉 a = 1 ma ∑ b 2µaµb µa+µb ( 1 θ2a + 1 θ2 b )rab ·∇aWab(h) r2 ab +0.01h2 uab (2.10) where uab = ua −ub. The key feature of the approach proposed by Hu and Adams (2006) is the treatment of the continuity equation. Instead of rewriting Eq. (2.7) into the SPH language, density is found from 〈̺〉a =ma ∑ b Wab(h) =maθa (2.11) This allows thedensity field tobe representedonlyby the spatial distributionof particles andnot by their masses. Thanks to this, densities of particles near the interface are not affected by the other fluid, which is important in multiphase flow modelling. Note that Eq. (2.11) requires the wholedomain tobefilledwithSPHparticles, hence it is not suitable formultiphase computations inwhich the lighter phase is neglected, e.g. see theapproachproposedbyOrdoubadiet al. (2017). In thiswork,weuse theWeaklyCompressible SPHapproach (WCSPH).The set of equations is closed with the state equation p= s2̺0 γ [( ̺ ̺0 )γ −1 ] (2.12) where s is the artificial speed of sound, ̺0 is the reference (initial) density and γ is a numerical parameter. Values of c and γ are chosen to ensure density fluctuations at a level of 1% or below. Inmultiphase flowmodelling, it is a common practice to treat the liquid as incompressible with γ=7, and the gas as compressible with γ=1.4. We follow this approach in this study. 2.3. Surface tension The influence of surface tension is modelled with the Continuum Surface Force method (CSF), originally proposed by Brackbill et al. (1992), with SPH implementation due to Morris (2000). In this approach, surface tension forces are converted into a force per unit volume Fs = fsδs (2.13) 678 M. Olejnik, K. Szewc where δs is a suitably chosen surface delta function and fs =σκn̂ (2.14) is the force per unit area, σ is the surface tension coefficient, κ is the local curvature of the interface and n̂ is the unit vector normal to the interface. Using the so-called color function c (say, c=0 for the first phase and c=1 for the second one) n̂ can be calculated using the formula n̂= n |n| = ∇c |∇c| (2.15) The vector n is obtained from na = ̺a ∑ b (c̃b− c̃a)∇aWab(h)Ωb (2.16) where c̃ stands for the smoothed color function, i.e. c̃a = ∑ b cbWab(h)Ωb (2.17) The local curvature is obtained from κ=−∇· n̂ (2.18) Assuming that δs = |n|, the influence of surface tension can be included in Eq. (2.6) by adding the term aa = σ ̺a κana (2.19) Following Morris (2000), we also exclude from calculations of surface tension effects particles for which |na| is below the threshold of 0.01/h. This greatly improves accuracy of curvature estimation on fringes of the interfacial area. 2.4. Interface correction As already mentioned, SPH does not require any special treatment of the interface. In some cases, however, particles of two immiscible phases can penetrate into the bulk of the other phase. This is particularly visible in problems involving high density ratios. To prevent this unphysical phenomenon, we enforce a small repulsive interaction between different phases (Szewc et al., 2013), by adding to Eq. (2.6) the term Ξa = ε ma ∑ b cb 6=ca ( 1 θ2a + 1 θ2 b ) ∇aWab(h) (2.20) where ε is a numerical parameter. This correction, however, may change the characteristics of the flow if used improperly. Detailed analysis of this approach and the spurious interface fragmentation in the case of gas bubbles rising in the liquid can be found in the work of Szewc et al. (2015). Smoothed particle hydrodynamics modelling of the Rayleigh-Plateau instability 679 3. Results In order to validate the SPH approach for modelling of the break-up of liquid jets, we decided to perform simulations of the case presented by Dai and Schmidt (2005). We consider a fully periodic, cubic domain of the edge of length L containing a column of liquid l surrounded by the gaseous phase g. The column is placed in the centre of the domain according to ( y− L 2 )2 + ( z− L 2 )2 ¬ r20 (3.1) where r0 = 0.1L and x∈ [0,L]. The density and viscosity ratios are respectively ̺l/̺g = 1000 and µl/µg =100. The initial perturbation of the liquid velocity field is given as ux =u0 sin 2πx L uy =0 uz =0 (3.2) where u0 is the initial velocity amplitude. Dimensionless numbers describing this case are the Weber number We= ̺lr0u 2 0 σ (3.3) and the Reynolds number Re= ̺lu0r0 µl (3.4) In our study, we takeWe=1.4 and Re=18. The time is made dimensionless with tc = √ ̺l(r0) D σ (3.5) where D stands for the number of spatial dimensions. In the 2D case, we simply consider the central slice of the domain in the xy plane. In the present study, instead of rawparticles dataweanalyse data interpolated onto a regular grid. The reason for this is that the studied phenomena involve very thin liquid ligaments. Even for relatively fine resolution positions of SPH particles may be slightly perturbed, which creates impression of a misshapen interface. The grid used for post-processing has the cell size of ∆r, which is the initial spacing between the particles. The interpolation is performed with SPH formulation with the same weighting function, i.e. the Wendland kernel, that has been employed in calculations. This way, without distorting any information, detailed analysis is made significantly easier. Since in this work the shape of the interface is of utmost interest, we treat it as the isoline or isosurface of the color function c=0.5. To obtain the isoline/isosurface Python libraries Matplotlib (2D cases, (Hunter, 2007)) andMayavi (3D cases, (Ramachandran and Varoquaux, 2011)) have been used. Figure 1 shows an exemplary result of such treatment. 3.1. 3D case TheSPH simulations have been performed forL/h=64 andL/h=128withh/∆r=1.5625 which resulted in 106 particles filling the whole domain, while 31600 of them were forming the liquid cylinder for the lower resolution, and 8 · 106 to 252800 respectively for the higher one. Figure 2 presents results obtained from SPH simulation with L/h = 64 for different values of the interface sharpness correction term, Eq. (2.20), against the reference material (Dai and Schmidt, 2005). The general agreement is good, however, some discrepancies can be observed at later stages of simulations. The ligament break-up for SPH calculations occurs around t+ =5.5, 680 M. Olejnik, K. Szewc Fig. 1. Positions of SPH particles representing the liquid vs. isosurface of the color function c=0.5 interpolated onto a uniform grid while for the reference material it is t+ = 6.49. The main reason for this discrepancy is that in the SPH computations we do not use the Adaptive Mesh Refinement (AMR) or any similar techniques, contrary to calculations by Dai and Schmidt (2005). We should also mention that in the cited work, due to the use of AMR, the break-up time was defined as the moment when the liquid ligament between droplets reaches radius of 0.05r0. Note that AMR in grid- -based methods is a standard technique, while in SPH it is still a novelty in development, see Vacondio et al. (2016) or Olejnik et al. (2017). Furthermore, with an increase in the value of ε we observe a decrease in the diameter of the liquid bridge between the droplets at t+ = 4.49. This tendency is the outcome of additional interfacial pressure repelling phases from each other, see Eq. (2.20). The quantitative comparison is shown in Fig. 3. It shows the disturbance growth process quantified as rmax(t)−r0 r0 (3.6) where rmax(t) is themaximumdistance of the liquid phasemeasured from the axis of symmetry. The results show high agreement with the reference data, independently of the value of ε used. The higher resolution also does not yield significantly different results. We can conclude that even a coarse resolution is sufficient to correctly predict the size of formed droplets. As shown in Fig. 4, the resolution does influence their shape. For L/h=128, the curvature of the central droplet is lower than for L/h= 64. All in all, the quality of SPH simulation of the considered case proves that the method is suitable for simulations of the capillary jet break up provided that the resolution is high enough. The issue of the interface correction term and its influence on the flow requires further investigation. 3.2. 2D case Since SPH simulations in 3D are costly, for the purpose of a more detailed analysis we decided to settle on the 2D setup. To determine the influence of the resolution on the moment of break-up, the calculations have been performed for h/∆r = 2 and different values of the Smoothed particle hydrodynamics modelling of the Rayleigh-Plateau instability 681 Fig. 2. Results of SPH simulation of the R-P instability without the interface sharpness correction term (1st row), with ε=0.5 (2nd row) and ε=0.75 (3rd row). The 4th row presents the reference data of Dai and Schmidt (2005) reprinted with the permission fromElsevier 682 M. Olejnik, K. Szewc Fig. 3. Disturbance growth in time; SPH results compared to the referencematerial from Dai and Schmidt (2005) Fig. 4. Results of SPH simulation of the R-P instability without the interface sharpness correction term forL/h=64 (top row) andL/h=128 (middle row). In comparison to the reference data of Dai and Schmidt (2005) reprinted with the permission fromElsevier Smoothed particle hydrodynamics modelling of the Rayleigh-Plateau instability 683 smoothing length, i.e. L/h = 64, 128, 256, 384 and 512. Figure 5 shows the evolution of the interface in the simulation for L/h= 512. We defined the break-up time as the moment when the color function drops to the value c=0.5 in any of the points situated on the symmetry line of the liquid column. An example of the interface shape and the values of the color function at such a moment are presented in Fig. 6. As shown in Fig. 7, the outcome is highly dependent on the resolution. For the highest values of h the dependency is almost linear, i.e. the lower resolution, the earlier moment of break-up. Beginning from L/h = 256, the growth is barely visible, and it is safe to assume that this resolution is sufficient to fully resolve the flowwithout computational overhead. Analysis of the disturbance growth in time, defined with Eq. (3.6), confirms this statement. We see that the three highest resolutions tested are in agreement and hard to distinguishwhile the three lowest ones tend to diverge from them, especially in the later stage of simulation. Fig. 5. Evolution of 2DR-P instability. Result of SPH simulation forL/h=512 Fig. 6. Definition of the break-up time in SPH simulations; example forL/h=256. Shape of the interface (left) and color function profile (right) at t+ =1.684 684 M. Olejnik, K. Szewc Fig. 7. Dependency of the dimensionless moment of break-up on the resolution Fig. 8. Disturbance growth in time in the 2D case 4. Conclusions In the present study, we have successfully applied SPH to simulations of the R-P instability. Results obtained in the 3D case showed that SPH can predict droplets size with a comparable accuracy asameshbasedmethodusingAMR,despite the relatively low resolutionused.Analysis of the break-up time in 2D, however, showed that resolution needs to be on a relatively high level to correctly predict this moment. It is worth to note that thanks to the GPU parallelisation, for which SPH is exceptionally suitable (Szewc, 2014), it is still affordable for a desktop class computer. Themethod proposed in this paper can naturally be extended for other situations involving generations of microdroplets in the so-called lab-on-a-chip devices (Guzowski et al., 2013). The presence of solid boundaries in such devices is not an issue since reliable implementations of boundary conditions in SPH already exist (Adami et al., 2012). Recent works also show that thewetting phenomena and contact angles can be properly treatedwithin SPH framework (Das and Das, 2010a; Yeganehdoust et al., 2016). This makes SPH an interesting alternative to the traditional mesh based codes as a tool for engineering simulations. Optimisation of devices for e.g. precise and repeatable delivery of microdroplets sequences (Abate et al., 2009) can readily be performed with SPH. Acknowledgments K.S. was supported by Électricité de France (EDF) R&D project No. 8160-5910131187. M.O. was supported by the EUFP7 project Nugenia Plus No. 604965. The authors wish to thank Professor Jacek Pozorski for his encouragement and insightful discussions. Smoothed particle hydrodynamics modelling of the Rayleigh-Plateau instability 685 References 1. Abate A.R., Poitzsch A., Hwang Y., Lee J., Czerwinska J., Weitz D.A., 2009, Impact of inlet channel geometry onmicrofluidic drop formation,Physical Review E, 80, 026310 2. Adami S., Hu X.Y., Adams N.A., 2012, A generalized wall boundary condition for smoothed particle hydrodynamics, Journal of Computational Physics, 231, 7057-7075 3. BoeckT., Li J., López-Pagés E.,YeckoP., Zaleski S., 2007,Ligament formation in sheared liquid-gas layers,Theoretical and Computational Fluid Dynamics, 21, 59-76 4. Brackbill J.U., Kothe D.B., Zemach C., 1992, A continuum method for modeling surface tension, Journal of Computational Physics, 100, 335-354 5. DaiM., Schmidt S.P., 2005,Adaptive tetrahedralmeshing in free-surface flow, Journal of Com- putational Physics, 208, 228-252 6. Das A.K., Das P.K., 2010a, Equilibrium shape and contact angle of sessile drops of different volumes –Computation bySPHand its further improvementbyDI,Chemical Engineering Science, 65, 4027-4037 7. Das A.K., Das P.K., 2010b, Incorporation of diffuse interface in smoothed particle hydrodyna- mics: Implementation of the scheme and case studies, International Journal for NumericalMethods in Fluids, 67, 671-699 8. Dehnen W., Aly H., 2012, Improving convergence in smoothed particle hydrodynamics simula- tions without pairing instability,Monthly Notices of Royal Astronomical Society, 425, 1068-1082 9. Guzowski J., Jakiela S., Korczyk P.M., Garstecki P., 2013, Custom tailoring multiple droplets one-by-one,Lab on a Chip, 13, 4308-4311 10. Hu X.Y., Adams N.A., 2006, Amulti-phase SPHmethod for macroscopic andmesoscopic flows, Journal of Computational Physics, 213, 844-861 11. Hunter J.D., 2007,Matplotlib: A 2DGraphics Environment,Computing in Science and Engine- ering, 9, 90-95 12. Monaghan J.J., 1992, Smoothed Particle Hydrodynamics, Annual Review of Astronomy and Astrophysics, 30, 543-574 13. Monaghan J.J., 2012, Smoothed Particle Hydrodynamics and its diverse applications, Annual Review of Fluid Mechanics, 44, 323-346 14. Morris J.P., 2000, Simulating surface tension with smoothed particle hydrodynamics, Interna- tional Journal for Numerical Methods in Fluids, 33, 333-353 15. Olejnik M., Szewc K., Pozorski J., 2016, Modelling of the flow regime transition with the Smoothed Particle Hydrodynamics, Proceedings of 9th International Conference on Multiphase Flow, Florence, Italy, paper No. 1037 16. Olejnik M., Szewc K., Pozorski J., 2017, SPHwith dynamical smoothing length adjustment based on the local flow kinematics, Journal of Computational Physics, 348, 23-44 17. Ordoubadi M., Yaghoubi M., Yeganehdoust F., 2017, Surface tension simulation of free surfaceflowsusing smoothedparticle hydrodynamics,Scientia Iranica, Transactions B:Mechanical Engineering, 24, 2019-2033 18. RamachandranP.,VaroquauxG., 2011,Mayavi: 3Dvisualizationof scientificdata,Computing in Science and Engineering, 13, 40-51 19. Sirotkin F.V., Yoh J.J, 2012, A new particle method for simulating breakup of liquid jets, Journal of Computational Physics, 231, 1650-1674 20. Szewc K., 2014, Smoothed ParticleHydrodynamics simulations usingGraphics ProcessingUnits, TASK Quarterly, 18, 67-80 686 M. Olejnik, K. Szewc 21. SzewcK., Pozorski J.,Minier J.-P., 2012a,Analysis of the incompressibility constraint in the SPHmethod, International Journal for Numerical Methods in Engineering, 91, 343-369 22. SzewcK., Pozorski J.,Minier J.-P., 2013, Simulations of single bubbles rising throughviscous liquids using SmoothedParticleHydrodynamics, International Journal ofMultiphase Flow,50, 98- 105 23. Szewc K., Pozorski J., Minier J.-P., 2015, Spurious interface fragmentation in multiphase SPH, International Journal for Numerical Methods in Engineering, 103, 625-649 24. SzewcK.,TanièreA., Pozorski J.,Minier J.-P., 2012b,A study on application of Smoothed Particle Hydrodynamics to multi-phase flows, International Journal of Nonlinear Sciences and Numerical Simulation, 13(6), 383-395 25. Tryggvason G., Scardovelli R., Zaleski S., 2011, Direct Numerical Simulations of Gas- Liquid Multiphase Flows, Cambridge University Press 26. Wendland H., 1995, Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree,Advances in Computational Mathematics, 4, 389-396 27. Wieth L., Kelemen K., Braun S., Koch R., Bauer H.J., Schuchmann H.P., 2016, Smo- othed Particle Hydrodynamics (SPH) simulation of a high-pressure homogenization process, Mi- crofluidics and Nanofluidics, 20, 42 28. Vacondio R., Rogers B.D., Stansby P.K., Mignosa P., 2016, Variable resolution for SPH in three dimensions: Towards optimal splitting and coalescing for dynamic adaptivity, Computer Methods in Applied Mechanics and Engineering, 300, 442-460 29. Violeau D., 2012,Fluid Mechanics and the SPH Method, Oxford University Press, Oxford 30. VioleauD., Rogers B.D., 2016, Smoothed particle hydrodynamics (SPH) for free-surface flows: past, present and future, Journal of Hydraulic Research, 54, 1-26 31. Yeganehdoust F., YaghoubiM., EmdadH., Ordoubadi M., 2016,Numerical study ofmul- tiphase droplet dynamics and contact angles by smoothed particle hydrodynamics,AppliedMathe- matical Modelling, 40, 8493-8512 Manuscript received July 12, 2017; accepted for print October 27, 2017