An ISPH method for the simulation of three-phase °ows is presented in this

article. The proposed method is investigated through the simulation of a droplet located at the interface of two immiscible °uids as well as diamond droplet deformation. The extendibility of the proposed surface tension formulations for three-phase fows to two-phase fows is also investigated. It is observed that the results obtained from the numerical simulations are in good agreement with the analytical ones

Tofighi, Nima

Yıldız, Mehmet

Yildiz, Mehmet

Zainali, Amir

English

Sabanci University Research Database

II International Conference on Particle-based Methods - Fundamentals and ApplicationsPARTICLES 2011E. On˜ate and D.R.J. Owen (Eds)NUMERICAL INVESTIGATION OF THREE-PHASE FLOWSUSING INCOMPRESSIBLE SMOOTHED PARTICLEHYDRODYNAMICSAmir Zainali†, Nima Tofighi† and Mehmet Yildiz††Faculty of Engineering and Natural Sciences,Advanced Composites and Polymer Processing Laboratory,Sabanci University, Orhanli 34956, Tuzla, Istanbul, Turkeye-mail: meyildiz@sabanciuniv.edu,http://people.sabanciuniv.edu/meyildiz/Key words: SPH, multiphase flows, meshfree, phase specific surface tensionAbstract. An ISPH method for the simulation of three-phase flows is presented in thisarticle. The proposed method is investigated through the simulation of a droplet locatedat the interface of two immiscible fluids as well as diamond droplet deformation. Theextendibility of the proposed surface tension formulations for three-phase flows to two-phase flows is also investigated. It is observed that the results obtained from the numericalsimulations are in good agreement with the analytical ones.1 IntroductionTwo-phase flows are investigated extensively numerically in the literature during thepast decades while the physics behind the three and more phase flow problems are rela-tively unknown despite the fact that Multi-phase flows play an important role in manyindustrial and physical phenomena. Typical examples with many applications in chemi-cal and biomedical industries include mixture of three or more fluids (e.g. oil-water-gasmixture in oil industry) and compound droplets. The main reason for the scarcity ofnumerical simulations for three-phase flows in comparison to two-phase flows can be dueto the complications associated with the modeling of hydrodynamical interactions at theinterface for three-component systems [3, 4, 6].In this work, we have developed a 2-D multiphase incompressible smoothed particlehydrodynamics (ISPH) model based on the standard projection method initially used byCummins and Rudman [2] for SPH, while the continuum surface force model proposedby Brackbill et al. [1] is used to model surface tension forces. The method proposed bySmith et al. [6] is implemented to take into account different surface tension coefficientsbetween the phases in a three-component system.1A. Zainali, N. Tofighi and M. YildizThis paper is organized as follows. In section 2, the mathematical formulations as wellas numerical scheme are briefly described. In section 3, numerical results are presented.Finally in section 4 concluding remarks are provided.2 Governing EquationsAll constituents of the multiphase system are considered to be viscous, Newtonianand incompressible liquids with constant material properties, DΓ/Dt = 0, where D/Dt =∂/∂t+(),k vk is the material time derivative operator. The set of equations governing themotion of the flow under consideration are conservation of mass and linear momentum.In their local form for volume and discontinuity surface, these equations may be specifiedasDρ/Dt = −ρvk,k, on V-σ (1)ρDvl/Dt = τkl,k − ρbl, on V-σ (2)‖ρ (vk − uk)‖ nˆk = 0, on σ (3)‖ρvl (vk − uk)− τkl‖ nˆk = fσl , on σ (4)where V-σ denotes total volume, excluding the points lying on the surface of discontinuity,σ. Equations (1) and (2) are valid within the body excluding σ while equations (3) and(4) are valid only on the surface of discontinuity and represent the jump condition acrossσ. τkl = −pδkl + 2µdkl is the symmetric total stress tensor, dkl = 0.5 (vk,l + vl,k) isthe symmetric deformation tensor, δkl is the Kronecker delta, p is the thermodynamicpressure, bl is the body force, fσl is the surface force per unit area on the interface due tosurface tension for the mixture, vk is the divergence-free velocity, ρ and µ are the densityand viscosity of the mixture. The symbol ‖φ‖ indicates the jump of the enclosed quantityacross the discontinuity surface, φ+ − φ− where φ+ and φ− are the values of φ on thepositive and the negative sides of the discontinuity surface. Here, uk denotes the velocityof the discontinuity surface, nˆk is the unit normal to the discontinuity surface.One can write the local form of the jump condition for the momentum balance asfσl = ‖p‖nl = γκnˆl (5)where surface tension coefficient, γ, is assumed to be constant. Here, κ = −nm,m is thecurvature. For the sake of computational efficiency and simplification, it is preferable toexpress local surface force as an equivalent volumetric force (the force per unit volume).This can be achieved by using the continuum surface force model (CSF) proposed in [1].Using this approach equation (5) can be expressed as a volume force by multiplying thelocal surface tension force with a one-dimensional delta function defined on the interface,δ,fσvl = σκnˆlδ (6)2A. Zainali, N. Tofighi and M. YildizIn this study, to track the interface each particle is assigned to a color function,Cαi ={1, if i belongs to fluid α0, else(7)and finally, 1-D delta function on the interface is chosen to be δ = |C,m|.In the treatment of multiphase immiscible fluids with two constituents, the implemen-tation of the continuum surface force model is straightforward in that there is only oneinterface between the constituents. On the other hand, in the multiphase systems withmore than two constituents, an interface fluid particle may see two or more interfaces atthe same time and in what follows a decision should be made as to which surface tensioncoefficient is to be used. To be able resolve this difficulty, in [6] the surface tension coeffi-cient was decomposed into phase specific surface tension such that γαβ = γα + γβ whereγαβ is the physical surface tension coefficient between the phase α and the phase β, and γαand γβ are the phase specific surface tension coefficients for αth and βth phases respectively.For three-phase system, γ12 = γ1 + γ2, γ13 = γ1 + γ3, and γ23 = γ2 + γ3. Upon solving alinear system of equations for γ1, γ2 and γ3, one can write the following relations amongthe physical and phase specific surface tension coefficients as γ1 = 0.5 (γ12 + γ13 − γ23),γ2 = 0.5 (γ12 − γ13 + γ23), and γ3 = 0.5 (−γ12 + γ13 + γ23). As such, for each particle,three phase specific normals, curvatures, and the surface tension forces are computed. Thetotal surface tension force acting on the given interface particle is then the sum of threephase specific surface tension forces such that fσvil =∑3α=1 f(σv)αil where the subscript iis a particle identifier, and the adjacent Latin letter is the lth component of the surfaceforce.2.1 SPH preliminaries and the numerical methodThe numerical method used in this work for linearizing the governing equations andassociated interface and boundary conditions is based on the Incompresibble SmoothedParticle Hydrodynamics method presented in [5], which is briefly described in what fol-lows.SPH particles interact with each other by means of an interpolation (weighting, smooth-ing, kernel) function W (rij, h), concisely designated as Wij for a constant h. Here, rij isthe magnitude of the distance vector ~rij = ~ri−~rj between the particle of interest i and itsneighboring particles j, and h is referred to as the smoothing length which controls theinteraction length among particles.In the SPH method, an arbitrary function (i.e, scalar f(~ri), vectorial fp(~ri), or tensorialfps(~ri)) my be approximated as fpi =∫ΩfpjWijd3~rj where fp(~ri) is briefly denoted by fpi .Upon replacing the integral operation with a summation sign over all the particles withinthe cut-off distance (i.e. ~rij < Kh, where K is a constant coefficient for the selectedkernel function), and approximating the infinitesimal volume element by reciprocal of thenumber density ψj, defined for the particle j, one can write the SPH approximation of fpi3A. Zainali, N. Tofighi and M. Yildiz0 0.2 0.4 0.6 0.8 100.20.40.60.81xyFluid 1Fluid 2Fluid 3Figure 1: Domain size and fluid numbers at the initial condition for different simulations (The bubblehas a diameter of 0.3 with its center coordinates at 0.5, and 0.5.)asf pi = fp (ri) =∑j1ψjf pjWij (8)where the number density for particle i is defined as ψi =∑jWijUpon using a differentiable smoothing function, the spatial first order derivative of thefunction fpi evaluated at the location of particle i can simply be approximated throughmultiplying the gradient of the kernel function ∂Wij/∂xsi (taken with respect to the particlei) by the field variables and volumes of neighbor particles. There are several ways todiscretize the first order spatial derivative of a field variable. Throughout this work, theone used is of the following form∂f pi∂xkiaksij =∑j1ψj(fpj − f pi) ∂Wij∂xsi(9)where aksij is a corrective second-rank tensor given by,aksij =∑jrkijψj∂Wij∂xsi(10)This form is referred to as the corrective SPH gradient formulation that can be used toeliminate particle inconsistencies. It should be noted that the corrective term aksij is ideallyequal to Kronecker delta, δks for a continuous function. As for SPH approximation to4A. Zainali, N. Tofighi and M. Yildizthe second order derivative of a vector-valued function, it can be written in two differentways as∂2fpi∂xki ∂xkiapmij = 8∑j(fpi − f pj) rpijr2ij∂Wij∂xmi, (11)∂2fpi∂xki ∂xki(2 + assij)= 8∑j(fpi − fpj) rsijr2ij∂Wij∂xsi. (12)For the time integration in the ISPH approach, two-stage predictor-corrector methodwith a first-order Euler time step scheme is used. The algorithm starts with the predic-tor step where the intermediate positions ~r∗i for all particles are calculated through theknowledge of preceding particle positions ~r(n)i and the previous correct velocity field ~v(n)i as~r∗i = ~r(n)i +~v(n)i ∆t. The intermediate velocity field ~v∗i is calculated on the intermediate par-ticle locations by solving the momentum balance equations with forward time integrationwithout the pressure gradient term as ~v∗i = ~v(n)i +~f(n)i ∆t . The pressure Poisson equationis solved to obtain the pressure, p(n+1)i , which is required to enforce the incompressibilitycondition. The actual velocity field at time step (n+ 1), ~v(n+1)i ,can be obtained by usingthe computed pressure. Finally, with the correct velocity field for time-step (n + 1), allfluid particles are moved to their new positions ~r(n+1)i = ~r(n)i + 0.5(~v(n)i + ~v(n+1)i)∆t.3 Results and DiscussionIn this section, results of the conducted simulations are presented. Computationaldomain for every simulation is taken to be a square with the side length of 1. Eachrow and column consists of 100 equally spaced particles for test cases evolving from aCartesian particle arrangement. Test cases with non-Cartesian particle arrangement havethe same number of particles as with their Cartesian counterparts. All fluid properties areset to unity for every phase involved in simulations while binary surface tension coefficientassumes a value of 0 or 1, depending on the test cases considered. Different phases areordered as shown in figure 1. No-slip boundary conditions for the top and bottom, andside walls are applied. Zero pressure boundary condition is applied for all boundaries.In order to test the applicability of the three-phase surface tension treatment to two-phase flow systems, the deformation of the initially diamond-shaped droplet under theeffect of the surface tension force is considered. Particles of fluid 3 are initially positionedon a diamond surrounded by fluid 1 and 2 particles. The binary surface tension coefficientbetween fluids 1 and 2 is set to zero so as to make them act like a single fluid. Figure 2provides the sequences of shape-evolution of the diamond bubble in the aforementionedtest case in detail. The initial diamond evolves in time hence becoming a circle as thesimulation continues. The circle obtained in this way is also used for the initial positioningof particles in some other test cases. As expected, the new surface tension treatmentmethod presented for SPH of three-phase fluid flow is also applicable to a two-phase fluid5A. Zainali, N. Tofighi and M. Yildiz0.2 0.4 0.6 0.80.20.40.60.80.2 0.4 0.6 0.80.20.40.60.8x0.2 0.4 0.6 0.80.20.40.60.80.2 0.4 0.6 0.80.20.40.60.8t = 0y0.2 0.4 0.6 0.80.20.40.60.8t = 0.15y0.2 0.4 0.6 0.80.20.40.60.8t = 0.45yxFigure 2: The time evolution of a diamond shaped bubble to a circle shape bubble; Left: particle positions,Right: color contour showing bubble boundary.6A. Zainali, N. Tofighi and M. Yildizxy0.2 0.4 0.6 0.80.20.40.60.8x0.2 0.4 0.6 0.80.20.40.60.8(a) (b)Figure 3: Contours of surface tension force magnitude for different initial particle positioning at t =0.015 s; (a) crude circle marked in cartesian distribution, (b) circle obtained from evolution of a diamond.system.In this study, two different methods are used for the initial arrangement of particles in acircular bubble positioned in a rectangular domain for the liquid lens. In the first method,all the particles are placed on a Cartesian grid and the bubble is formed through assigningdifferent phases to each particle in accordance with their positions. The circular bubbleobtained by this method is rather crude thereby presenting many irregularities in shapewhich eventually result in an intermittent surface tension force as shown in figure 3-a. Inthe second method, an initial diamond shaped bubble is allowed to deform into a circularbubble in a two phase flow as described previously. The simulation using an initial circularbubble generated by the second method leads to improvements in the surface tension forceand particle positions, as seen in figure 3-b. The initial diameter for the crude bubbleused here is 0.293 while the bubble obtained from the diamond shape has a diameter of0.299. The reason behind this difference is that it is almost impossible to create a circularbubble shape of desired initial diameter (in this work, the initial diameter of interest is0.3) from limited number of equally spaced particles on a Cartesian grid. The crude circleexperiences an abrupt change in particle positions at early stages of the simulation due tothe aforementioned irregular surface tension force. This affects the initial reading for thediameter, rendering it difficult to compare the results with analytical data. This problemis more evident in figure 4-a, where lens diameter versus time is shown for both testcases. While both test cases show the same characteristics, the crude bubble has smallerelongation rate along the major axis. The results are more evident when the transversediameter of the lens is normalized by the analytically obtained equilibrium diameter, asshown in 4-b.7A. Zainali, N. Tofighi and M. Yildiz0 0.5 1 1.50.70.750.80.850.90.951D/Dat(s)0 0.5 1 1.50.250.30.350.40.45Dt(s)(a) (b)Figure 4: Change in transverse lens diameter versus time for different initial conditions; Solid line: crudecircle marked in Cartesian distribution, Dashed line: circle obtained from evolution of a diamond.0.2 0.4 0.6 0.80.20.40.60.8xyFigure 5: The time evolution of the liquid lens visualized by color contours.The shape of the droplet is controlled by the values of the surface tension forces. There-fore, under the effect of these surface tension forces, the initially circular droplet deformsinto an elliptic or lens shape. The equilibrium three-phase contact angle is formulated assinϕ1/γ23 = sinϕ2/γ13 = sinϕ3/γ12. The transverse diameter of the lens along the majoraxis (the distance between triple junctions) is computed through [4]d =(2 (pi − ϕ1)− sin (2pi − 2ϕ1)8A sin2 (pi − ϕ1)+2 (pi − ϕ3)− sin (2pi − 2ϕ3)8A sin2 (pi − ϕ3))(13)8A. Zainali, N. Tofighi and M. Yildizwhere ϕα is the contact angle of the αth phase, and A is the area of the liquid lens.0.2 0.4 0.6 0.80.20.40.60.8t = 0xy0.2 0.4 0.6 0.80.20.40.60.8t = 0.375xy0.2 0.4 0.6 0.80.20.40.60.8t = 0.75xy0.2 0.4 0.6 0.80.20.40.60.8t = 1.5xyFigure 6: Particle positions for liquid lens at different simulation times.Based on the initial diameter of circle obtained from a diamond, equation (13) predictsa diameter of 0.415 while the result obtained from the simulation is 0.413. It is obvious infigure 4-a that the difference between theory and simulation is much greater in the casewhere a crude circle is used as the initial condition for the simulations. The analyticalresult for the crude circle is 0.406 and the simulation results in an equilibrium diameterof 0.389. It is obvious in figure 4-b that the difference between theory and simulationis much greater in the case where a crude circle is used as the initial condition for the9A. Zainali, N. Tofighi and M. Yildizsimulations.Having chosen the second method for generating an initial circular bubble geometry,figure 5 shows the time evolution of a fluid lens between two layers of fluid. In this case,all fluid properties are set to 1 while binary surface tension coefficients, γij, are equalbetween each pair of surfaces. The time snapshots for particle positions inside the lens isalso provided in figure 6.4 ConclusionsIn this study, a three-phase ISPH method is developed and presented for the simulationof the three-phase flows. Phase specific surface tension coefficient [6] is used to obtainthe surface tension force for the systems with three components. The applicability ofthis formulations for two-phase flows is demonstrated. In addition, It is observed thatchoosing appropriate initial condition is crucial to obtain accurate results. And finally,it is shown that the equilibrium three-phase transverse lens diameter obtained is in goodagreement with the analytical one.5 AcknowledgementFunding provided by the European Commission Research Directorate General underMarie Curie International Reintegration Grant program with the Grant Agreement num-ber of PIRG03- GA-2008-231048 is gratefully acknowledged.References[1] J. U. Brackbill, D. B. Kothe, and C. Zemach. A continuum method for modelingsurface tension* 1. Journal of computational physics, 100(2):335–354, 1992.[2] S. J. Cummins and M. Rudman. An sph projection method. Journal of computationalphysics, 152(2):584–607, 1999.[3] J. Kim. Phase field computations for ternary fluid flows. Computer methods in appliedmechanics and engineering, 196(45-48):4779–4788, 2007.[4] J. S. Kim and J. Lowengrub. Phase field modeling and simulation of three-phase flows.Interfaces and free boundaries, 7(4):435, 2005.[5] M. S. Shadloo, A. Zainali, S. H. Sadek, and M. Yildiz. Improved incompressiblesmoothed particle hydrodynamics method for simulating flow around bluff bodies.Computer methods in applied mechanics and engineering, 200(9-12):1008–1020, 2011.[6] K. A. Smith, F. J. Souls, and D. L. Chopp. A projection method for motion of triplejunctions by level sets. Interfaces and free boundaries, 4(3):263–276, 2002.10

