A multi-phase SPH model based on Riemann solvers for simulation of jet breakup

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Contents lists available at ScienceDirect

Engineering Analysis with Boundary Elements journal homepage: www.elsevier.com/locate/enganabound

A multi-phase SPH model based on Riemann solvers for simulation of jet breakup Qiuzu Yang a,b, Fei Xu a,b,∗, Yang Yang a,b, Lu Wang a,b a b

School of Aeronautics, Northwestern Polytechnical University, 710072 Xi’an, PR China Institute for Computational Mechanics and Its Application, Northwestern Polytechnical University, 710072 Xi’an, PR China

a r t i c l e Keywords: Multiphase flow Discontinuous interface SPH Steady jet Pulsating jet Jet breakup length

i n f o

a b s t r a c t The discontinuous physical properties of the multiphase flow interface lead to the numerical instability, especially for the multiphase flow problems with large density or viscosity ratios. Based on Riemann solution, a novel multiphase SPH model is presented. The momentum equation of the Riemann form is improved to calculate the physical viscosity of multiphase flow and decrease the Riemann dissipation. The solid-wall boundary is imposed by one-sided Riemann problem. In addition, a repulsive force and the surface tension are applied on both sides of the multiphase interface to keep the interface sharpness and consider the small-scale interface effect. Three cases, squared droplet oscillating, single bubble rising (with different density and viscosity ratios) and two-bubble rising examples, are simulated to demonstrate the accuracy and effectiveness of the proposed method in dealing with the multiphase flow problems with the wide range of density and viscosity ratios and with the complex interfaces. Finally, the jet breakup problems of steady jets and pulsating jets are studied using the present model. A detailed process of jet breakup is observed, the surface wave structure of liquid jet could be identified. The influence of jet parameters (surface tension, velocity amplitude and velocity period) on physical information of jet is analyzed.

1. Introduction Jet breakup problems have been widely applied in many fields, such as pesticide spray, inkjet printing, nuclear reactor, fuel injection and so on. In the field of internal combustion engine [1], the manner and form of the fuel jet breakup directly affect the combustibility of the liquid fuels for internal combustion engine and the power performance of the engine. And, in the field of spray forming [2], the processing technology parameters such as the exit velocity and pulsation period of nozzle injection have a significant impact on the production efficiency and the product property of the material. Therefore, the study of jet breakup is of great theoretical value and practical significance. Early the characteristics of jet breakup had been studied through theory and experiment. From the literature [3], Savart (1833) was the first to discover the phenomenon of formation of planar liquid film through impinging jet experiments. Rayleigh [4] proposed the first jet theoretically model in 1878, and Tomotika [5] (1935) extended the Rayleigh’s model. Later, Meeister [6] (1969) presented an improved theory to study the jet breakup and droplet formation. Lin and Woods [7] (1991) also performed an experiment and showed that the excitation of jet could accelerate the jet frature at linear instability frequency. In recent years, due to the rapid development of computer technology, the computa∗

tional fluid dynamics has been widely used in fluid dynamics, which resulted in the development of numerical simulation of jet breakup problems. Richards et al. [8] used Volume-of-Fluid method (VOF) to solve the jet breakup problems under the axisymmetric Rayleigh model, and discussed the influence of jet parameters on the breakup length and droplet size of the jet. And, Pan and Suga [9] applied Level-Set method (LS) to simulate laminar jet breakup, and studied the development of the jet surface wave and the jet breakup mechanisms. Saito et al. [10] also studied the jet breakup problems by LBM (Lattice Boltzmann Modeling) method. Ménard et al. [11] combined LS, VOF and GFM (Ghost Fluid Method) methods to show the dynamic process of three-dimensional jet entering static air. To reduce the mesh size and save computing resources, Fuster et al. [12] used AMR (Adaptive Mesh Refinement) technology to simulate the primary jet breakup process. Nevertheless, the above methods are based on Eulerian grid-base numerical methods, they need specific interface tracking algorithms and grid reconstruction process to obtain valid numerical results, so the accuracy of capturing interfaces is low. In order to solve this problem, high resolution grids is needed, which increases the calculation cost. Compared with the methods above, particle methods are an ideal solution for the jet breakup problems with liquid breakup. Such as Smoothed Particle Hydrodynamics (SPH) [13], Moving Particle Semi-implicit method (MPS) [14] and

Corresponding author at: School of Aeronautics, Northwestern Polytechnical University, 710072 Xi’an, PR China. E-mail address: [email protected] (F. Xu).

https://doi.org/10.1016/j.enganabound.2019.10.015 Received 12 May 2019; Received in revised form 26 July 2019; Accepted 29 October 2019 0955-7997/© 2019 Elsevier Ltd. All rights reserved.

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

so on. SPH is initially introduced by Lucy and Gingold & Monaghan on dealing with astrophysics problems. Because of its programming simplicity and phase-interface tracing nature without introducing interface tracking algorithms, SPH is widely used in computational fluid dynamics in recent years [15–17]. However, there are only a few researches about the jet breakup problems using particle methods. Shibata et al. [18] (2004) applied MPS method to simulate the jet breakup of inviscid liquid released from a nozzle, and analyzed the effects of the Weber number and the Froude number on the size distribution of jets. Ganzenmüller et al. [19] (2007) used the SPH method to model the diesel injection by considering the influence of liquid surface tension. Later, Sirotkin and Yoh [20] (2012) applied a corrected SPH method to study the viscous jet breakup. They reproduced the behaviors of the jet breakup from the liquid’s jets to dripping faucets by reducing the Weber number to below the critical value. Hoefler et al. [21] (2013) also adopted the SPH method to simulate the evolution process of low and high viscosity spray drived by gravity. They observed the Rayleigh breakup by decreasing the spray’s Reynold number. Pourabdian et al. [22] (2016) carried out the simulation of breakup of liquid jet using SPHysics open source code, and investigated the relationship between length of liquid breakup in Rayleigh regime and Reynolds and Weber numbers, which is agreement with experimental results. The above numerical simulations using particle methods neglect the existence of air around the liquid jet. However, in practical engineering applications, there are other liquids around the jet, and the surface tension between the jet and the liquid cannot be ignored. Although the surface tension is exerted on the jet surface by searching free surface particle of single-phase SPH [20], the accuracy is low due to the particle lack in the computational domain of free surface particles. Therefore, it is necessary to find a multi-phase SPH model (MSPH) for simulating the jet breakup problems with the large density and viscosity ratios. As to MSPH models, there are great challenges to deal with the physical field discontinuity (such as density, viscosity and so on) across the interface, especially the multiphase flow problems with high density ratio and viscosity ratio. The large density and viscosity drop contribute to the acceleration discontinuity across the interface, which is liable to cause spurious oscillations in the pressure and velocity field. To address this issue, many improvements have been made by researchers. In 2003, Colagrossi and Landrini [23] adopted density reinitialization at intermediate time steps to maintain the consistency of particle mass, density and volume and reduce the pressure oscillation, and the artificial viscosity and XSPH technique were used in multiphase SPH model. In 2009, Grenier et al. [24] proposed a MSPH model, in which a renormalization of the kernel function was introduced to recover the interpolation accuracy of free surface flows, but it would increase additional computational cost. In 2010, Adami et al. [25] used a density-weight average pressure to solve the discontinuity of pressure gradient across the multiphase flow interface, but this did not guarantee interface stability in calculating the large density ratio multiphase flows. In 2013, Monaghan and Rafiee [26] proposed a simple MSPH model which can simulate large density ratio multiphase flow. However, in order to stabilize the interfaces, it was necessary to artificial treat the low-density interface particles as a rigid body during the simulation of large deformation interface problems (i.e. Rayleigh-Taylor instability). In 2015, Chen et al. [27] pointed out the problem existed in the literature [25] and presented a new MSPH model. This model assumed the pressure continuity across the interface and regarded other phases’ information as the same phase except for density and mass, and the artificial viscosity was considered to alleviate pressure oscillations. In 2016, Lind et al. [28] also proposed a new multiphase SPH algorithm to simulate large density ratio multiphase flow by coupling weak compressible SPH and incompressible SPH, in which the artificial viscous and the particle shifting algorithm were used to reduce numerical oscillation of the compressible phase and particle aggregation of incompressible phase. In 2017, Rezavand et al. [29] proposed an incompressible SPH multiphase flow model which did not involve velocity divergence term, the density and viscosity were smoothed using

an arithmetic mean interpolation. Recently, Zheng et al. [30] presented a low-diffusion multiphase SPH model, a new switch-function-based artificial viscosity term is used to reduce numerical diffusion and maintain stability. All above MSPH models focus on improving the pressure oscillation and stabilizing the interface, the non-physical artificial technology such as artificial viscosity, XSPH and so on are adopted in most of them, and multi variables are introduced by the artificial technology. But the value of variables is not easy to determine, which would introduce excessive numerical dissipation. By reconstructing SPH with convolution integral, Inutsuka [31] proposed a new SPH algorithm with no artificial viscosity based on Riemann solution (which is called GSPH). Later, Monaghan [32] combined Riemann solution with SPH to simulate the shock tube problem, which not only effectively solved the numerical oscillation of the discontinuity, but also improved the calculation accuracy. Therefore, Parshikov [33] used the Riemann solution of sound wave approximation to describe the interaction between particles, which solved the problem of strong shock oscillation and broaden the scope of SPH method at the same time [34–35]. Later, Cai et al. [36] presented a multiphase Godunov-type smoothed-particle hydrodynamics to simulate the multiphase Riemann problems, including air-helium shok tube, stiff airhelium shock tube and so on. Puri [37] established an equivalence between the disspative terms in GSPH and the signal based SPH artificial viscosity and gave the advice to the modification of the dissipation term in GSPH. Recently, based on a low-dissipation Riemann solver (which is obtained by modifying the Riemann dissipative term) and the transportvelocity formulation, Rezavand [38] proposed a SPH method to simulate violent multi-phase problems. Taking full advantage of Riemann solution in dealing with contact discontinuity problems, the present work proposes a multiphase SPH model based on Riemann solvers (R-MSPH) to study the jet breakup problems with complex interfaces and large density ratios. To decrease the Riemann dissipation and calculate the physical viscosity of multiphase flow, the SPH momentum equation of the Riemann form is improved, one-sided Riemann problem is used to impose multiphase flow solid-wall boundary, and interface sharpness is maintained using surface tension and a repulsive force on both sides of the interface. No artificial dissipation technology is added to the model. The layout of the paper is as follows: Section 2 outlines the governing equations of multiphase flow multiphase, and describes the presented multiphase SPH model based on Riemann solvers in detail. A number of multiphase test cases are presented in Section 3. In Section 4, the present model is applied to simulation of the jet breakup problems. A sufficient validation is presented for the steady and pulsating jets. The influencing factors of the steady and pulsating jet breakup length, such as surface tension, velocity amplitude and period, are investigated. The concluding remarks are reported in Section 5.

2. Numerical models 2.1. Governing equation of the multiphase flow The governing equations of multi-phase flow considered here can be described by the Navier–Stokes equations for the conservation of mass and momentum. Expressed in Lagrangian form the governing equations are as follow: ⎧ d𝜌 = −𝜌∇ ⋅ 𝐯 ⎪ d𝑡 ⎨ d𝐯 1 ⎪ = − ∇𝑝 + 𝐠 ⎩ d𝑡 𝜌

(1)

Where, 𝜌 is density, t is the time, p is pressure, and v is velocity vector, d 𝜕 = 𝜕𝑡 + 𝐯 ⋅ ∇ is the material derivative, g is the gravitational accelerad𝑡 tion. 135

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

To solve the Eq. (1), the following equation of state is applied in this paper [23]: ] [( )𝛾 𝜌 𝑐2 𝜌 𝑝= 0 − 1 + 𝑝0 (2) 𝛾 𝜌0 Where, 𝜌0 is the initial density, 𝛾 is the polytropic constant, and c is the numerical sound of speed. 𝛾 and c are used to keep the variation of particle density within 1% and to fulfill the incompressible limit. When the density ratio meets this condition 𝜌a /𝜌b ≥10, the heavier phase 𝛾=7 is used and the lighter phase 𝛾 = 1.4 is used in all cases. For the multiphase flow cases with smaller density ratios, the parameter 𝛾 of the lighter phase can be set to 3.5 or 7, which is depend on the complexity of the multiphase flow problems (such as the viscosity ratios). To ensure the initial pressure is smooth at the interface, the different phase have the same coefficient 𝜌0 c2 /𝛾. p0 is the background pressure which can avoid negative particle pressure and guarantee the computational stability [39]. When considering inertia force, gravity, and surface tension, the selection of the numerical sound speed can refer to the literature [39].

Fig. 1. Wave solution of Riemann solver.

Where, q=|ri -rj |/h, one-dimensional problems: 𝛼=3/64 h, twodimensional problems: 𝛼=7/64h2 , three-dimensional problems: 𝛼=21/256h3 .

2.2. The conventional multiphase SPH model 2.3. The R-MSPH model The basic idea of the SPH method is to discretize a continuous material field by a set of moving particles. Each particle has its own physical characteristics (such as velocity, pressure, energy, and mass), and represents an interpolation point of the known physical field. The particle’s mass remains invariant throughout the numerical solution, while other characteristics will be updated by interpolating the physical quantities over the neighboring particles using a weighting function W at each time step. By calculating the property values of all the particles, the solution of the entire material field problem is obtained. Based on this method, a function and a derivative representing the physical quantity of the particle i are expressed as the following forms [40]: ∑ 𝑓 ( 𝒙𝑖 ) = 𝑓 (𝒙′ )𝑊 (𝐱𝑖 − 𝐱′ , ℎ) 𝑑 𝐱′ ≈ 𝑓 (𝐱𝑗 )𝑊 (𝐱𝑖 − 𝐱𝑗 , ℎ)𝑉𝑗 (3) ∫Ω 𝑗 ∇𝑓 (𝐱𝑖 ) =

∫Ω

𝑓 (𝐱′ )∇𝑊 (𝐱𝑖 − 𝐱′ , ℎ) 𝑑 𝐱′ ≈

) − 𝑓 (𝐱𝑖 ) ∇𝑊 (𝐱𝑖 − 𝐱𝑗 , ℎ)𝑉𝑗