Numerical investigation of three-phase flows using incompressible smoothed particle hydrodynamics

An ISPH method for the simulation of three-phase ﬂows is presented in this article. The proposed method is investigated through the simulation of a droplet located at the interface of two immiscible ﬂuids as well as diamond droplet deformation. The extendibility of the proposed surface tension formulations for three-phase ﬂows to twophase ﬂows is also investigated. It is observed that the results obtained from the numerical simulations are in good agreement with the analytical ones

Toﬁghi, Nima

UPCommons. Portal del coneixement obert de la UPC

II International Conference on Particle-based Methods - Fundamentals and ApplicationsPARTICLES 2011E. Oñate and D.R.J. Owen (Eds)NUMERICAL INVESTIGATION OF THREE-PHASE FLOWSUSING INCOMPRESSIBLE SMOOTHED PARTICLEHYDRODYNAMICSAmir Zainali†, Nima Tofighi† and Mehmet Yildiz††Faculty of Engineering and Natural Sciences,Advanced Composites and Polymer Processing Laboratory,Sabanci University, Orhanli 34956, Tuzla, Istanbul, Turkeye-mail: meyildiz@sabanciuniv.edu,http://people.sabanciuniv.edu/meyildiz/Key words: SPH, multiphase flows, meshfree, phase specific surface tensionAbstract. An ISPH method for the simulation of three-phase flows is presented in thisarticle. The proposed method is investigated through the simulation of a droplet locatedat the interface of two immiscible fluids as well as diamond droplet deformation. Theextendibility of the proposed surface tension formulations for three-phase flows to two-phase flows is also investigated. It is observed that the results obtained from the numericalsimulations are in good agreement with the analytical ones.1 IntroductionTwo-phase flows are investigated extensively numerically in the literature during thepast decades while the physics behind the three and more phase flow problems are rela-tively unknown despite the fact that Multi-phase flows play an important role in manyindustrial and physical phenomena. Typical examples with many applications in chemi-cal and biomedical industries include mixture of three or more fluids (e.g. oil-water-gasmixture in oil industry) and compound droplets. The main reason for the scarcity ofnumerical simulations for three-phase flows in comparison to two-phase flows can be dueto the complications associated with the modeling of hydrodynamical interactions at theinterface for three-component systems [3, 4, 6].In this work, we have developed a 2-D multiphase incompressible smoothed particlehydrodynamics (ISPH) model based on the standard projection method initially used byCummins and Rudman [2] for SPH, while the continuum surface force model proposedby Brackbill et al. [1] is used to model surface tension forces. The method proposed bySmith et al. [6] is implemented to take into account different surface tension coefficientsbetween the phases in a three-component system.1843A. Zainali, N. Tofighi and M. YildizThis paper is organized as follows. In section 2, the mathematical formulations as wellas numerical scheme are briefly described. In section 3, numerical results are presented.Finally in section 4 concluding remarks are provided.2 Governing EquationsAll constituents of the multiphase system are considered to be viscous, Newtonianand incompressible liquids with constant material properties, DΓ/Dt = 0, where D/Dt =∂/∂t+ (),k vk is the material time derivative operator. The set of equations governing themotion of the flow under consideration are conservation of mass and linear momentum.In their local form for volume and discontinuity surface, these equations may be specifiedasDρ/Dt = −ρvk,k, on V-σ (1)ρDvl/Dt = τkl,k − ρbl, on V-σ (2)‖ρ (vk − uk)‖ n̂k = 0, on σ (3)‖ρvl (vk − uk) − τkl‖ n̂k = fσl , on σ (4)where V-σ denotes total volume, excluding the points lying on the surface of discontinuity,σ. Equations (1) and (2) are valid within the body excluding σ while equations (3) and(4) are valid only on the surface of discontinuity and represent the jump condition acrossσ. τkl = −pδkl + 2µdkl is the symmetric total stress tensor, dkl = 0.5 (vk,l + vl,k) isthe symmetric deformation tensor, δkl is the Kronecker delta, p is the thermodynamicpressure, bl is the body force, fσl is the surface force per unit area on the interface due tosurface tension for the mixture, vk is the divergence-free velocity, ρ and µ are the densityand viscosity of the mixture. The symbol ‖φ‖ indicates the jump of the enclosed quantityacross the discontinuity surface, φ+ − φ− where φ+ and φ− are the values of φ on thepositive and the negative sides of the discontinuity surface. Here, uk denotes the velocityof the discontinuity surface, n̂k is the unit normal to the discontinuity surface.One can write the local form of the jump condition for the momentum balance asfσl = ‖p‖nl = γκn̂l (5)where surface tension coefficient, γ, is assumed to be constant. Here, κ = −nm,m is thecurvature. For the sake of computational efficiency and simplification, it is preferable toexpress local surface force as an equivalent volumetric force (the force per unit volume).This can be achieved by using the continuum surface force model (CSF) proposed in [1].Using this approach equation (5) can be expressed as a volume force by multiplying thelocal surface tension force with a one-dimensional delta function defined on the interface,δ,fσvl = σκn̂lδ (6)2844A. Zainali, N. Tofighi and M. YildizIn this study, to track the interface each particle is assigned to a color function,Cαi ={1, if i belongs to fluid α0, else(7)and finally, 1-D delta function on the interface is chosen to be δ = |C,m|.In the treatment of multiphase immiscible fluids with two constituents, the implemen-tation of the continuum surface force model is straightforward in that there is only oneinterface between the constituents. On the other hand, in the multiphase systems withmore than two constituents, an interface fluid particle may see two or more interfaces atthe same time and in what follows a decision should be made as to which surface tensioncoefficient is to be used. To be able resolve this difficulty, in [6] the surface tension coeffi-cient was decomposed into phase specific surface tension such that γαβ = γα + γβ whereγαβ is the physical surface tension coefficient between the phase α and the phase β, and γαand γβ are the phase specific surface tension coefficients for αth and βth phases respectively.For three-phase system, γ12 = γ1 + γ2, γ13 = γ1 + γ3, and γ23 = γ2 + γ3. Upon solving alinear system of equations for γ1, γ2 and γ3, one can write the following relations amongthe physical and phase specific surface tension coefficients as γ1 = 0.5 (γ12 + γ13 − γ23),γ2 = 0.5 (γ12 − γ13 + γ23), and γ3 = 0.5 (−γ12 + γ13 + γ23). As such, for each particle,three phase specific normals, curvatures, and the surface tension forces are computed. Thetotal surface tension force acting on the given interface particle is then the sum of threephase specific surface tension forces such that fσvil =∑3α=1 f(σv)αil where the subscript iis a particle identifier, and the adjacent Latin letter is the lth component of the surfaceforce.2.1 SPH preliminaries and the numerical methodThe numerical method used in this work for linearizing the governing equations andassociated interface and boundary conditions is based on the Incompresibble SmoothedParticle Hydrodynamics method presented in [5], which is briefly described in what fol-lows.SPH particles interact with each other by means of an interpolation (weighting, smooth-ing, kernel) function W (rij, h), concisely designated as Wij for a constant h. Here, rij isthe magnitude of the distance vector rij = ri −rj between the particle of interest i and itsneighboring particles j, and h is referred to as the smoothing length which controls theinteraction length among particles.In the SPH method, an arbitrary function (i.e, scalar f(ri), vectorial fp(ri), or tensorialfps(ri)) my be approximated as fpi =∫Ωfpj Wijd3rj where fp(ri) is briefly denoted by fpi .Upon replacing the integral operation with a summation sign over all the particles withinthe cut-off distance (i.e. rij < Kh, where K is a constant coefficient for the selectedkernel function), and approximating the infinitesimal volume element by reciprocal of thenumber density ψj, defined for the particle j, one can write the SPH approximation of fpi3845A. Zainali, N. Tofighi and M. Yildiz0 0.2 0.4 0.6 0.8 100.20.40.60.81xyFluid 1Fluid 2Fluid 3Figure 1: Domain size and fluid numbers at the initial condition for different simulations (The bubblehas a diameter of 0.3 with its center coordinates at 0.5, and 0.5.)asfpi = fp (ri) =∑j1ψjfpj Wij (8)where the number density for particle i is defined as ψi =∑j WijUpon using a differentiable smoothing function, the spatial first order derivative of thefunction fpi evaluated at the location of particle i can simply be approximated throughmultiplying the gradient of the kernel function ∂Wij/∂xsi (taken with respect to the particlei) by the field variables and volumes of neighbor particles. There are several ways todiscretize the first order spatial derivative of a field variable. Throughout this work, theone used is of the following form∂fpi∂xkiaksij =∑j1ψj(fpj − fpi) ∂Wij∂xsi(9)where aksij is a corrective second-rank tensor given by,aksij =∑jrkijψj∂Wij∂xsi(10)This form is referred to as the corrective SPH gradient formulation that can be used toeliminate particle inconsistencies. It should be noted that the corrective term aksij is ideallyequal to Kronecker delta, δks for a continuous function. As for SPH approximation to4846A. Zainali, N. Tofighi and M. Yildizthe second order derivative of a vector-valued function, it can be written in two differentways as∂2fpi∂xki ∂xkiapmij = 8∑j(fpi − fpj) rpijr2ij∂Wij∂xmi, (11)∂2fpi∂xki ∂xki(2 + assij)= 8∑j(fpi − fpj) rsijr2ij∂Wij∂xsi. (12)For the time integration in the ISPH approach, two-stage predictor-corrector methodwith a first-order Euler time step scheme is used. The algorithm starts with the predic-tor step where the intermediate positions r∗i for all particles are calculated through theknowledge of preceding particle positions r(n)i and the previous correct velocity field v(n)i asr∗i = r(n)i +v(n)i ∆t. The intermediate velocity field v∗i is calculated on the intermediate par-ticle locations by solving the momentum balance equations with forward time integrationwithout the pressure gradient term as v∗i = v(n)i +f(n)i ∆t . The pressure Poisson equationis solved to obtain the pressure, p(n+1)i , which is required to enforce the incompressibilitycondition. The actual velocity field at time step (n + 1), v(n+1)i ,can be obtained by usingthe computed pressure. Finally, with the correct velocity field for time-step (n + 1), allfluid particles are moved to their new positions r(n+1)i = r(n)i + 0.5(v(n)i + v(n+1)i)∆t.3 Results and DiscussionIn this section, results of the conducted simulations are presented. Computationaldomain for every simulation is taken to be a square with the side length of 1. Eachrow and column consists of 100 equally spaced particles for test cases evolving from aCartesian particle arrangement. Test cases with non-Cartesian particle arrangement havethe same number of particles as with their Cartesian counterparts. All fluid properties areset to unity for every phase involved in simulations while binary surface tension coefficientassumes a value of 0 or 1, depending on the test cases considered. Different phases areordered as shown in figure 1. No-slip boundary conditions for the top and bottom, andside walls are applied. Zero pressure boundary condition is applied for all boundaries.In order to test the applicability of the three-phase surface tension treatment to two-phase flow systems, the deformation of the initially diamond-shaped droplet under theeffect of the surface tension force is considered. Particles of fluid 3 are initially positionedon a diamond surrounded by fluid 1 and 2 particles. The binary surface tension coefficientbetween fluids 1 and 2 is set to zero so as to make them act like a single fluid. Figure 2provides the sequences of shape-evolution of the diamond bubble in the aforementionedtest case in detail. The initial diamond evolves in time hence becoming a circle as thesimulation continues. The circle obtained in this way is also used for the initial positioningof particles in some other test cases. As expected, the new surface tension treatmentmethod presented for SPH of three-phase fluid flow is also applicable to a two-phase fluid5847A. Zainali, N. Tofighi and M. Yildiz0.2 0.4 0.6 0.80.20.40.60.80.2 0.4 0.6 0.80.20.40.60.8x0.2 0.4 0.6 0.80.20.40.60.80.2 0.4 0.6 0.80.20.40.60.8t = 0y0.2 0.4 0.6 0.80.20.40.60.8t = 0.15y0.2 0.4 0.6 0.80.20.40.60.8t = 0.45yxFigure 2: The time evolution of a diamond shaped bubble to a circle shape bubble; Left: particle positions,Right: color contour showing bubble boundary.6848A. Zainali, N. Tofighi and M. Yildizxy0.2 0.4 0.6 0.80.20.40.60.8x0.2 0.4 0.6 0.80.20.40.60.8(a) (b)Figure 3: Contours of surface tension force magnitude for different initial particle positioning at t =0.015 s; (a) crude circle marked in cartesian distribution, (b) circle obtained from evolution of a diamond.system.In this study, two different methods are used for the initial arrangement of particles in acircular bubble positioned in a rectangular domain for the liquid lens. In the first method,all the particles are placed on a Cartesian grid and the bubble is formed through assigningdifferent phases to each particle in accordance with their positions. The circular bubbleobtained by this method is rather crude thereby presenting many irregularities in shapewhich eventually result in an intermittent surface tension force as shown in figure 3-a. Inthe second method, an initial diamond shaped bubble is allowed to deform into a circularbubble in a two phase flow as described previously. The simulation using an initial circularbubble generated by the second method leads to improvements in the surface tension forceand particle positions, as seen in figure 3-b. The initial diameter for the crude bubbleused here is 0.293 while the bubble obtained from the diamond shape has a diameter of0.299. The reason behind this difference is that it is almost impossible to create a circularbubble shape of desired initial diameter (in this work, the initial diameter of interest is0.3) from limited number of equally spaced particles on a Cartesian grid. The crude circleexperiences an abrupt change in particle positions at early stages of the simulation due tothe aforementioned irregular surface tension force. This affects the initial reading for thediameter, rendering it difficult to compare the results with analytical data. This problemis more evident in figure 4-a, where lens diameter versus time is shown for both testcases. While both test cases show the same characteristics, the crude bubble has smallerelongation rate along the major axis. The results are more evident when the transversediameter of the lens is normalized by the analytically obtained equilibrium diameter, asshown in 4-b.7849A. Zainali, N. Tofighi and M. Yildiz0 0.5 1 1.50.70.750.80.850.90.951D/Dat(s)0 0.5 1 1.50.250.30.350.40.45Dt(s)(a) (b)Figure 4: Change in transverse lens diameter versus time for different initial conditions; Solid line: crudecircle marked in Cartesian distribution, Dashed line: circle obtained from evolution of a diamond.0.2 0.4 0.6 0.80.20.40.60.8xyFigure 5: The time evolution of the liquid lens visualized by color contours.The shape of the droplet is controlled by the values of the surface tension forces. There-fore, under the effect of these surface tension forces, the initially circular droplet deformsinto an elliptic or lens shape. The equilibrium three-phase contact angle is formulated assin ϕ1/γ23 = sin ϕ2/γ13 = sin ϕ3/γ12. The transverse diameter of the lens along the majoraxis (the distance between triple junctions) is computed through [4]d =(2 (π − ϕ1) − sin (2π − 2ϕ1)8A sin2 (π − ϕ1)+2 (π − ϕ3) − sin (2π − 2ϕ3)8A sin2 (π − ϕ3))(13)8850A. Zainali, N. Tofighi and M. Yildizwhere ϕα is the contact angle of the αth phase, and A is the area of the liquid lens.0.2 0.4 0.6 0.80.20.40.60.8t = 0xy0.2 0.4 0.6 0.80.20.40.60.8t = 0.375xy0.2 0.4 0.6 0.80.20.40.60.8t = 0.75xy0.2 0.4 0.6 0.80.20.40.60.8t = 1.5xyFigure 6: Particle positions for liquid lens at different simulation times.Based on the initial diameter of circle obtained from a diamond, equation (13) predictsa diameter of 0.415 while the result obtained from the simulation is 0.413. It is obvious infigure 4-a that the difference between theory and simulation is much greater in the casewhere a crude circle is used as the initial condition for the simulations. The analyticalresult for the crude circle is 0.406 and the simulation results in an equilibrium diameterof 0.389. It is obvious in figure 4-b that the difference between theory and simulationis much greater in the case where a crude circle is used as the initial condition for the9851A. Zainali, N. Tofighi and M. Yildizsimulations.Having chosen the second method for generating an initial circular bubble geometry,figure 5 shows the time evolution of a fluid lens between two layers of fluid. In this case,all fluid properties are set to 1 while binary surface tension coefficients, γij, are equalbetween each pair of surfaces. The time snapshots for particle positions inside the lens isalso provided in figure 6.4 ConclusionsIn this study, a three-phase ISPH method is developed and presented for the simulationof the three-phase flows. Phase specific surface tension coefficient [6] is used to obtainthe surface tension force for the systems with three components. The applicability ofthis formulations for two-phase flows is demonstrated. In addition, It is observed thatchoosing appropriate initial condition is crucial to obtain accurate results. And finally,it is shown that the equilibrium three-phase transverse lens diameter obtained is in goodagreement with the analytical one.5 AcknowledgementFunding provided by the European Commission Research Directorate General underMarie Curie International Reintegration Grant program with the Grant Agreement num-ber of PIRG03- GA-2008-231048 is gratefully acknowledged.References[1] J. U. Brackbill, D. B. Kothe, and C. Zemach. A continuum method for modelingsurface tension* 1. Journal of computational physics, 100(2):335–354, 1992.[2] S. J. Cummins and M. Rudman. An sph projection method. Journal of computationalphysics, 152(2):584–607, 1999.[3] J. Kim. Phase field computations for ternary fluid flows. Computer methods in appliedmechanics and engineering, 196(45-48):4779–4788, 2007.[4] J. S. Kim and J. Lowengrub. Phase field modeling and simulation of three-phase flows.Interfaces and free boundaries, 7(4):435, 2005.[5] M. S. Shadloo, A. Zainali, S. H. Sadek, and M. Yildiz. Improved incompressiblesmoothed particle hydrodynamics method for simulating flow around bluff bodies.Computer methods in applied mechanics and engineering, 200(9-12):1008–1020, 2011.[6] K. A. Smith, F. J. Souls, and D. L. Chopp. A projection method for motion of triplejunctions by level sets. Interfaces and free boundaries, 4(3):263–276, 2002.10852

https://upcommons.upc.edu/bitstream/2117/189892/1/Particles_2011-79_Numerical%20investigation%20of%20three-phase.pdf

Numerical investigation of three-phase flows using incompressible smoothed particle hydrodynamics

Abstract

Similar works

Full text

Available Versions

Sabanci University Research Database

UPCommons. Portal del coneixement obert de la UPC