The traditional SPH method describes the interaction between particles by means of the average densities, velocities and stresses of the surrounding particles. Now we assume that the interaction exists in the middle of two particles. The Riemann problem between each pair of particles is used to describe the interaction between the central particle i and the neighboring particle j and the interface in normal direction is along a unit vector eij = -rij /|rij | [37,38]. The physical quantities of the two related particles constitute the initial left and right states of the Riemann problem. The left and right states are: ) ( ) {( 𝜌𝐿 , 𝑢𝐿 , 𝑃𝐿 = 𝜌𝑖 , 𝐯𝑖 ⋅ 𝐞𝑖𝑗 , 𝑃𝑖 (7) ( ) ( ) 𝜌𝑅 , 𝑢𝑅 , 𝑃𝑅 = 𝜌𝑗 , 𝐯𝑗 ⋅ 𝐞𝑖𝑗 , 𝑃𝑗 The first-order Riemann solvers using four states divided by three waves is used to approximate the Riemann problem. As shown in Fig. 1. Assuming that the velocity and pressure of left intermediate state are equal to the right, the linear Riemann solvers based on initial variables can be given as [42]: 𝐯𝑖 + 𝐯𝑗 1 𝑝𝐿 − 𝑝𝑅 ⎧ ∗ ∗ ∗ + 𝐞 ⎪𝑢𝐿 = 𝑢𝑅 = 𝐯𝑖𝑗 = 2 2 𝜌̄𝑖𝑗 𝑐̄𝑖𝑗 𝑖𝑗 (8) ⎨ 𝑝𝑖 + 𝑝𝑗 1 ⎪𝑃 ∗ = 𝑃 ∗ = 𝑃 ∗ = + 𝜌 ̄ 𝑐 ̄ ( 𝑢 − 𝑢 ) 𝑖𝑗 𝑖𝑗 𝐿 𝑅 𝑖𝑗 𝑅 ⎩ 𝐿 2 2 Where 𝜌̄𝑖𝑗 =(𝜌i +𝜌j )/2 and 𝑐̄𝑖𝑗 =(ci +cj )/2 are the mean values of the density and the speed of sound between particle i and j respectively. The interaction between the two particles is replaced by the Riemann solution as follows: ⎧ 1 (𝑝 + 𝑝 ) → 𝑃 ∗ 𝑗 𝑖𝑗 ⎪2 𝑖 (9) ⎨ ( ) 1 ∗ ⎪ 𝐯𝑖 + 𝐯𝑗 → 𝐯 𝑖𝑗 ⎩2 Through the transformation of Eq. (9), the continuity and momentum equation in Eq. (5) can be written as

∑( 𝑓 (𝐱𝑗 ) 𝑗

(4)

Where W(x − x′, h) and ∇W(x − x′, h) are smooth function and gradient respectively. h is the smooth length of the area affected by the kernel function, Ω is the support domain of the smooth function, Vj =mj /𝜌j is the volume of particle j, mj and 𝜌j represent the mass and density of particle j respectively. For the inviscid flow, the governing equations of the multi-phase flow are discretized in SPH method. Following the literature [26], the discrete forms of the governing Eq. (1) are as follows: ⎧ 𝑑 𝜌𝑖 ∑ = 𝜌𝑖 (𝐯𝑖 − 𝐯𝑗 ) ⋅ ∇𝑖 𝑊𝑖𝑗 (𝐫𝑖 − 𝐫𝑗 , ℎ)𝑉𝑗 ⎪ 𝑗 ⎪ 𝑑𝑡 ∑ ⎪ = 2𝜌𝑖 (𝐯𝑖 − 𝐯̄ 𝑖𝑗 ) ⋅ ∇𝑖 𝑊𝑖𝑗 (𝐫𝑖 − 𝐫𝑗 , ℎ)𝑉𝑗 ⎪ 𝑗 ⎨ 𝑑 𝐯𝑖 1 ∑ ⎪ =− (𝑝𝑖 + 𝑝𝑗 )∇𝑖 𝑊𝑖𝑗 (𝐫𝑖 − 𝐫𝑗 , ℎ)𝑉𝑗 + 𝐠 𝜌𝑖 𝑗 ⎪ 𝑑𝑡 ∑ ⎪ = − 𝜌2 𝑝̄𝑖𝑗 ∇𝑖 𝑊𝑖𝑗 (𝐫𝑖 − 𝐫𝑗 , ℎ)𝑉𝑗 + 𝐠 ⎪ 𝑖 𝑗 ⎩

(5)

⎧ d 𝜌𝑖 ∑ = 2𝜌𝑖 (𝐯𝑖 − 𝐯∗𝑖𝑗 ) ⋅ ∇𝑖 𝑊𝑖𝑗 𝑉𝑗 ⎪ 𝑗 ⎪ d𝑡 ⎨ d 𝐯 2 ⎪ 𝑖 = − ∑ 𝑝∗ ∇ 𝑊 𝑉 + 𝐠 𝑖 𝑖𝑗 𝑗 ⎪ d𝑡 𝜌𝑖 𝑗 ⎩

Where, 𝒗̄ 𝑖𝑗 is the average velocity and 𝑝̄𝑖𝑗 represents the average pressure between particle i and j. The selection of the kernel function is important for the SPH calculation. Some traditional kernel functions are liable to cause tensile instability, such as cubic spline and Gaussian kernel function. In order to avoid this problem, this paper uses the kernel function proposed by Wendland [41]: { 𝛼(2 − 𝑞 )4 (2𝑞 + 1), 0≤𝑞<2 𝑊 (𝐫, ℎ) = (6) 0, 𝑞≥2

(10)

In the process of hydrodynamic numerical simulation, the traditional SPH method with artificial viscous term has some disadvantages, such as non-physical viscous dissipation, poor calculation accuracy and poor numerical stability. The SPH method based on Riemann solution can optimize the non-physical pressure oscillation and introduce fewer empirical parameters, but which can introduce larger numerical viscosity. 136

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

In order to overcome these shortcomings, a limiter is introduced into Eq. (8) by Zhang and Hu [43], which can effectively study the impact problems of single-phase flow but cannot simulate the viscous multiphase flow problems due to the limit of its application scope (it cannot express the real physical viscosity of the fluid). This paper not only solves the numerical dissipation of SPH model based on Riemann solution, but also extends its application for the simulation of the viscous multiphase flow problems. To describe the physical viscous properties of multiphase flow in calculation, the second term on the right hand of Eq. (8) is improved to make it suitable for solving viscous multi-phase flows with large density ratios. The following conversion is made as: ⎧ 𝜇̄ ⎪𝐾 𝑖𝑗 (𝑢𝐿 − 𝑢𝑅 ) 1 𝜌̄ 𝑐̄ (𝑢 − 𝑢𝑅 ) → ⎨ ℎ 2 𝑖𝑗 𝑖𝑗 𝐿 ⎪0 ⎩

𝑟𝑖𝑗 > 0 𝑟𝑖𝑗 ≤ 0

2.5.1. The interface repulsive force In this paper, the interface force is added in the momentum equation to avoid the penetration of particles of different phase which is proposed by Monaghan [45]. The force is applied only in two particles coming from different fluids and acts in a band of width 2 h on either side of the interface. So far, it has been further extended to the problems of multiphase flow with all density ratios. The form of the force in the present is: 1 ∑ 𝐹𝑖𝐶 = 𝜁 (𝑃𝐿 + 𝑃𝑅 )∇𝑖 𝑊𝑖𝑗 𝑉𝑗 (17) 𝜌𝑖 𝑗 Where the factor 𝜁 is the interface force coefficient and presets between 0.01 and 0.1, following Grenier et al. [24] suggestion. The choice of the interface force coefficient value is related to the complexity of multiphase flow problems. Too small interface force coefficient would induce the interface to be penetrated and make the numerical results unreasonable. If the coefficient is much larger, it would affect the momentum conservation.

(11)

where, K=2.5(2+d) is the correction coefficient of viscous term according to the numerical examples, and d is the dimension; 𝜇̄ 𝑖𝑗 = 2𝜇 i 𝜇 j /(𝜇 i +𝜇 j ) is the inter-particle average dynamic viscosity, and 𝜇i is the dynamic viscosity of particle i. So the final pressure of intermediate state is determined as: 𝑃𝑖𝑗∗ =

𝑝𝑖 + 𝑝𝑗 2

+

𝐾 𝜇𝑖 𝜇𝑗 (𝑢 − 𝑢𝑅 ) ℎ 𝜇𝑖 + 𝜇𝑗 𝐿

2.5.2. The surface tension force In the present work, the continuum surface tension model (CSF) model [46] is used as a volumetric force and applied to different phase fluid particles near the interface. According to Morris et al. [47], the surface tension force per unit volume is described as follows:

(12)

⌢

𝐹 𝑠 = 𝜎𝜅 𝐧 𝛿 𝑠

Therefore, we convert the Riemann numerical dissipation term into the Riemann viscous term on the right hand of Eq. (12). The Riemann viscous term is similar to the discretized viscous formulae in the literature [39].

Where 𝜎 is the surface tension coefficient, 𝜅 is the curvature of the in⌢ terface, 𝐧 is the unit external normal of interface, and 𝛿 s is the surface delta function. In order to calculate the normal vector of the interface, the color function is used to label particles: { 0, 𝑖 ∈ 𝑙𝑖𝑞𝑢𝑖𝑑 𝐴 𝐶𝑖 = (19) 1 𝑖 ∈ 𝑙𝑖𝑞𝑢𝑖𝑑 𝐵

2.4. Solid wall conditions The distribution of the fixed virtual particles is independent of fluid particles. So this paper combines Adami et al. [44] virtual particle solidwall conditions with one-sided Riemann problem to impose solid-wall boundary of multiphase flow. Assuming that particle i is fluid particle and particle j is solid-wall boundary particle. The initial left and right states are written as (

( ) = 𝜌𝑖 , − 𝐯 𝑖 ⋅ 𝐧 , 𝑃 𝑖

(13)

( ) ( ) 𝑈𝑗 = 𝜌𝑅 , 𝑢𝑅 , 𝑃𝑅 = 𝜌𝑗 , −𝐯𝑗 ⋅ 𝐧, 𝑃𝐿 + 𝜌𝑖 (𝐚𝑤 − 𝐠) ⋅ 𝐫𝑖𝑗

(14)

𝑈 𝑖 = 𝜌𝐿 , 𝑢 𝐿 , 𝑃 𝐿

)

The color function is smoothed and the interface normal vector is obtained by applying the variational principle. 𝐧𝑖 =

𝑁 ∑ 𝑗=1

(𝐶̄𝑗 − 𝐶̄𝑖 ) ⋅ ∇𝑖 𝑊𝑖𝑗 𝑉𝑗

(20)

∑ Where 𝐶̄𝑖 = 𝑗 𝐶𝑗 𝑊𝑖𝑗 𝑉𝑗 . In order to avoid the influence of edge effect, the interface unit normal vector is as follow: { 𝐧𝑖 ∕||𝐧𝑖 ||, if ||𝐧𝑖 || > 𝜀, ̂ (21) 𝐧𝑖 = 0, 𝑜𝑡ℎ𝑒𝑟𝑤𝑖𝑠𝑒 Where, the parameter 𝜀 is set to 0.01 h, which is obtained from numerical experiments. The curvature is expressed by normal divergence and modified by CSPM [48] as: ∑𝑁 ̂ ̂ ( ) 𝑗=1 min(𝑁𝑖 , 𝑁𝑗 )(𝐧𝑗 − 𝐧𝑖 ) ⋅ ∇𝑖 𝑊𝑖𝑗 𝑉𝑗 𝜅 =− ∇⋅̂ 𝐧 𝑖= (22) ∑𝑁 𝑗=1 min(𝑁𝑖 , 𝑁𝑗 )𝑊𝑖𝑗 𝑉𝑗

where, n is the normal unit vector of the solid-wall boundary, vj is the wall velocity and aw is the wall acceleration. Both free-slip and no-slip boundary conditions can be implemented by the choice of the Riemann viscous term in Eq. (12). In the former case, the viscous term is ignored. While, in the latter case, the right state velocity uR in Eq. (12) is evaluated as follows: 𝑢𝑅 = 𝐯𝑖 ⋅ 𝐧 + 2𝐯𝑗 ⋅ 𝐧

(18)

(15)

𝛿 s is used to characterize the distribution of the surface tension. In order to satisfy the continuity condition, 𝛿 s =|n|. Finally, the surface tension of particle i is obtained as:

The right state density is obtained by solving the state equation in Eq. (2). As for the fluid pressure extension, the pressure of wall particles Pw can be calculated as follows: ∑ 𝑓 ∈𝑓 𝑙𝑢𝑖𝑑 𝑝𝑓 𝑊𝑤𝑓 𝑃𝑤 = ∑ (16) 𝑓 ∈𝑓 𝑙𝑢𝑖𝑑 𝑊𝑤𝑓

⌢

𝐹𝑖𝑠 = −𝜎(∇ ⋅ 𝐧 )𝑖 𝐧𝑖 ∕𝜌𝑖

(23)

For the viscous multiphase flow, the final continuous equation and momentum equation are written as ⎧ d 𝜌𝑖 ∑ = 2𝜌𝑖 (𝐯𝑖 − 𝐯∗𝑖𝑗 ) ⋅ ∇𝑖 𝑊𝑖𝑗 𝑉𝑗 ⎪ 𝑗 ⎪ d𝑡 ⎨ ⌢ d 𝐯 ∑ ∗ 𝜁 ∑ 2 𝜎 𝑖 ⎪ =− 𝑝 ∇𝑖 𝑊𝑖𝑗 𝑉𝑗 − (∇ ⋅ 𝐧 )𝑖 𝐧𝑖 + (𝑃 + 𝑃𝑅 )∇𝑖 𝑊𝑖𝑗 𝑉𝑗 + 𝐠 ⎪ d𝑡 𝜌𝑖 𝑗 𝜌𝑖 𝜌𝑖 𝑗 𝐿 ⎩ (24)

2.5. The treatment for the multiphase interface In order to keep sharp interface and consider the influence of the interface scale, there are a small interface force and the surface tension force at the interface. 137

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 2. Oscillation of Square-droplet: (a) the initial setup of square bubble. The evolution of square droplet (b) t==0.0, (c) t==2.0 s and (d) t==12.0 s.

2.6. The time step In order to ensure the stability and accuracy of numerical calculation, the time step needs to meet three restrictions: limitations of the CourantFriedrichs-Lewy (CFL) condition, the viscosity model condition, and the surface tension model condition [39,49]. ( ) ⎧ ℎ ⎪Δ𝑡 ≤ min 𝑐 + 𝑣max ⎪ ( ⎪ 2) ⎪Δ𝑡 ≤ 0.125 min 𝜌𝑖 ℎ (25) ⎨ 𝜇𝑖 ⎪ ⎪ ( )1∕2 𝜌𝑖 ℎ 3 ⎪ Δ𝑡 ≤ 0 . 5 min ⎪ 2𝜋𝜎 ⎩ Where, vmax is the maximum particle velocity in the calculation. The velocity-Verlet time-stepping scheme [44] is applied in the present work.

Fig. 3. Time evolution of pressure different between the inside and outside droplet pressure. Table 1 Physical parameters used in the numerical single bubble rising.

3. Validation of the new model to solve multi-fluid problems In this section, several cases of multi-fluid problems are presented to test the accuracy, stability and applicability of the proposed R-MSPH model. The first case is the 2D square-droplet deformation which is used to verify the effectiveness of the interface force and the surface tension in the model. The second case is the single bubble rising which is devoted to demonstrate the ability of the R-MSPH model for the simulation of multiphase flows with different density and viscosity ratios. The third case is the double bubbles rising, in which the ability to deal with complex interface deformation can be verified. The SPH results are compared with the existing reference results in the literature.

case

𝜌𝑎 (𝑘𝑔∕𝑚3 )

𝜌𝑏 (𝑘𝑔∕𝑚3 )

𝜂𝑎 (Pa ⋅ s)

𝜂𝑏 (Pa ⋅ s)

𝑔 (𝑚∕s2 )

𝛽 (𝑁∕m)

Re

Bo

1 2

1000 1000

100 1

10 10

1 0.1

0.98 0.98

24.5 1.96

35 35

10 125

subjected to surface tension force. Driven by surface tension force, the square droplet gradually becomes smooth (Fig. 2(c)), and finally evolves into a stable circle droplet (Fig. 2(d)). Based on the Laplace law, when the droplet reaches the steady equilibrium, the pressure difference between phase A and phase B is a cer√ tain value 𝑝̄𝐴 − 𝑝̄𝐵 = 𝜎 𝜋∕𝐿. Fig. 3 shows the evolution process of the pressure difference over time under different particle spacing. With the evolution of square droplet, the pressure difference oscillates up and down, and the amplitude √ becomes smaller and smaller. Finally, it tends to be a stable value 𝜎 𝜋∕𝐿 = 1.77, which conforms to the Laplace’s law. It is found that the numerical solutions are convergent and the present surface tension is accurate and effective in simulating multiphase flows with large density and viscosity ratios.

3.1. Oscillation of 2D square-droplet The oscillation of 2D square-droplet problem is a typical benchmark for validating the model of the interface repulsive force and the surface tension force in zero gravity condition. The case has been adopted by Adami et al. [25]. The initial condition is shown as Fig. 2(a). A 2D square-droplet L=1 m each side filled with phase A is statically placed in the square liquid tank 2 L each side filled with phase B. The density radio of the both fluid is 𝜌B /𝜌A =1000 (𝜌B =1000 kg/m3 ). The viscosity radio is 𝜇 B /𝜇 A =100 (𝜇 B =0.5Pa ⋅s). The surface tension coefficient is 𝜎=0.02 N/m. Three kinds of particle configuration are 50 × 50, 100 × 100 and 150 × 150 respectively to simulate the evolution of square droplet. The background pressure is p0 =50 Pa and the interface force coefficient is 𝜁 =0.08. The all sides of the square tank is set as free-slip boundary condition. The evolution snapshots of the oscillating droplet over time are shown in Fig. 2(b)–(d). Because the surface tension is proportional to the curvature of the interface, the four corners of the droplet are first

3.2. Single bubble rising through a vertical column of water Two cases of a single bubble rising (the same research see e.g. [24,29]) with different density and viscous ratios are simulated to verify the ability of the R-MSPH model to deal with multiphase problems with different density and viscosity ratios. The bubble (R=0.25 m) placed in a cylindrical tank floats up under buoyancy. The initial setup is shown as Fig. 4(a). Table 1 summarizes the physical parameters of the cases. 138

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

√ Fig. 4. Single bubble rising: (a) the initial setup. The pressure distributions of the present method and the comparisons of the numerical results at 𝑡 𝑔∕𝐷 = 3: (b) the case 1; (c) the case 2 (left half represents the R-MSPH results and right half shows the FEM results [50]).

Fig. 5. Comparisons about the variation of bubble’s center with time: (a) the case 1; (b) the case2.

The Reynolds number and Bond number are 𝑅𝑒 =

√ 8𝜌2𝑎 𝑔 𝑅3 ∕𝜂𝑎 and

Bo = 4g𝜌a R2 /𝛽. The interface force coefficient is 𝜁 =0.08. The two cases are carried out with the same particle resolution (100×200). Fig. 4(b) and (c) show the comparisons of the simulation results of case 1 and case 2 between the present R-MSPH and the FEM in literature [50] respectively. The left half each figure represents the cloud chart of pressure distribution and the bubble shape predictions of √ the present method at 𝑡 𝑔∕𝐷 = 3, and the right half show the predictions of the FEM method. It is observed that the pressure distribution is uniform and stable, and the present results are in good agreement with the FEM results [50]. In order to further analyze and study the results of this paper, Fig. 5(a) and (b) show the curves with time of the center of mass of case 1 and case 2 respectively, which are compared with those calculated by FEM method [50]. It is found that the numerical solutions are also convergent and are in good coincidence with FEM results. 3.3. Interaction of two rising bubbles through a fluid column In this section, the case of two rising bubbles presented in Rezavand and Taeibi-Rahni [29] is simulated, which involves a series of complex

Fig. 6. The schematic of two rising bubbles.

139

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Table 2 Physical parameters and gravitational constant of two rising bubbles.

Table 3 Calculation parameters of the steady jet simulation.

𝜌𝑎 (kg∕m3 )

𝜌𝑏 (kg∕m3 )

𝜂𝑎 (Pa ⋅ s)

𝜂𝑏 (Pa ⋅ s)

𝑔 (𝑚∕s2 )

1000

100

0.156

0.078

9.8

dynamic evolution processes such as tearing, merging and re-tearing of bubbles. The two bubble rising case is a challenging test case to validate the ability of multiphase flow model to deal with the multiphase flow problems with complex interface deformation. Grenier et al. [49] also used this case to validate their two-phase ISPH scheme. The schematic of this problem is illustrated in Fig. 6 and physical parameters are showed in Table 2. The Reynolds number using small bubble radius is Re=633, the interface force coefficient is 𝜁 =0.08, and the 320 × 480 particles discrete the computational domain. Considering the influence of surface tension on the interface of bubbles, Bond number Bo=∞ and Bo=80 are used to simulate the case. We observe the dynamic evolution process of two rising bubbles. Figs. 7 and 8 show the numerical predictions of bubble particle positions in six different time instants for Bo=∞ and Bo=80 respectively, and the Level-Set results come from Grenier et al. [49]. The right side of each image is the prediction of bubble shape in this paper, and the left side is the Level-Set results. At first, two bubbles rise simultaneously under the action of buoyancy and evolve into mushroom shape. Then the lower bubble is lengthened and swallowed down by the upper one. Two tail-toe bubbles separate from the lower bubble and rise with the lower one, and the liquid film thickness between two bubbles becomes thinner and thinner. The two bubbles merge and tear apart from the symmetry at last. Comparing Fig. 7 with Fig. 8, it can be seen that the surface tension has a major impact on the evolution of rising bubble shape. The surface tension causes the interface more rounded and restrains the generation of large vortices of upper bubble. Though local differences are visible in the tail toes of the lower bubble, a good agreement is observed with Level-Set results. Level-Set method is restricted by grids, and the tail-toe bubbles keep connection with the lower bubble. While the results of this paper agree very well with those obtained by SPH method from reference [29,49]. The simulation of two rising bubbles in this section confirms the applicability of the R-MSPH model in dealing with complex multiphase flow problems. The complex evolution process of the interface can also be clearly and steadily captured.

We

Re

Fr

Oh

1 2 3 4

12.56 17.76 21.75 25.12

1.0 2.0 3.0 4.0

11.09 15.68 19.20 22.18

0.57 0.8 0.98 1.13

0.09 0.09 0.09 0.09

Where Re = 𝜌𝑣𝐷 is the Reynolds number, which represents the ratio of 𝜇 inertial force to viscous force. And, the dimensionless Froude number expresses the ratio of inertial force to gravity: 𝑣 𝐹𝑟 = √ (28) 𝑔𝐷 The period and amplitude of the length of the jet liquid column are stochastic variable, but they can fluctuate within a certain range. When the internal force remains the dominant over the gravity, an predicted formula of the jet breakup length by experimental modification is referred in the literature [49] as follow: √ 𝐿𝑏 = 13 𝑊 𝑒(1 + 3𝑂ℎ) (29) 𝐷 Where, Lb is the jet breakup length, and the definition of jet breakup length is the average value of liquid column length for successive breakup periods. The recommended formula is provided for reference in the simulation results. 4.1. Steady jet The inlet velocity of the steady jet does not change with time. In literature [20], the jet breakup problems with different inlet velocity are simulated and analyzed using the SPH method, but the influence of the surrounding air is not considered. In this section, using the RMSPH model presented in this paper, the four numerical cases in literature [20] are simulated again. The parameters of the cases are shown in Table 3. The liquid jet with density 𝜌l =1120 kg/m3 injected into the air-filled container through the nozzle with diameter D=0.005 m. The jet viscosity is 𝜐l =0.0316 kg/(ms), the air density is 𝜌a =1 kg/m3 , the air viscosity is 𝜐a =0.0177 × 10−3 kg/(ms), the surface tension coefficient is 𝜎=44mN/m, and the interface force coefficient is 𝜁 =0.08. The gravity direction is the same as that of the jet entering into container and the gravitational acceleration is g=9.81 m/s2 . The initial particle spacing is Δx=0.0002 mm, and the independence check of particle spacing has been tested. At present, it is the largest and most stable particle spacing. Figs. 10–13 illustrate the behaviors of the jet breakup with different inlet velocity which are obtained by using the R-MSPH model. The liquid jet injects the static air at a certain speed, and the liquid column is stretched by inertia force, while the viscous force in liquid and the surface tension at the interface hinder the increase of the liquid column length. As time goes on, the jet head expands, the liquid column length elongates and the inertia force increases. When the inertia force is about to overcome the effect of viscous force and surface tension, the volume of the jet head is the largest and the neck connected to the upstream is the thinnest. At that time, the jet core length reaches the maximum in the period. The jet head falls off into a droplet, the liquid column shrinks temporarily under the action of surface tension and viscous force. It shrinks to the minimum value in the period, and then

In this section, the R-MSPH model proposed in this paper is adopted to study the jet breakup problems of steady jets or pulsating jets. The influencing factors of jet breakup length, such as jet velocity, gravity, surface tension, amplitudes and period of pulse velocity, are comprehensively compared and studied. A liquid jet is injected into air from the nozzle of the width D, and the inlet velocity is v0 . In order to reduce the influence of the calculation domain boundary on the jet breakup, the numerical examples indicate that the width of the calculation domain is 20D.The length of the domain is long enough to ensure the integrity of the jet breakup. Ignoring the viscous friction between the nozzle and the jet, the upper and lower boundaries are set as the entrance and exit boundaries respectively. The computational setup is shown in Fig. 9. Under the gravity, inertia force, viscous force and surface tension, the jet breaks up and produces droplets when the jet reaches a certain length. Weber proposed a dimensionless Weber number (the ratio of aerodynamic force to surface tension) characterizing the jets to analyze the mechanism of the jet breakup: 𝜌𝑣2 𝐷 2𝜎

V0(m/s)

Ignoring aerodynamic force and considering the effects of liquid density, viscosity, surface tension and nozzle diameter on the jet breakup, Ohnesorge proposed a dimensionless Ohnesorge number to characterize the jet: √ 𝑊𝑒 𝑂ℎ = (27) Re

4. Results and discussion of jet breakup problems

𝑊𝑒=

case

(26) 140

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 7. Two rising bubble shape evolution for Bo =∞ (marked by red line on the left half of each panel: Level-Set simulation [49], right half of each panel particle-filled results: present R-MSPH model). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.) Table 4 Comparison of the breakup length between the present work.

begins the next period. Later the peak-to-valley difference of liquid column length maintains in a certain range. In the Figs. 10–13, the surface wave structure of the liquid jet can be seen clearly and the evolution process of jet is consistent with the jet theory. The above results are compared with those of literature and theory, as shown in the Table 4. The jet breakup length increases with the increase of the inlet velocity, which is same as the literature [20]. In this paper, considering the existence of air around jet, there is no lack of particles in the support domain of the interface particles, and the aerodynamic force is also taken into account. The calculations of the surface tension are more accurate, and the results are closer to the predicted by the linear theory results [51] than the literature [20]. To further research the influence of surface tension on the stable jet breakup problem, four other cases with different surface tension coeffi-

case

Lb /D(R-MSPH)

Lb /D [20]

Lb /D [51]

1 2 3 4

14.4 22.9 28.0 32.1

13.2 21.1 27.5 31.2

16.5 23.4 28.6 33.0

cient are simulated. The parameters setting are same as case 2 except the value of surface tension coefficient. We change the surface tension coefficient in the simulation (the properties of liquid jet remain unchanged, the surface tension can be changed by adding a small amount of surfactant to the jet fluid in the experiment), and the 𝜎 is 33mN/m, 88 mN/m, 132 mN/m and 176 mN/m separately. 141

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 8. Two rising bubble shape evolution for Bo=80 (marked by red line on the left half of each panel: Level-Set simulation [49], right half of each panel particle-filled results: present R-MSPH model). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

the above discussion shows with the inlet velocity increasing, the Weber number and Reynolds number increase and the jet breakup length increases. And, the Weber number is inversely proportional to the surface tension, the result shows with the surface tension increasing, the Weber number decreases and the jet breakup length also decreases. In conclusion, the breakup length of the steady jet is closely connected with the Weber number. With the increase of the Weber number, the jet breakup length has a trend of increase by degrees. It is in good agreement with the linear theory [51].

Fig. 14 shows the primary breakup behaviors, and Fig. 15 shows the jet breakup behaviors during the stable stage with five different surface tension coefficients. It is observed that the primary breakup length and the jet breakup length decrease, primary droplet radius increases as the surface tension increases. The jet breakup length increases tends to gently and the uniform drops are produced at 𝜎=176 mN/m. The direction of surface tension is same as the normal direction of jet surface, and the surface tension enables the liquid surface to contract inward. It can be concluded that the surface tension prevents the liquid jet extending downward, increasing the surface tension of liquid jet is beneficial to obtaining smooth and uniformly distributed droplets. Accurate application of surface tension plays an important role in jet breakup behaviors. Because the Weber number is proportional to the inlet velocity and the Reynolds number is proportional to the square of the inlet velocity;

4.2. Pulsating jet For the jet breakup problems, the inlet velocity has a great influence on the liquid column breakup length. However, among all liquid 142

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 9. The setup of jet breakup problem. Fig. 11. Computed snapshots of the steady jet case 2.

Fig. 10. Computed snapshots of the steady jet case 1. Dark blue particles are liquid and the surrounding green particles are air. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

jet researches, there are far more steady jet studies than the pulsating jets. In this section, the numerical simulation of pulsating jets is carried out. The inlet velocity of pulsating jet varies periodically with time as following: Fig. 12. Computed snapshots of the steady jet case 3.

𝑣 = 𝑣0 + 𝐴𝑚 sin(2𝜋𝑡∕𝑇 )

(30)

Where v0 is the average velocity, Am is the amplitude of inlet velocity, and T is the period of inlet velocity. We change the period and frequency of inlet velocity and study the pulsating jet breakup behaviors in this section. First of all, the case of the pulsating jet of the literature [52] is simulated by the R-MSPH model. The same jet parameters as the jet experiment are used as the following parameters at 13°C: density is 𝜌l =1000 kg/m3 , kinetic viscosity is 𝜐 =1.20 × 10−6 m2 /s, surface tension coefficient is 𝜎 =73.77 × 10−3 N/m, and the interface force coefficient is 𝜁 =0.08. The physical properties of the surrounding air are same as the Section 4.1. The jet average velocity is v0 =0.39 m/s (Re=890), velocity amplitude is Am =0.1 m/s, and pulsating period is T0 =50 ms. The nozzle width of D =3 mm and the gravitational acceleration of g = 9.81 m/s2 are adopted.

Fig. 16 shows the simulation results of the pulsating jet injecting into the static air in the direction of gravity at Am =0.1 m/s and T0 =50 ms. Each image in Fig. 16 represents a snapshot of the jet breakup, there is no fixed time interval. In order to have a better comparison, the time t/T0 =0 is the starting time of the literature [52]. At the beginning, gravity is the dominant mechanism for jet breakup, and the liquid column breakup length is getting longer and longer over time. The first breakup length of liquid column is Lb /D =16.4 for this case and the largest droplet is formed, the second breakup length is Lb /D=32.9. As time goes by, the driving effect of the nozzle becomes more and more important, the end-pinching-off behaviors of pulsating jet display periodic variation, and the axisymmetric surface waves and satellite droplets are clearly observed. In the Fig. 16, three jet breakup periods are shown. 143

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 15. The breakup behaviors of the steady jet during the stable stage.

tained by the SPH methods are larger than the experimental results. But the peak-to-valley difference of liquid column length in the SPH results nearly equals the experiment: Lb /D=14.4 in the no-air SPH, Lb /D=13.5 in the R-MSPH and Lb /D=12.7 in the experiment. Because the inlet velocity of the pulsating jets varies periodically as the time, as well as the values of the Weber number and the Reynolds number can be increased by increasing inlet velocity. In the following sections, the influence of the amplitude and the period of inlet velocity on pulsating jet breakup length is studied.

Fig. 13. Computed snapshots of the steady jet case 4.

4.2.1. Effect of the amplitude of velocity on pulsating jet breakup length In this section, combined with the case of Am =0.1 m/s, four more cases, the average velocity of Am =0.025 m/s, 0.05 m/s, 0.075 m/s and 0.125 m/s (the other parameters remain the same), are simulated and the influence of the amplitude on the jet breakup behaviors is analyzed. Fig. 18 shows the time histories of liquid column length with different amplitude of inlet velocity. It can see that the liquid column length of five cases displays periodic variation, and their periods are approximately equal to T0 =50 ms. Along with the velocity amplitude increasing, the end-pinching-off length of liquid jet shortens and the values are Lb /D=37.13, 28.45, 20.05, 20.03 and 17.78 respectively. By increasing the amplitude from Am =0.025 m/s to Am =0.05 m/s, it is easy to obtain the jet that the peak-to-valley difference of liquid column length decreases and the droplets produced are evenly distributed. When the velocity amplitude is greater than 50 m/s, the end-pinching-off period is strictly equal to T0 =50 ms. And, as the velocity amplitude continue to increase, the end-pinching length increases slightly and the peak-to-valley length is almost equal to Lb /D =12.1. To study the breakup characteristics of pulsating jets in detail, Fig. 19 shows the breakup behaviors of pulsating jets with different amplitude of inlet velocity at T0 =50 ms during the jet stable stage. The case of Am =0.05 m/s generates smooth and uniform droplets. When the velocity amplitude is greater than 0.05 m/s, the pulsating jet breaks up and forms larger diameter droplet accompanied by satellite droplets. With the drop of large droplet, the droplets change from sphere to oblate form.

Fig. 14. The primary breakup behaviors of the steady jets.

Fig. 17 illustrates the history of liquid column length at Am =0.1 m/s and T0 =50 ms. The breakup behaviors can be seen in Fig. 16. Results from present R-MSPH simulation are compared with the existing experimental results and no-air SPH results from Takashima et al. [52]. It is can be seen that the change periods of the history curves of liquid column length from 3 approaches in general are agreeable. The endpinching-off period is t/T0 =1 in the steady time, which is equal to the driving period of the nozzle. The end-pinching-off length of liquid jet from long to short is the no-air SPH (Lb /D=22.8), the present work (Lb /D =15.1) and the experiment (Lb /D = 10.1). This can also demonstrate again that using the SPH method to simulate the jet breakup problems, it is important to consider the influence of the air around the jet. Since the results of no-air SPH and the present are obtained by twodimensional simulation, the values of the end-pinching-off length ob144

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

Fig. 16. Computed snapshots of the pulsating jet into air at Am =0.1 m/s and T0 =50 ms.

Fig. 17. Comparisons of the time history of liquid column length at Am =0.1 m/s and T0 =50 ms.

Fig. 19. The breakup behaviors of pulsating jets at T0 =50 ms.

The velocity amplitude can change the length of the surface wave of the pulsating jet, thus affecting the jet breakup length and making the jets prone to parametric resonance phenomenon. The research conclusion will help to obtain the droplets with low temporal resolution and uniform distribution or with high temporal resolution and breakup regulation by changing the amplitude of inlet velocity. 4.2.2. Effect of the period of velocity on pulsating jet To investigate the effect of the velocity period on pulsating jet, the problem is solved with five different velocity period cases namely, T=10 ms, 25 ms, 50 ms, 75 ms and 100 ms (other parameters remain

Fig. 18. Time histories of liquid column length at T0 =50 ms.

145

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

lite droplets is increasing, and the length of jet surface waves becomes longer and longer under the influence of the driving force of the jet inlet. Therefore, under a given value of velocity amplitude, the breaking behaviors of pulsating jet are controlled by the length of the surface waves and pulsating disturbance of the jet, the end-pinching-off length could decrease first and then increase with the increasing of the velocity period. 5. Conclusion In this paper, a new SPH model for multiphase fluids with large density and viscosity ratios based on Riemann solvers is proposed. In order to simulate the viscous multiphase flows and reduce the Riemann dissipation in the model, we convert the Riemann numerical dissipation term into a Riemann viscous term that can describe the physical viscosity of the multiphase flow. In order to ensure the interfacial stability, the model includes surface tension and repulsive force between different phases near the interface. The numerical stability is improved even if artificial dissipation techniques are not added, and the test of the tuning parameters of artificial dissipation technology for different cases is reduced. Different multiphase tests have been conducted to verify the stability and accuracy of the R-MSPH model. After validation, the jet breakup problems of steady jets or pulsating jets are simulated and analyzed. The numerical results show the following:

Fig. 20. Time histories of liquid column length at Am =0.1 m/s.

(1) In the numerical study of the liquid jet breakup problems, the influence of the surrounding air cannot be neglected, especially the SPH method is used to do it. (2) The droplets with low temporal resolution and uniform distribution can be obtained by changing velocity amplitude of pulsating jet. (3) The velocity period has a great influence on the end-pinching breaking length and the peak-to-valley breaking length. In conclusion, the present multi-phase SPH model can simulate the jet breakup problems very well. The results are benefit to design and analyze predictive control in the field of jet application. By controlling jet velocity, surface tension, driving velocity amplitude and velocity period, the ideal injection results can be achieved. However, due to the limitation of computational scale, we only study the low-speed twodimensional jet breakup problems, the studies of the breakup mechanism and characteristic of the high-speed three-dimensional liquid jet are open in future works. Adaptive particle refinement and derefinement are necessary to reduce the computational cost. Fig. 21. The breakup behaviors of the pulsating jets at Am =0.1 m/s.

Acknowledgments This work was supported by the National Nature Science Foundations of China (Grant no. 11972309), the National Natural Science Foundation for Young Scientists of China (Grant no. 11702220), the Aeronautical Science Foundation of China (Grant no. 20162353022), and the Fundamental Research Funds for the Central Universities (Grant no. 3102017zy066).

the same as Section 4.2). Fig. 20 shows time histories of liquid column length for different velocity period at Am =0.1 m/s. The results of five cases show that the peak-to-valley periods of pulsating jets are almost equal to their driving periods. As the velocity period increases, the peak-to-valley length increases, the end-pinching-off length first decreases then increases, and the end-pinching-off length has a minimum value near the period T0 =25 ms. It should be noted that the pulsating jet T=100 ms breaks up several times in one period, which is beneficial to the formation of satellite droplets, as shown the Fig. 21. Fig. 21 shows the breakup behaviors of pulsating jets with different period of inlet velocity at Am =0.1 m/s during the jet stable stage. From the picture, it can be seen that the droplets produced by the case with velocity period T = 10 ms are small and uniformly distributed. When the velocity period is T = 25 ms, the droplet diameter and center spacing become larger than that of T = 10 ms, and the surface waves generated by jet driving are clearly observed. When the velocity period is larger than T = 25 ms, the generation of main droplets is accompanied by some satellite droplets. When continue to increase the velocity period, the diameter of main droplets is getting larger and larger, the number of satel-

References [1] Som S, Aggarwal SK. Effects of primary breakup modeling on spray and combustion characteristics of compression ignition engines. Combust Flame 2010;157(6):1179– 93. doi:10.1016/j.combustflame.2010.02.018. [2] Marmottant P, Villermaux E. On spray formation. J Fluid Mech 2004;498(498):73– 111. doi:10.1017/S0022112003006529. [3] Huang JCP. The break-up of axisymmetric liquid sheets. J Fluid Mech 1970;43(2):305–19. doi:10.1017/S0022112070002392. [4] Rayleigh L. On the instability of jets. Proceedings of the London Mathematical Society; 1878. [5] Tomotika S. On the instability of a cylindrical thread of a viscous liquid surrounded by another viscous liquid. Proc R Soc London, Ser A 1935;150(870):322–37. doi:10.1098/rspa.1935.0104. [6] Meister BJ. Prediction of jet length in immiscible liquid systems. AIChE J 1969;15(5):689–99. doi:10.1002/aic.690150512. 146

Q. Yang, F. Xu and Y. Yang et al.

Engineering Analysis with Boundary Elements 111 (2020) 134–147

[7] Lin SP, Woods DR. A branching liquid jet. Phys Fluids 1994;6(8):2671. doi:10.1063/1.868156. [8] Richards JR, Lenhoff AM, Beris AN. Dynamic breakup of liquid–liquid jets. Phys Fluids 1994;6(8):2640. doi:10.1063/1.868154. [9] Pan Y, Suga K. A numerical study on the breakup process of laminar jets into a gas. Phys Fluids 2006;18(5):052101. doi:10.1063/1.2194936. [10] Saito S, Abe Y, Koyama K. Lattice boltzmann modeling and simulation of liquid jet breakup. Phys Rev E 2017;96(1):013317. doi:10.1103/PhysRevE.96.013317. [11] Ménard T, Tanguy S, Berlemont A. Coupling level set/VOF/ghost fluid methods: Validation and application to 3D simulation of the primary break-up of a liquid jet. Int J Multiph Flow 2007;33(5):510–24. doi:10.1016/j.ijmultiphaseflow.2006.11.001. [12] Fuster D, Bagué Anne, Boeck T, et al. Simulation of primary atomization with an octree adaptive mesh refinement and VOF method. Int J Multiph Flow 2009;35(6):550– 65. doi:10.1016/j.ijmultiphaseflow.2009.02.014. [13] Gingold RA, Monaghan JJ. Smoothed particle hydrodynamics: theory and application to non-spherical stars. Mon Not R Astron Soc 1977;181(3):375–89. doi:10.1093/mnras/181.3.375. [14] Koshizuka S, Oka Y. Moving-Particle semi-implicit method for fragmentation of incompressible fluid. Nucl Sci Eng 1996;123(3):421–34. doi:10.13182/NSE96-A24205. [15] Liu MB, Liu GR. Smoothed particle hydrodynamics (SPH): an overview and recent developments. Arch Comput Methods Eng 2010;17(1):25–76. doi:10.1007/s11831-010-9040-7. [16] Liu MB, Li SM. On the modeling of viscous incompressible flows with smoothed particle hydro-dynamics. J Hydrodyn Ser B 2016;28(5):731–45. doi:10.1016/s1001-6058(16)60676-5. [17] Zhang AM, Sun PN, Ming FR, et al. Smoothed particle hydrodynamics and its applications in fluid-structure interactions. J Hydrodyn Ser B (English Ed) 2017;29(2):187–216 CNKI:SUN:SDYW.0.2017-02-002. [18] Shibata K, Koshizuka S, Oka Y. Numerical analysis of jet breakup behavior using particle method. J Nucl Sci Technol 2004;41(7):715–22. doi:10.1080/18811248.2004.9715538. [19] Ganzenmüller S, Nagel A, Holtwick S, et al. Object-Oriented SPH-simulations with surface tension. High Perform Comput Sci Eng 2007;6:69–82. doi:10.1007/978-3-540-36183-1_6. [20] Sirotkin FV, Yoh JJ. A new particle method for simulating breakup of liquid jets. J Comput Phys 2012;231(4):1650–74. doi:10.1016/j.jcp.2011.10.020. [21] Hoefler C, Braun S, Koch R, et al. Modeling spray formation gas turbines — a new meshless approach. J Eng Gas Turbines Power 2013;135(1):011503. doi:10.1115/1.4007378. [22] Pourabdian M, Omidvar P, Morad MR. Numerical simulation of liquid jet breakup using smoothed particle hydrodynamics (SPH). Modares Mech Eng 2016;16(3):55– 66. doi:10.1142/S0129183117500541. [23] Colagrossi A, Landrini M. Numerical simulation of snterfacial flows by smoothed particle hydrodynamics. J Comput Phys 2003;191(2):448–75. doi:10.1016/S0021-9991(03)00324-3. [24] Grenier N, Antuono M, Colagrossi A, et al. An hamiltonian interface sph formulation for multi-fluid and free surface flows. J Comput Phys 2009;228(22):8380–93. doi:10.1016/j.jcp.2009.08.009. [25] Adami S, Hu XY, Adams NA. A new surface-tension formulation for multi-phase sph using a reproducing divergence approximation. J Comput Phys 2010;229(13):5011– 21. doi:10.1016/j.jcp.2010.03.022. [26] Monaghan JJ, Rafiee A. A simple sph algorithm for multi-fluid flow with high density ratios. Int J Numer Methods Fluids 2013;71(5):537–61. doi:10.1002/fld.3671. [27] Chen Z, Zong Z, Liu M B, et al. An sph model for multiphase flows with complex interfaces and large density differences. J Comput Phys 2015;283:169–88. doi:10.1016/j.jcp.2014.11.037. [28] Lind SJ, Stansby PK, Rogers BD. Incompressible-compressible flows with a transient discontinuous interface using smoothed particle hydrodynamics(SPH). J Comput Phys 2016;309:129–47. doi:10.1016/j.jcp.2015.12.005. [29] Rezavand M, Taeibi-Rahni M, Rauch W. An isph scheme for numerical simulation of multiphase flows with complex interfaces and high density ratios. Comput Math Appl 2018;75(8):2658–77. doi:10.1016/j.camwa.2017.12.034.

[30] Zheng BX, Chen Z. A multiphase smoothed particle hydrodynamics model with lower numerical diffusion. J Comput Phys 2019;382:177–201. [31] Inutsuka S. Godunov-type SPH. Memorie Della Societa Astronomica Italiana 1994;65:1027–31. [32] Monaghan JJ. SPH and Riemann solvers. J Comput Phys 1997;136(2):298–307. doi:10.1006/jcph.1997.5732. [33] Parshikov AN. Application of a solution to the riemann problem to the SPH method. Comput Math Math Phys 1999;39(7):1173–82. [34] Sirotkin FV, Yoh JJ. A smoothed particle hydrodynamics method with approximate riemann solvers for simulation of strong explosions. Comput Fluids 2013;88(12):418–29. doi:10.1016/j.compfluid.2013.09.029. [35] Rogers D, Dalrymple A, Stansby K. Simulation of caisson breakwater movement using 2-D SPH. J Hydraul Res 2010;48:135–41. doi:10.1080/00221686.2010.9641254. [36] Cai Z W, Zong Z, Chen Z, et al. Multiphase Godunov-type smoothed particle hydrodynamics method with approximate Riemann solvers. Int J Comput Methods 2018;1846010. [37] Puri K, Ramachandan P. Approximate riemann solvers for the Godunov SPH (GSPH). J Comput Phys 2014;270(8):432–58. doi:10.1016/j.jcp.2014.03.055. [38] Rezavand M, Zhang C, Hu XY. A weakly compressible sph method for violent multi-phase flows with high density ratio. J Comput Phys 2019:109092. doi:10.1016/j.jcp.2019.109092. [39] Ming FR, Sun PN, Zhang AM. Numerical investigation of rising bubbles bursting at a free surface through a multiphase SPH model. Meccanica 2017:1–20. doi:10.1007/s11012-017-0634-0. [40] Monaghan JJ. Smoothed particle hydrodynamics. Rep Prog Phys 2005;68(8):1703– 59. doi:10.1088/0034-4885/68/8/R01. [41] Wendland H. Piecewise polynomial, positive definite and compactly supported radial functions of minimal degree. Adv Comput Math 1995;4(1):389–96. doi:10.1007/BF02123482. [42] Toro EF. Riemann solvers and numerical methods for fluid dynamics. Berlin: Springer; 2013. p. 87–114. doi:10.1007/b79761_3. [43] Zhang C, Hu XY. A weakly compressible sph method based on a low-dissipation Riemann solver. J Comput Phys 2017;335:605–20. doi:10.1016/j.jcp.2017.01.027. [44] Adami S, Hu XY, Adams NA. A generalized wall boundary condition for smoothed particle hydrodynamics. J Comput Phys 2012;231(21):7057–75. doi:10.1016/j.jcp.2012.05.005. [45] Monaghan JJ. SPH without a tensile instability. J Comput Phys 2000;159:290–311. doi:10.1006/jcph.2000.6439. [46] Lafaurie B, Nardone C, Scardovelli R, et al. Modelling merging and fragmentation in multiphase flows with surfer. J Comput Phys 1994;113(1):134–47. doi:10.1006/jcph.1994.1123. [47] Morris JP. Simulating surface tension with smoothed particle hydrodynamics. Int J Numer Methods Fluids 2000;33(3):333–53 3<333::aid-fld11>3.0.co;2-7. doi:10.1002/1097-0363(20000615)33. [48] Chen JK, Beraun JE. A generalized smoothed particle hydrodynamics method for nonlinear dynamic problems. Comput Methods Appl Mech Eng 2000;190:225–39. doi:10.1016/s0045-7825(99)00422-3. [49] Grenier N, Le Touzé D, Colagrossi A, et al. Viscous bubbly flows simulation with an interface SPH model. Ocean Eng 2013;69:88–102. doi:10.1016/j.oceaneng.2013.05.010. [50] Hysing S, Turek S, Kuzmin D, et al. Quantitative benchmark computations of twodimensional bubble dynamics. Int J Numer Methods Fluids 2009;60(11):1259–88. doi:10.1002/fld.1934. [51] Ashgriz Nasser. Handbook of atomization and sprays. Springer; 2011. doi:10.1007/978-1-4419-7264-4. [52] Takashima T, Ito T, Shigeta M, et al. Simulation of liquid jet breakup process by three-dimensional incomressible SPH method. In: Proceedings of the seventh international conference on computer fluid dynamics (ICCFD7), Big Island, Hawaii; July 9–13; 2012.

147