Numerical simulation of Rayleigh-Taylor Instability with periodic boundary condition using MPS method

Progress in Nuclear Energy 109 (2018) 130–144

Contents lists available at ScienceDirect

Progress in Nuclear Energy journal homepage: www.elsevier.com/locate/pnucene

Numerical simulation of Rayleigh-Taylor Instability with periodic boundary condition using MPS method

T

Kailun Guoa,b, Ronghua Chena,b,∗, Yonglin Lia,b, Wenxi Tiana,b, Guanghui Sua,b, Suizheng Qiua,b a b

School of Nuclear Science and Technology, Xi'an Jiaotong University, No. 28, Xianning West road, Xi'an, 710049, China State Key Laboratory of Multiphase Flow in Power Engineering, Xi'an Jiaotong University, No. 28, Xianning West road, Xi'an, 710049, China

A R T I C LE I N FO

A B S T R A C T

Keywords: Multiphase moving particle semi-implicit (MMPS) method Rayleigh-Taylor instability(RTI) Multiphase flow Periodic boundary condition Initial velocity disturbance arrangement

The Multiphase Moving Particle Semi-implicit(MMPS) method is adopted to simulate the Rayleigh-Taylor Instability (RTI) process in an incompressible viscous two-phase immiscible fluid. To eliminate influences of the initial disturbance arrangement and boundary conditions in RTI simulations, the initial velocity distribution is deduced and the periodic boundary condition is developed in this paper. The RTI cases in this paper include those with single-mode disturbance and multimode disturbance. For the single-mode RTI process, cases with different initial disturbance wavelengths, different Atwood numbers and different surface tension coefficients are simulated. Results are quantitatively compared with analytical solutions in a linear growth stage and the simulation results fit well with the theoretical results. The long-term development of simulations is presented for investigating changes in the topology of rising bubbles and falling spikes in RTI, and the steady velocity conforms well with that calculate by theoretical solution. What's more, the multimode RTI processes are also simulated in this paper, and bubbles' competition and mergence process can be observed.

1. Introduction The Rayleigh-Taylor Instability (RTI) occurs at the interface of different fluids with density difference. If there is a pressure gradient from lighter fluid towards heavier fluid at the normal direction of the interface (equivalent to the gravity from heavier fluid towards lighter fluid), small disturbance will grow over time, this kind of interface instability is defined as the Rayleigh-Taylor Instability. When shock waves pass through the interface of two fluids with density difference, fluid at the interface will obtain an instantaneous acceleration, the interface instability will also happen, and this kind of interface instability is defined as the Rayleigh-Meshkov Instability (RMI). In nature and industries, the interface instability of stratified fluids are extremely important processes, RTI and RMI play important roles in some stages of stellar evolution, fragmentation of the gas film, crystallization of underground salt dome and the ignition stage of Inertial Confinement Fusion (ICF) (Meyer-ter-Vehn, 2009). In nuclear reactor severe accidents, RTI is an important influence factor (Kim and Corradini, 1988). In vapor explosion accident, the pressure vessel has broken, molten corium with high temperature, high density and high viscosity leaks into the cooling water tank. During the interaction of molten corium and coolant, the interface instability happens, these will lead to the fragmentation of the gas film near the ∗

molten corium, then the molten corium bursts forth and fragment. Heat transfer between molten corium and coolant increases during these processes which finally lead to the vapor explosion. The uncertainty of severe accidents increases due to the interface instability, thus the interface instability becomes an important influence factor in Fuel Coolant Interaction (FCI) process, and it is necessary to do some researches on these phenomena. Studies on RTI have been carried out for decades, especially in recent years. Due to the ICF ignition and FCI process have become hot topics of research, the study of RTI has entered a high growth stage, this phenomenon has obtained more and more attentions. Through experimental and theoretical studies, development of RTI can be broadly divided into the initial stage, linear stage and nonlinear stage (Kull, 1991). In the initial stage of RTI, the lighter fluid and heavier fluid invade into each other, but the amplitude of growing disturbance is small, there are only weak disturbances at the interface, so the influence of RTI in this stage can be ignored. With the development of the instability, more and more fluid exchange area, but the interface can be still observed clearly. The amplitude of the disturbance at the interface changes, and the growth of the amplitude follows the exponential law. This stage is considered as the linear growth stage of RTI. When the amplitude of the disturbance is larger than the 10% of the initial disturbance wavelength, the RTI is generally considered to be

Corresponding author. School of Nuclear Science and Technology, Xi'an Jiaotong University, No. 28, Xianning West road, Xi'an, 710049, China. E-mail address: [email protected] (R. Chen).

https://doi.org/10.1016/j.pnucene.2018.08.008 Received 11 December 2017; Received in revised form 22 May 2018; Accepted 6 August 2018 0149-1970/ © 2018 Elsevier Ltd. All rights reserved.

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Nomenclature

Greek letters

C

ϕ α κ λ μ ξ ρ σ

D ⇀ F G H L N ⎯⇀ N W ⇀ g h k l0 n p t u ⇀ u r ⇀ r x,y,z

dimensionless modified matrix/color function/Courant number/constant number dimension number force vector Gaussian kernel function hyperbolic kernel function/height, m length, m number of neighboring particles degree of irregularity of particles distribution generalized kernel function/width, m gravity acceleration, m/s2 amplitude, m wave number of the disturbance particle average distance, m particle number density pressure, Pa time, s particle velocity, m/s particle velocity vector distance between particles, m particle location vector position in each direction

physical quantity compressibility of fluid local curvature wavelength, m dynamic viscosity, Pa•s dynamic parameter density, kg/m3 surface tension coefficient, N/m

Superscript/subscript 0 1 2 i j S ave min *

initial condition fluid phase 1/constant in Eq. (32) fluid phase 2/constant in Eq. (33) No. of target particles No. of neighboring particles surface tension average value minimum value temporal value

Front Tracking, VOF and LBM methods are the most common approaches in RTI simulations. In recent years, the DNS is also adopted in RTI simulation with the development of computer technology (Reckinger et al., 2016). There are also a series numerical research about RTI using meshless method, which are mainly refer to those particle methods such as Smoothed Particle Hydrodynamics (SPH) and Moving Particle Semi-implicit method (MPS). But the RTI simulations carried out by these particle methods are mainly adopted to validate different multiphase flow algorithms of each researcher. Simulation results in these researches lack of rigorous theoretical analysis and always show poor agreement with the existing RTI theory. For example, Jie Liu et al. (2005) used MPS-FVM method to simulate the RTI with single-mode disturbance, he compared his simulation results with the theoretical solution in linear growth stage. He concluded that when the surface tension force was not involved in the simulation, simulation results were far from the theoretical solution, but when the surface tension force was considered, the situation can be improved. In Shakibaeinia's paper (Shakibaeinia and Jin, 2012), he adopted the explicit multiphase algorithm similar to that in SPH to calculate the fluid field pressure, and a RTI simulation with upwards disturbance was simulated. The results were compared with those calculated by VOF, but deviations were shown between these two methods. Jeong used the pressure gradient model based on sum of neighbor particles' pressure (Jeong et al., 2013) to simulate the same RTI case as Liu's in reference (Liu et al., 2005), but the interface deformation in his paper was significantly slower than Liu's results. There is a widespread problem in all of the studies above: the initial disturbances in previous simulations of the RTI processes using particle methods are all displacement disturbances formed by particles along interfaces. The resolution of the interface can't be high enough to describe irregular shape (e.g. a Cosine wave shape) due to the characteristic of the Lagrange particle method. This leads to a large different between the results of particle methods and those of grid methods or the theoretical solution, especially in the early stage of RTI. The shape of the disturbance in the linear stage of RTI is closer to the square wave rather than the Cosine wave (Liu et al., 2005; Shakibaeinia and Jin, 2012). What's more, the boundary conditions in the literature above are all set as the no-slip boundary, but in many other RTI numerical researches with grid methods, boundary conditions are often set as the periodic boundary and free-slip

the non-linear growth stage (Jacobs and Catton, 2006). The non-linear stage can be further divided into steady growth stage, irregular nonlinear growth stage and turbulent mixing stage (Dimonte, 2000). In the steady growth stage, the initial disturbances at the interface rapidly develop into bubbles (lighter fluids flow into heavier fluids) and spikes (heavier fluids invade into lighter fluids) which both develop with the terminal velocity. Another interface instability, Kelvin-Helmholtz Instability (KHI) will happen during this process, thus the KHI can be regarded as a kind of phenomenon that along with the RTI. The top of bubbles and spikes will roll up when the KHI appears, and the shape of them transforms like “mushrooms”. In the multimode disturbance RTI, bubbles will compete and merge with each other, “mushrooms” will break, the interface becomes complicated and hardly to predict. The RTI finally grows into the turbulent mixing stage. Experimental researches laid the foundation of the field of interface instability. Taylor first observed the single bubble's rising at a constant speed in a RTI experiment (Taylor, 1950). In 1973, Ratafia observed the formation of bubbles and spikes, he first found the rolling up of spikes' top (formation of “mushrooms”) due to the KHI (Ratafia, 1973). Since 1980s, many elaborate and creative researches has been carried out with the demand of the high-tech field such as ICF. Ken Read adopted a vertical acceleration vessel to study the variation of the interface during RTI in 1984 (Read, 1992). Dimonte & Schneider used straight line engine to accelerate the interface and they got the characteristic of the RTI under the complicate acceleration (Dimonte et al., 2007). Dalziel et al. (Dalziel, 2000) isolated two fluids with plate, and the process of RTI under the gravity was studied. In the field of nuclear reactor severe accidents field, lots of researches on FCI have also been carried out since 1990s. Dinh et al. (1999) carried out some researches on the jet breaking in FCI process. In recent years, Hyo Heo et al. (2015) poured molten Wood metal into water to study different influential factor in FCI. The liquid-liquid injection experiment was carried out by Shimpei Satio (Saito et al., 2015) to study the influences of density ratio and viscosity ratio on the FCI process. All these researches have shown that the interface instability is an important impact factor in jet breaking behavior of FCI process, RTI can be observed in these experiments. In numerical simulations of RTI, Youngs (1984), Glimm (Glimm and Li, 1988) and Li Xiaolin (Li, 1993) used Interface Tracking, MAC and Level Set method to simulate the RTI in 2D and 3D, respectively. The 131

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

part of Source term (ECS) is adopted into the PPE:

boundary. The no-slip boundary in particle method always leads to an incomplete bubble or spike near the calculation boundary. Shuai Zhang (Zhang et al., 2004) firstly introduced the periodic boundary condition in MPS method to simulate the RTI process, but his MPS code included the artificial buoyancy model which is not widely adopted by other researchers and his simulation results were short of quantitative analysis. As comparison, the RTI numerical studies of the SPH method are very similar to those of the MPS method. There has been little research on the long-term development of RTI, most of them are adopted to validate researchers' own multiphase flow algorithms (Cummins and Rudman, 1999; Tartakovsky and Meakin, 2005; Hu and Adams, 2007; Shao and Lo, 2003; Grenier et al., 2009). Besides that, there are also problems about boundary conditions and arrangement of initial disturbance in SPH as mentioned above. So there is little research about the RTI phenomenon using meshless methods. Especially, the MPS method is more immature due to its appearance is later than the SPH. To the author's knowledge, there is no detailed simulation and quantitative analyses on the long-term development RTI process using the MPS method. In this paper, the Multiphase MPS method (MMPS) is chosen as the research tool, limitations of the boundary conditions and the shape of the interface in previously meshless methods' simulations will be overcome. Each stage of RTI will be simulated and analyzed. What's more, bubble's competition and mergence in multimode disturbance RTI will be also simulated.

Δt 2 k + 1 ∇P ρ

j

μi + μ j

ni nj2

9 9r 2 exp ⎛⎜− 2 ⎞⎟ πre2 ⎝ re ⎠

∂2Wij D − 1 ∂Wij ⎞ ⎞ ⎛ ∂φij ∂Wij + φij ⎛⎜ 2 + ⎟ ∂ rij ∂rij ⎠ ⎟ r ∂ r ⎝ ij ⎠ ⎝ ij ij

∑ ⎜ ∂r j≠i

(6)

e ⎧ r − 1 (r < re ) ⎨ 0 (r ≥ re ) ⎩

(7)

The left hand side of Eq. (5) transfers into:

∇2 P k + 1 i =

1 n0

∑ ⎛⎜ j≠i

3Pij re ⎞ ⎟

3 ⎝ rij ⎠

(8)

Then we can get the complete discretized scheme of PPE:

1 n0

∑ j≠i

1 ⎛ 3Pij re ⎞ nk − n0 nk − n0 ⎞ = ∇⋅⇀ u∗+ u k ) + ∇⋅⇀ uk ⎛ (∇⋅⇀ ⎜ 3 ⎟ 0 0 ρi ⎝ rij ⎠ n ⎠ ⎝ n 1 k+1 + α 2 Pi Δt ⎜

⎟

(9)

Eq. (9) can be solved by the Incomplete Choleskey Conjugate Gradient (ICCG) method, and the pressure of whole fluid flied can be got. The gradient model as Eq. (10) developed by Duan (Duan et al., 2017) has been validated accurate and stable enough in most multiphase flow simulations, so Eq. (10) is also adopted in this paper: 1 ∇P ρ

+

d n0

=

∑j ≠ i

d n0

∑j ≠ i

⇀ ⇀ 2 ⎧ (Pj − Pi)( rj − ri ) Cij H ( ⇀ ρi + ρj ⎨ rj − ⇀ ri 2

⎩

⎧ ⎜⎛Pi − Pi′,min ⎟⎞ (⇀ rj − ⇀ ri ) 1 ⎪⎝ ⎠ Cij H ( ⇀ ρi ⎨ rj − ⇀ ri 2 ⎪ ⎩

⇀ rj − ⇀ ri ) ⎫ ⎬ ⎭

⎫ ⎪ ⇀ rj − ⇀ ri ) ⎬ ⎪ ⎭

(10)

where Cij is a dimensionless modified matrix introduced by Khayyer & Gotoh which can make the gradient model more accurate (Khayyer and Gotoh, 2011). The CCSF model (Duan et al., 2015) as follows is also adopted in this paper:

⇀ FS ρ

= i

2σκi ∇C ρ1 + ρ2

(11)

where σ is the surface tension coefficient, κi is the local curvature at target particle i, C is the color function, ρ1 and ρ2 are constant density for different phase. The detailed solution process can be refered in Dr. Duan's research (Duan et al., 2015). What's more, the Dynamic Displacement Modification (DDM) technique (Zainali, 2013) is adopted to make the calculation more stable:

uij Gij ×⇀ (3)

where re is the radius of the interaction zone, ⇀ uij = ⇀ uj − ⇀ ui , n is the Particle Number Density (PND), Gij is the Gaussian kernel function (Duan and Chen, 2013):

G (r ) =

(5)

r

H (r ) =

D⇀ u ⇀ u + ρ⇀ g + FS = −∇p + μ∇2→ (2) Dt ⇀ ⇀ ⇀ where u , g , p , ρ and FS are velocity, gravity, pressure, density and surface tension, respectively. In MPS, the governing equations are solved by particle interaction models such as gradient model and Laplacian model. The momentum conservation equation is dived into viscosity, pressure gradient, gravity and surface tension forces terms, pressure is solved by PPE through implicit method while rest parts of Eq. (2) are addressed explicitly. If the MPS method wants to be extended into the MMPS method, the first problem should be solved is that the different physical properties between two fluids. The sharply density/viscosity difference makes the acceleration of particles at interface discontinuous, which will always lead to unstable simulation results. So the multi-viscosity model as follows are adopted to make properties continuous across the interface (Hu, 2006):

2μi μj (ni2 + nj2)

1 k+1 Pi Δt 2

where Wij is the generalized kernel function. By considering the classical hyperbolic function as the kernel function:

ρ

∑

1 n0

∇2 φ i =

(1)

18 re2

⎟

The first order accuracy Laplace model introduced by Khayyar (Khayyer and Gotoh, 2010) is used to discrete the left hand side of Eq. (5):

The MMPS method which is developed based on Koshizuka's classical MPS method (Koshizuka and Oka, 1996) is initially introduced by Dr. Duan (Duan et al., 2017). The main improvements of MMPS method are primarily about the Pressure Poisson Equation (PPE), pressure gradient model and viscosity model. In addition, Duan also developed Contoured Continuum Surface Force model (CCSF) (Duan et al., 2015) for the MMPS method. The governing equations are showed as follows:

μ⋅∇2 ⇀ u =

⎜

i

+α

2. MMPS method

∇⋅⇀ u =0

nk − n0 nk − n0 ⎞ (∇⋅⇀ u k ) + ∇⋅⇀ uk ⎛ 0 0 n ⎠ ⎝ n

u∗− = ∇⋅⇀

Δri = ξ ∑ j≠i

⇀ ri − ⇀ rj rij3

(riave )2u max Δt (12)

where riave is the average distance between neighboring particles j and target particle i which is calculated as follows:

(4)

riave =

∑ (rij/Ni) j

where r is the distance between particle j and i. The Error Compensating 132

(13)

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

ξ in Eq. (12) is the displacement modification parameter which is calculated as follows:

⎯⇀ ξ = 0.01 + 0.2 × N

when a Lagrange method is adopted to describe an irregular initial interface, the resolution won't be very high, which resulting in the shape of interface deviates Cosine wave farther. This problem exists in a large number of articles of simulation RTI with particle methods, but does not seem to cause enough attention of researchers. It is because they mainly carried out RTI calculation examples to detect the stability of their algorithms. These simulation results were compared with other researchers', but the influence of initial interface on the development of RTI were less considered. M. S. Shadloo found the influence of initial interface arrangement on RTI development in his SPH article (Shadloo et al., 2013), so he proposed the cubic arrangement, staggered arrangement, radially-centered arrangement and radially-off-centered arrangement to study the influence of initial displacement disturbance, but interfaces formed by these ways are still irregular. Although in some articles, the consequences of convergence tests show that the particle configuration has reached the optimum, the simulation results are always in error with the theoretical solution or the mesh method. In fact, the irregular initial displacement disturbance interface is equivalent that there are many tiny local disturbances along the disturbance interface, which cause deviations in development of RTI, especially in the initial linear growth stage.

(14)

⎯⇀ where N stands for the degree of irregularity of particles distribution (Tamai and Koshizuka, 2014):

1 ⎯⇀ N = ni

∑

⇀ rj − ⇀ ri rij

j≠i

Hij (15)

In Ref. (Zainali, 2013), the ξ is equal to a constant number which can't guarantee it is suitable for the whole calculation domain or the whole calculation process. So in this paper, the degree of irregularity of particles distribution shown in Eq. (15) is adopted to modify the parameter ξ. Eq. (15) means that more irregularly particles distribute, ⎯⇀ the larger N and ξ will be, and a dynamic adjustment parameter ξ will be generated according to the degree of irregularity of particles distribution in each time step, which can make sure the adjustments in this paper will be always more proper than those in Ref. (Zainali, 2013). The density smoothing scheme is also involved in this paper. Here we adopt the First order accuracy Density Smoothing (FDS) developed by Khayyer et al. (Khayyer and Gotoh, 2011)

ρi =

1 ∑ Wij

∑ ⎛⎜ρj − ⎝

∂ρ ⎞ ∂ρi x ij − i yij ⎟ Wij ∂yij ∂x ij ⎠

3.2. Initial velocity disturbance (16)

The initial velocity disturbance refers to an initial Cosine velocity distribution of the flow field. This kind of initial arrangement can be seen in some simulations using grid methods, but in particle methods, the RTI simulation arranged in the initial velocity disturbance has not been found. When the initial velocity is assigned to the flow field, if the interface is simply assigned v0=Bcos(kx), the initial velocity field will not guaranteed to obey the continuity equation (Eq. (1)). Therefore, the following initial velocity field is considered in this paper (initial disturbance velocity is up along the Y axis at x = 0): 2D:

When Eq. (4) is adopted as the kernel function in Eq. (16), the FDS finally transfer as:

1 ρi = ∑j Gij

∑ j

⎛ ⎛ ⎞⎞ ⎜ 18rij ⎜ ∑j rij ρj Gij × ∑j Gij − ∑j ρj Gij × ∑j rij Gij ⎟ ⎟ ⎜ρj + 2 ⎜ ⎟ ⎟ Gij 2 re ⎜ ⎛⎜∑ G ⎞⎟ ⎜ ⎟⎟ j ij ⎜ ⎟⎟ ⎜ ⎝ ⎠ ⎝ ⎠⎠ ⎝ (17)

The calculation process of the MMPS method in this paper is summarized as follows: The governing equations are shown as Eqs. (1) and (2), the viscosity term is firstly explicitly solved by the multi-viscosity model as Eq. (3). The surface tension is also solved in an explicit way. Then we can get a temporal velocity ⇀ u * for each particle. The particles' pressure is solved ' u can by the PPE as Eq. (5) in an implicit way, then the correct velocity ⇀ be got. Finally, particles' real velocity in next time step can be calculated as ⇀ u =⇀ u*+⇀ u ′.

−g y ⋅sign(y ) ⎧u 0 (x , y ) = B sin(kx )⋅e ⎪ Bk v0 (x , y ) = g cos(kx )⋅e−g y ⎨ ⎪ w0 (x , y ) = 0 ⎩

3D: −g y ⋅sign (y ) ⎧u 0 (x , y, z ) = B sin(kx )⋅e ⎪ Bk v0 (x , y, z ) = g (cos(kx ) + cos(kz ))⋅e−g y ⎨ ⎪ w0 (x , y, z ) = B cos(kz )⋅e−g y ⋅sign (y ) ⎩

3. Initial condition

(20)

The initial velocity field obtained by above equations can guarantee the conservation of the continuity equation in Eq. (1).

3.1. Initial displacement disturbance In the numerical simulation of RTI, the initial condition is mainly the arrangement way of initial disturbance at the interface, and most of the studies adopt the initial disturbance shape as Cosine wave. There are two kinds of disturbance arrangement: displacement disturbance and velocity disturbance. When displacement disturbance is adopted, the shape of Cosine wave should be initially arranged at the interface. Consider the layout of the coordinates as shown in Fig. 1, take the center of calculation domain as the origin, and the central axis as Y axis. The width of calculation area is same as the wavelength of initial disturbance, then the Cosine arrangement of initial interface can be written as mathematical expression:

h (x ) = h0cos (kx )

(19)

(18)

where k=2πλ/2 is the wave number of the disturbance, and h0 is the amplitude of the initial disturbance. Take h0=0.02λ, the interface arranged by the above formula is shown below in Fig. 2 (a), t*=t(g/H)0.5 is dimensionless time scale adopted in this paper. As mentioned above,

Fig. 1. Diagram of calculation domain. 133

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 2. Interface shape under different initial arrangement.

(a) and (b) are equivalent, except the initial disturbance. Results show that the development of RTI has been greatly different under the two work conditions. Because of the low resolution of the interface, the RTI of the initial displacement disturbance has been completely divorced from the development trend of the cosine wave. The bubbles formed by the lighter fluid are also irregular, and the cusp of the heavier fluid is more likely to form the square wave. But the initial velocity disturbance does not have this problem, and always develops in a smooth sinusoidal form. It should be noted, however, that the difference caused by the arrangement of such initial disturbance is closely related to the initial wavelength. Fig. 5 shows comparison of the calculation results in different initial interface arrangements when λ = 0.8m. It can be concluded that although the initial displacement disturbance still cause deviations, the difference between calculation results under these two kinds of initial disturbances is smaller. As shown in Fig. 6, when the initial wavelength is further reduced to 0.5m, the results of the initial displacement disturbance and the initial velocity disturbance is very similar, but there are still small wavy local disturbances along the initial displacement disturbance interface, it is more likely to form a sharp spike roll. Fig. 7 shows the comparison between initial displacement disturbance interface under different wavelengths. It can be concluded that when the wavelength is smaller, local disturbances along the interface (marked by white circles) are also reduced correspondingly, so their influences on the RTI development process are smaller. Through the above discussion of the initial interface arrangement, it can be concluded that the initial velocity arrangement shown in Eqs. (19) and (20) should be the first choice for the RTI simulation with the initial disturbance as Cosine wave. When the wavelength is large, calculation results of the displacement disturbance and velocity disturbance have large deviation. While when the wavelength is small, the

Fig. 3. Initial velocity field under initial velocity disturbance.

Taking two-dimensional disturbance as example, set the wavelength λ=1.0m, the initial disturbance velocity direction is downward at x = 0, and the coefficient B=0.1λ. The density of heavier fluid ρ1 = 3 kg/m3, lighter fluid's density ρ2 = 1 kg/m3, the dynamic viscosity for two fluids is μ1 = 0.03 Pa s and μ2 = 0.01 Pa s, respectively. The average particle distance l0 is set as 0.01m, and the height of the calculation domain is 2.0m, so the total particle number is 20000 for this case. The initial velocity field distribution nephogram can be obtained as shown in Fig. 3. Fig. 4 shows the different calculation results under different initial disturbance arrangement, and all the calculation parameters in Fig. 4

Fig. 4. Comparison of the calculation results in different initial interface arrangements (λ = 1.0m) (a) initial displacement disturbance (b) initial velocity disturbance. 134

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 5. Comparison of the calculation results in different interface arrangements (λ = 0.8m). (a) initial displacement disturbance; (b) Initial velocity disturbance.

Fig. 6. Comparison of the calculation results in different initial interface arrangements (λ = 0.5m). (a) initial displacement disturbance (b) initial velocity disturbance.

deviation will be weak. Through several calculations, it is found that the initial displacement disturbance won't cause much deviations when the height of the initial amplitude h0 is not larger than twice of the average particle distance l0. So in the paper, the criteria when the initial displacement disturbance can be adopted to simulate the RTI process is shown as follows:

Fig. 8. Diagram of periodic boundary condition.

4. Boundary condition

λ ≤ 100l 0

(21) The setting of boundary conditions has great influence on simulation results of RTI. Many literature of RTI numerical simulation, especially those using meshless method, tends to adopt the no-slip condition around the boundary. While it is difficult to obtain the complete bubble

The displacement disturbance interface arrangement is only suggested when the initial disturbance wavelength meets the above criteria.

Fig. 7. Initial Displacement Disturbance Interface. (a)λ = 1.0m; (b) λ = 0.8m; (c) λ = 0.5m. 135

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 9. Diagram of free slip boundary. Table 1 Summary of single-mode disturbance operating condition. Para Case

ρ1/(kg/ m3)

ρ2/(kg/ m3)

A

μ1/Pa·s

μ2/Pa·s

Φ

g/(m/ s2)

λ/m

1 2 3 4 5 6 7 8 9

3 3 11 13 9 3 3 3 3

1 1 9 7 1 1 1 1 1

0.5 0.5 0.1 0.3 0.8 0.5 0.5 0.5 0.5

0.03 0.03 0.11 0.13 0.09 0.03 0.03 0.03 0.03

0.01 0.01 0.09 0.07 0.01 0.01 0.01 0.01 0.01

0.0 0.0 0.0 0.0 0.0 0.125 0.25 0.50 1.0

10 10 10 10 10 10 10 10 10

0.5 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0

Fig. 11. Comparison of bubble rising velocity in case 1 under different particle Configurations.

method firstly (Zhang et al., 2004), but didn't describe it in detail. For the calculation example with random disturbance (usually not arrange the initial disturbance), if the boundaries are set as the no-slip boundary condition, RTI phenomenon will not happen (Zhang et al., 2004), unless the upper side boundary is set as the free surface boundary. This is not consistent with the actual physical phenomenon. In this paper, the periodic boundary condition of the particle method is developed with the idea of dummy buckets and dummy particles. The diagram of periodic boundary condition is shown in Fig. 8. The both sides of calculation domain are set as the periodic boundary conditions, the calculation domain is divided into several square sub domains by taking particle interact radius as side length, sub domains of each column are numbered as 1 to N from left to right. If particles are located in the sub domains of column iX, the particle can only react with the particles in column iX+1 on the right side and those in column iX-1 on the left side. When the column 1 on the most left side is iX, there is no bucket on its left side. In order to achieve periodic boundary condition, the column N on the most right side need to be spliced on the left side of column 1 as dummy buckets, the distance between particle in dummy bucket and actual particle is the width of the calculation domain. In addition, as shown in green areas in Fig. 8, the bucket width can't always guarantee buckets in column N is full with particles in the initial partition of buckets. At this time, column N-1 also needs to be spliced into the left side of the calculation area as the dummy buckets of column iX-2. Thus, the periodic boundary condition of the left calculation domain is implemented. For the right side periodic boundary condition, the principle is the same, but buckets in column 2 are not needed to be spliced, because the code ensures that buckets in column 1 are always full of particles. The advantage is the small calculation amount by setting the periodic boundary condition like this. Particles are not needed to be searched in whole calculation domain. It only needs to identify two columns of buckets near calculation boundaries at most, then particles in the buckets are implemented the above process. Moreover, this method has high calculation accuracy, there won't be a miscarriage of justice for boundary particles.

Fig. 10. Comparison of bubble rising height in Case1 under different particle Configurations.

or spike shape near the boundary. Therefore, in this paper, the periodic boundary condition is adopted on the left and right boundaries of the calculation area, and the free slip boundary condition for the upper and lower boundaries. 4.1. Periodic boundary conditions

4.2. Free slip boundary condition In numerical simulations of RTI taken in FVM or LBM, the setting of periodic boundaries is quite common. But for the MPS method, it is hardly to find the calculation example contains periodic boundary. Shuai Zhang implemented the periodic boundary condition in MPS

The achievement of free slip boundary condition is relatively simple, the diagram is as shown in Fig. 9. As long as neighbor wall particles in the interaction zone of target particle i are guaranteed to 136

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 12. Simulation results of Case1 at different times.

obvious feature is that the growth of disturbance amplitude is satisfying the exponential relation, and can be well described by a mathematical model. If the initial disturbance at the interface is arranged as the initial conditions mentioned in Section 4.2, Eq. (22) can be obtained by linear approximation:

d2h (t ) + n2h (t ) = 0 dt 2

(22)

where h(t) is the disturbance amplitude of the time t, and n is the linear growth rate of the RTI in the linear growth stage. The classical linear growth rate is calculated as follows (Danielson, 1962):

n=

k 2σ ⎤ kg ⎡A− ⎥ ⎢ g ( ρ 1 + ρ2 ) ⎦ ⎣

(23)

where A=(ρ1-ρ2)/(ρ1+ρ) is Atwood number, which is the dimensionless parameter to describe the density ratio of fluids. However, this linear growth rate has been proved by many researchers to be inaccurate, because it does not consider the influence of fluid viscosity, the growth rate equation with viscous influence is (Mikaelian, 1993):

k 2σ ⎤=0 n2 + 2νk 2n − kg ⎡A− ⎥ ⎢ g ( ρ 1 + ρ2 ) ⎦ ⎣

(24)

where ν=(μ1+μ2)/(ρ1+ρ2) represents the average kinematic viscosity fluids. The exact expressions for the linear growth rate can be obtained:

Fig. 13. Pressure nephogram and isobars at different time of Case1.

n = −νk 2 +

have same velocity in value and direction as target particle i. It can ensure there is no the relative speed between the boundary and fluid near the boundary.

k 2σ ⎤ + (νk 2 ) 2 kg ⎡A− ⎢ g (ρ1 + ρ2 ) ⎥ ⎦ ⎣

(25)

The differential equation as Eq. (22) can be solved by considering initial conditions as the displacement disturbance (initial velocity is 0 and initial amplitude is h0):

5. Conclusions and analysis

dh (t )

5.1. Theoretical analysis of RTI

⎧ dt t = 0 = 0 ⎨ h (t ) t =0 = h0 ⎩

5.1.1. Linear growth stage of RTI In the linear growth stage of RTI, as mentioned above, its most

The relationship between amplitude and time can be obtained with displacement as the initial disturbance at the linear growth stage of the 137

(26)

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 14. Simulation results of Case2 at different times.

scholars often use a simple formula (Liu et al., 2005; Shakibaeinia and Jin, 2012) in the form of h(t)∝exp(nt), while study in this paper found that the simple formula is not accurate, and can not exactly describe the amplitude change situation in linear growth stage, concrete results will be shown in the paper later. The RTI under the effect of surface tension is also studied in this paper. Assuming that the width of calculation domain is exactly equal to the initial disturbance wavelength, when the linear growth rate as shown in Eq. (25) is 0, critical surface tension coefficient σc can be obtained:

σc =

RTI:

h0 [exp (nt ) + exp (−nt )] 2

(30)

The critical surface tension coefficient can make the linear growth rate be 0, which means RTI will not happen under the limitation of surface tension. When the surface tension coefficient is 0, the linear growth rate is the biggest, and the development of RTI is the most severe. Therefore, it can be inferred that the surface tension is the force that maintains the interface shape, the development of RTI is slower when the surface tension coefficient is lager. In this paper, the following dimensionless numbers are adopted to represent the magnitude of surface tension:

Fig. 15. Velocity Distribution of Case2 at t* = 2.57.

h (t ) =

Ag (ρ1 + ρ2 ) k2

Φ= (27)

σ σc

(31)

In subsequent numerical simulations, the effect of surface tension on the development of RTI will be investigated by changing the different Φ.

If the velocity disturbance is adopted as initial condition, which means the initial velocity is v0(0,0), and the initial amplitude is 0: dh (t )

⎧ dt t = 0 = Bk / g ⎨ h (t ) t=0 = 0 ⎩

5.1.2. The nonlinear development stage of RTI The nonlinear growth stage of RTI is one of key research objects in this paper. Because the turbulent mixing stage is a very intense stage which is full of uncertainty, the main research work is often concentrated in the steady stage of nonlinear growth stage. This paper will also focus on this stage. The main research object is the velocity of bubbles and spikes at the steady stage. In the nonlinear growth stage, the development process of RTI is difficult to be described by a theoretical solution. Therefore, the main research approach is to determine the empirical relation formula of

(28)

The amplitude variation relation formula with velocity as the initial disturbance at the linear development stage of the RTI can be obtained:

h (t ) =

Bk / g [exp (nt ) − exp (−nt )] 2

(29)

Eq. (29) is also the relation of amplitude in linear growth stage adopted in the theoretical analysis in this paper. In many studies, 138

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 17. Comparison of bubbles' amplitude and velocity under different A.

Fig. 16. Quantitative analysis of the RTI development of Case1 and Case2. Table 2 Comparison of the value of C1 and analytic value under different conditions.

C1,b C1,s

Case1

Case2

Analytical

Error1

Error2

0.27 0.46

0.27 0.47

0.27 0.46

0.00% 0.00%

0.00% 2.17%

instability development through experiments. The most classical one is the bubble rising velocity relation at the steady stage proposed by Taylor (1950):

vb,terminal = C1 gλ / 2

(32)

By the method of experiment, he obtained that bubble rising velocity at the steady stage was proportional to (gλ)0.5, where C1 is a constant number changing along with A. This conclusion was also demonstrated by subsequent most studies (He et al., 1999; Abarzhi, 1998). Subsequently, further studies have shown that in the development of RTI, not only the rising of bubbles has a steady stage, but also the spike dropping has this process. When the buoyancy force and resistance of the fluid are balanced, the steady velocity of the bubble and spike can be calculated by the following equation (Reckinger et al.,

Fig. 18. Comparison of spikes' amplitude and velocity under different A.

139

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

2016):

Table 3 Steady velocity constants under different A.

C1,b Analytical Error C1,s Analytical Error

1/2

A = 0.1

A = 0.3

A = 0.5

A = 0.8

0.14 0.14 0.00% 0.15 0.15 0.00%

0.22 0.22 0.00% 0.30 0.30 0.00%

0.27 0.27 0.00% 0.47 0.46 2.17%

0.31 0.31 0.00% 0.92 0.92 0.00%

2A gλ ⎤ vb/s = ⎡ ⎢ (1 ⎦ ⎣ ± A) C2 ⎥

(33)

For two-dimensional cases, C2=6π, for three-dimensional cases, C2=2π. It should be noted that the results obtained by Eq. (33) conform well with the experimental results obtained in Eq. (32) by Taylor, and the above two equations both determine the relationship vb∝(gλ)0.5. In this paper, steady velocity of the bubble in the nonlinear stage of R-T instability will be studied based on the theoretical solution of bubble velocity and spike velocity in Eq. (33). 5.2. Numerical simulation of RTI with single-mode disturbance 5.2.1. Calculation parameters In single-mode disturbance cases, the development of RTI in initial disturbances with different wavelength, different Atwood numbers, and different dimensionless surface tension coefficients Φ are simulated in this paper successively. Calculation domain is shown in Fig. 1, for both sides of the domain, periodic boundary condition introduced in section 4.1 is adopted, the up and down of calculation domain adopt free slip boundary condition. Initial disturbance wavelength is exactly equal to the width of the calculation domain, the initial interface disturbance adopts the initial velocity arrangement, the parameter B is set as 0.1λ in Eq. (19). Table 1 summarizes the RTI operating conditions of all singlemode disturbance simulated in Section 5.2 and their corresponding parameter settings. As shown in Table 1, RTI simulation of single-mode disturbance includes 9 kinds of different conditions, where Case1 and Case2 are mainly adopted to study the influence of different initial disturbance wavelength on RTI development process, at the same time, Case1 is also used as the convergence analysis of RTI numerical simulation. Case2Case5 mainly focused on the influence of Atwood number on the development of RTI, and the Atwood number covered almost all the conditions from 0.1 to 0.8. Case2 and Case6-Case9 are mainly adopted to examine the influence of surface tension on the development of RTI. The actual surface tension coefficients are set as 1/8, 1/4, 1/2 and 1 time of the critical surface tension coefficients respectively.

Fig. 19. Comparison of bubbles' amplitude and velocity under different surface tension coefficients.

5.2.2. Convergence test Before the numerical simulation is carried out, the convergence analysis of the calculation results should be adopted to determine the reliability of the results. For numerical simulations, it is usually necessary to guarantee the independence of the space step and the time step of the result. In the MMPS code adopted in this paper, time step is controlled by following equations (Duan et al., 2017):

Δt ≤ 0.5⋅

Δt ≤

Δt ≤

ρmin l 03 2πσ

(34)

re2 18dνmax

(35)

C0 l 0 u max

(36)

Eq. (34) is the time step controlled by the surface tension, and Eq. (35) is the time step controlled by the viscosity accuracy condition, Eq. (36) is the CFL condition. The time step is determined by the minimum time step in the above three equations. In the calculation, we usually use the program to automatically assign the time step, so the verification of the time step independence is not involved in this section. In this section, Case1 in Table 1 is adopted as a convergence case, change the different particle configurations, and the results obtained under different particle sizes are shown in Figs. 10 and 11. Fig. 10 shows the contrasting of rising height of lighter fluid in RTI with different particle sizes, it can be concluded that, whether λ/l0 is 50, 100 or

Fig. 20. Comparison of spikes' amplitude and velocity under different surface tension coefficients.

140

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 21. Rti process with dual mode disturbance.

results of Case1 at different time are shown in Fig. 12. Simulation results in this section are normalized by SPLASH software (Price, 2007) to optimize its display effect. As shown in Fig. 12, it is possible to observe the variation of the interface disturbance at different times. The linear growth of the harmonic, the formation and deformation of bubbles and spikes can be observed. Fig. 12 also shows the streamline distribution in the process of RTI development. At the linear growth stage of the interface, it appeared two vortex which are symmetrical with Y axis, in the whole linear growth stage, the center of the vortex are at y = 0, but its action scope is increasing. With the continuous development of RTI process, when the spike had rolled, vortex will move with the emergence of spikes, the vortex center is always at the tail of RTI spikes. With the further development of the interface instability, new vortex forms at bubble tail. Fig. 13 shows the dimensionless pressure distribution and the isobar distribution of Case1 in several typical moments, it is the time t* = 0.43, t* = 1.09, t* = 2.46 respectively from left to right, where p is the average pressure of the whole flow field, p0 = 1Pa is the reference pressure used to make pressure dimensionless. It can be concluded from isobars in Fig. 13 that pressure changes rapidly in heavier fluid, while in the lighter fluid, pressure changes slowly. With the development of the interface instability, two fluids invaded into each other, pressure began to fluctuate, and the pressure fluctuation's boundary is the interface between two fluid. In the whole RTI process, the pressure difference between the bottom and the top of the calculation domain is always maintained at about 40, which conformed well with the theoretical pressure difference Δp=(ρ1+ρ2)gH/2/ P0 = 40, and the correctness of the pressure calculation in the RTI process is proved. Fig. 14 shows the RTI simulation results of Case2 with an initial wavelength of 1m. Compared with Case1, only the wavelength of initial disturbance is changed. In order to facilitate the display, snapshots in Fig. 14 is displayed as the same scale as those shown in Fig. 12. It can be observed that the RTI development process of Case2 is faster than that of Case1. When t* = 3.42, spike arrives at the bottom of the calculation domain, which is because that, the longer the initial interface disturbance wavelength, the more rapidly instability phenomenon develop. However, the development trend of RTI with different

Fig. 22. Amplitude change of bubble in RTI process with dual mode disturbance.

200, the rising height changes of lighter fluid conformed very well. Fig. 11 shows the variation of the bubble rising velocity, its change trend is also conforming under different particle configurations. Velocity has some fluctuations in the process of bubble steady rising, but its average value is conforming. By comparing the calculation results for the different particle configurations of Case1, it can be concluded that when the particle configuration λ/l0 is 50, the calculation has been converged. With this configuration, the total number of particles is less, which can speed up the calculation, so in the subsequent numerical simulations, we calculate by the particle configuration of λ/l0=50. 5.2.3. RTI simulation under different initial disturbance wavelengths In this section, the influence of initial disturbance wavelength on RTI is analyzed by simulation of Case1 and Case2 in Table 1. Simulation 141

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

Fig. 23. Rti process with seven mode disturbance.

growth stage earlier. In addition, the green dotted line and the pink dotted line represent the fitting curve of the rising of the bubble and the spike drop at the steady stage. In the RTI simulation of two cases, whether in the development of bubbles or spikes, a uniform rising stage can be both observed. The slope of the fitting curve is the steady speed of bubble or spike. Before the development of bubble or spike in steady velocity, there is a very short period of amplitude deviate the exponential development law, but has not yet reached the stage of steady velocity, this stage is called the weakly nonlinear growth stage. According to this, the nonlinear development stage in Fig. 16 can be further divided into ② and ③ two areas. In this paper, the study of steady velocity is mainly concentrated in the ③ area. Velocities of bubbles and spikes can be obtained through the slope of fitting curves: vb1 = 0.42 m/s, vs1 = 0.73 m/s, vb2 = 0.60 m/s, vs2 = 1.04 m/s. If the formula of steady velocity is unified into the form of Eq. (32), the C1 obtained by the simulation can be compared with that calculated by Eq. (33). The comparison results are shown in Table 2. It can be found in Table 2 that the initial disturbance wavelength of Case1 and Case2 is different, but the corresponding constants obtained are the same due to the same A of the two cases. This also conform that in Eq. (33), steady velocity constant is related to A. Moreover, in the case of A = 0.5, the steady velocity constant and analytical solution calculated in the two cases is very consistent. It is proved that the result of MMPS is correct in the steady growth stage of simulating single-mode disturbance RTI.

wavelengths is similar, all can observe the process such as interface instability development process, the formation and rolling process of spike and bubble, etc. Fig. 15 shows the velocity distribution nephogram of Case2 at the moment t* = 2.57, the velocity distribution nephogram of X axis direction shows the rolling up of spike and bubble are expanding in horizontal direction due to the formation of vortex in the corresponding position. The velocity on Y axis direction is in a symmetric form with Y axis, show the respective motion trend of bubble and spike under gravity and buoyancy. Then the quantitative analysis of RTI development is carried out. Fig. 16 shows the amplitude variations of lighter and heavier fluids during the whole simulation of Case1 and Case2. The red triangles represent the bubble rising, the orange circles represent the spike intrusion, black line represents the theoretical solution of the bubble development process at the linear growth stage, blue dashed line is the theoretical solution of spike development process at the linear growth stage. In the linear growth stage, the amplitude of bubble and spike is symmetrical around the interface. According to the nonlinear threshold theory of the development of RTI, when the interface disturbance amplitude is 0.1 times longer than the initial disturbance wavelength (i.e. the nonlinear threshold shown in black and blue dotted lines), RTI development have been into nonlinear growth stage (Jacobs and Catton, 2006). Therefore the quantitative development process can be divided into linear and nonlinear growth stage. The linear growth stage is as shown in the area ① of the figure, as shown in the figure, in the area ①, the simulation results conforms well with the theoretical solution, which proves the correctness of the simulation in the RTI linear growth process. By contrast of the linear growth stage under two conditions, it can be found that although the RTI development process of Case2 is faster, but the time that Case2 enter into the nonlinear growth stage is later than that of Case1, this is because Case2 has greater initial disturbance wavelength, but according to the linear stage development rate calculated from Eq. (25), n1 = 6.5, n2 = 5.22, which means the linear growth rate of Case1 is bigger than that of Case2, so Case1 leave from the linear

5.2.4. RTI simulation under different A In this section, the development of RTI under different Atwood numbers is simulated. Sttings of calculation domain under different operating conditions are the same as those of Case2. The other parameters are shown in Table 1. As shown in Figs. 17 and 18, the RTI development under different A numbers are shown, including the amplitude and velocity variation of bubbles and spikes. Simulation of each case has captured the moment of 142

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

the heavy fluid arriving at the bottom of calculation area. As shown in figures, the development of RTI become quicker with the increase of A, which is because the larger density difference between two fluids is, the stronger the force drives RTI emergence. In addition, with the increase of A, the time of spike arriving at the bottom of the calculation area is accelerated, so the duration of RTI in the steady growth stage is shorter, but we can still observe the steady development process under different conditions. The comparison of the steady velocity constant and analytical solutions under different A is shown in Table 3. According to Eq. (33), when the initial disturbance is the same, the steady velocity of the bubble and spike should be related to the density difference. The bigger A number is, the bigger the ratio of C1,s and C1,b is, which means the spike's steady velocity is larger than that of the bubble. This is same with the situation reflected in Figs. 17 and 18. From the comparison of simulation values and analytical values of steady velocity constant under different A numbers, it can be concluded that the MMPS method can achieve convincing simulation results for most cases, except in the case of A = 0.5, spikes steady velocity constant has the error of 2.17%.

with 0.3m wavelength is faster, and develop into the nonlinear development stage earlier. Then there is bubble mergence emerging. As shown in Fig. 22, the bubble with the smaller wavelength decreases gradually, and eventually merged by larger wavelength bubbles. Fig. 23 shows a more complex RTI development process with seven mode disturbance. The wavelengths of the seven initial disturbances are 0.08m, 0.1m, 0.12m, 0.14m, 0.16m, 0.18m and 0.22m, respectively. It can be observed clearly from Fig. 23 that, larger bubbles' develop velocity are faster and gradually merge bubbles with smaller wavelengths. In the final snapshot, it is difficult to distinguish their initial wavelengths. Therefore, we can conclude that due to the competition of bubbles in the interface instability development process of multimode disturbance, the disturbance wave became longer and longer, interface formed is developing accelerated, and the initial disturbance wavelengths were gradually forgotten. Our numerical simulation results are consistent with those of other researchers adopting Finite Difference Methods and the Interface Tracking Technique (Youngs, 1984; Glimm and Li, 1988). 6. Conclusion

5.2.5. RTI simulation under different surface tension coefficients For the Case2 shown in Table 1, when the linear growth rate is zero, the critical surface tension coefficient σc=0.507N/m can be calculated by Eq. (30), and conditions of Case6∼9 in Table 1 are designed accordingly. Figs. 19 and 20 show the development of RTI under different surface tension coefficients. As shown in figures, the greater the surface tension, the slower the development of RTI. Especially when Φ is 1, according to the RTI linear development stage theory, the linear development rate is 0, which means the RTI process will not happen. The green dashed lines shown in Figs. 19 and 20 indicate that whether the lighter fluid or heavier fluid, there is almost no change for its amplitude in the whole calculation process, while its velocity have shocks in a certain range due to the assignment of initial velocity in the simulation. But in the RTI simulation article by Shadloo M S et al. with SPH method (Shadloo et al., 2013), the emergence of RTI process is not prevented until they set Φ as the value of around 1.15, the deviation to the theoretical value is about 15%, which is mainly because the initial displacement disturbance adopted in their article, and there are obvious numerical errors in CSF model used in the article. In this paper, the initial velocity disturbance and the CCSF model are used to avoid these two problems. Therefore, the simulation results which conform well to the theoretical results are obtained. In addition, by comparing of the velocity of bubble and spike, it can be found that although the different surface tension coefficients lead to the differences of RTI development, steady velocity obtained is the same under different conditions, which means the steady development velocity is not affected by the surface tension under different conditions.

In this paper, the Rayleigh-Taylor instability is simulated by the MMPS method. Firstly, the influence of initial interface disturbance on the development process of RTI is discussed. It is found that the initial displacement disturbance will cause greater error in larger wavelength, therefore, this paper puts forward the RTI simulation by the initial velocity disturbance. In terms of boundary conditions, the achievement method of periodic boundary conditions based on the particle method is developed, taken as one of the boundary conditions of RTI simulation in this paper. Then the RTI process with different initial disturbance wavelengths, different density ratios and different surface tension coefficients were simulated, the linear growth stage and steady stage of the RTI process were analyzed, and the corresponding results were compared with the RTI theoretical solutions, which proved the correctness of MMPS method in RTI simulation. Finally, for the simulation of RTI with multimode disturbance, multi bubbles' competition and mergence process are observed in simulation results. In the future, we will continue our research about the RTI processes, especially the experimental research, to further study the theory of the development of RTI processes. Acknowledgements The present study is supported by the National Science Foundation of China (No. 11505134), the National Natural Science Foundation of China outstanding youth fund (No. 11622541) and by the China Postdoctoral Science Foundation (no. 2018M633522). I thank Dr. Guangtao Duan for his excellent work in the multiphase MPS method and a preprint of Ref. (Duan et al., 2017)

5.3. Numerical simulation of RTI with multimode disturbance References In the last part of this paper, RTI cases with initial multimode disturbance are simulated. Through the discussion about the initial condition in Section 3, when the initial wavelength meets the criteria shown in Eq. (21), initial displacement disturbance and initial velocity disturbance will not have too much error. So in the multimode disturbance simulation of this section, we adopt the arrangement of the initial displacement disturbance. The amplitude of each initial disturbance is the spatial scale l0 of the particle. Take the physical properties of different fluids as ρ1=1kg/m3, ρ2=3kg/m3, μ1=μ2=0.00334 Pa s, g=0.04m/s2。 The dual mode disturbance RTI process with initial disturbance wavelengths of 0.2m and 0.3m respectively is shown in Fig. 21. As shown in the figure, two disturbances develop in the form of Cosine at the initial linear growth stage. The disturbance with larger wavelength has larger linear growth rate, therefore the disturbance development

Abarzhi, S.I., 1998. Stable steady flows in Rayleigh-Taylor instability. Phys. Rev. Lett. 81 (2), 337–340. Cummins, S.J., Rudman, M., 1999. An SPH projection method. J. Comput. Phys. 152 (2), 584–607. Dalziel, S.B., 2000. Self-similarity and internal structure of turbulence induced by, Rayleigh–Taylor instability. J. Fluid Mech. 399 (399), 1–48. Danielson, R.E., 1962. Hydrodynamic and hydromagnetic stability by S. Chandrasekhar. Am. Sci. 4 428A,430A. Dimonte, G., 2000. Spanwise homogeneous buoyancy-drag model for Rayleigh–Taylor mixing and experimental evaluation. Phys. Plasmas 7 (6), 2255–2269. Dimonte, G., Ramaprabhu, P., Andrews, M., 2007. Rayleigh-Taylor instability with complex acceleration history. Phys. Rev. E - Stat. Nonlinear Soft Matter Phys. 76 (4 Pt 2) 046313. Dinh, T.N., Bui, V.A., Nourgaliev, R.R., et al., 1999. Experimental and analytical studies of melt jet-coolant interactions: a synthesis. Nucl. Eng. Des. 189 (1–3), 299–327. Duan, G., Chen, B., 2013. Stability and accuracy analysis for viscous flow simulation by the moving particle semi-implicit method. Fluid Dynam. Res. 45 (3), 282–294. Duan, G., Koshizuka, S., Chen, B., 2015. A contoured continuum surface force model for

143

Progress in Nuclear Energy 109 (2018) 130–144

K. Guo et al.

incompressible, multiphase flows. J. Comput. Phys. 202 (1), 65–93. Meyer-ter-Vehn, Jürgen, 2009. The Physics of Inertial Fusion. Oxford University Press. Mikaelian, K.O., 1993. Effect of viscosity on Rayleigh-Taylor and Richtmyer-Meshkov instabilities. Phys. Rev. E: Stat. Phys., Plasmas, Fluids, Relat. Interdiscip. Top. 47 (1), 375. Price, D.J., 2007. SPLASH: An interactive visualisation tool for smoothed particle hydrodynamics simulations. Publ. Astron. Soc. Aust. 24 (3), 159–173. Ratafia, M., 1973. Experimental investigation of Rayleigh-Taylor instability. 16 (8), 1207–1210. Read, K.I., 1992. Experimental investigation of turbulent mixing by Rayleigh-Taylor instability. Phys. Nonlinear Phenom. 12 (1), 45–58. Reckinger, S.J., Livescu, D., Vasilyev, O.V., 2016. Comprehensive numerical methodology for direct numerical simulations of compressible Rayleigh–Taylor instability. J. Comput. Phys. 313, 181–208. Saito, S., Abe, Y., Kaneko, A., et al., 2015. Experimental study on jet instability and breakup behavior in liquid-liquid system. In: International Conference on Nuclear Engineering. Shadloo, M.S., Zainali, A., Yildiz, M., 2013. Simulation of single mode Rayleigh—Taylor instability by SPH method. Comput. Mech. 51 (5), 699–715. Shakibaeinia, A., Jin, Y.C., 2012. MPS mesh-free particle method for multiphase flows. Comput. Meth. Appl. Mech. Eng. 229–232 (2), 13–26. Shao, S., Lo, E.Y.M., 2003. Incompressible SPH method for simulating Newtonian and non-Newtonian flows with a free surface. Adv. Water Resour. 26 (7), 787–800. Tamai, T., Koshizuka, S., 2014. Least squares moving particle semi-implicit method. Comput. Part. Mech. 1 (3), 277–305. Tartakovsky, A.M., Meakin, P., 2005. A smoothed particle hydrodynamics model for miscible flow in three-dimensional fractures and the two-dimensional Rayleigh–Taylor instability. J. Comput. Phys. 207 (2), 610–624. Taylor, G., 1950. The instability of liquid surfaces when accelerated in a direction perpendicular to their planes. I. Proc. Roy. Soc. Lond. 201 (1065), 192–196. Youngs, D.L., 1984. Numerical simulation of turbulent mixing by Rayleigh-Taylor instability. Phys. Nonlinear Phenom. 12 (1), 32–44. Zainali, A., Tofighi, N., Shadloo, M.S., et al., 2013. Numerical investigation of Newtonian and non-Newtonian multiphase flows using ISPH method. Comput. Meth. Appl. Mech. Eng. 254 (2), 99–113. Zhang, S., Morita, K., Shirakawa, N.[E.A., 2004. Simulation of the Rayleigh-Taylor instability with the MPS method. Mem. Facul. Eng. Kyushu Univ. 64 (4), 215–228.

particle methods. J. Comput. Phys. 298 (C), 280–304. Duan, G., Chen, B., Koshizuka, S., et al., 2017. Stable multiphase moving particle semiimplicit method for incompressible interfacial flow. Comput. Meth. Appl. Mech. Eng. 318, 636–666. Glimm, J., Li, X.L., 1988. Validation of the Sharp–Wheeler bubble merger model from experimental and computational data. 31 (8), 2077–2085. Grenier, N., Antuono, M., Colagrossi, A., et al., 2009. An Hamiltonian interface SPH formulation for multi-fluid and free surface flows. J. Comput. Phys. 228 (22), 8380–8393. He, X., Chen, S., Zhang, R., 1999. A lattice Boltzmann scheme for incompressible multiphase flow and its application in simulation of Rayleigh-Taylor instability. J. Comput. Phys. 152 (2), 642–663. Heo, H., Bang, I.C., Park, S.D., et al., 2015. Visualization study on molten metal jet breakup phenomena in horizontal and vertical injections for SFR with metallic fuel. In: International Conference on Nuclear Engineering. Hu, X.Y.,N.A., 2006. A multi-phase SPH method for macroscopic and mesoscopic flows. J. Comput. Phys. 213 (2), 844–861. Hu, X.Y., Adams, N.A., 2007. An incompressible multi-phase SPH method. J. Comput. Phys. 227 (1), 264–278. Jacobs, J.W., Catton, I., 2006. Three-dimensional Rayleigh-Taylor instability Part 1. Weakly nonlinear theory. J. Fluid Mech. 187 (187), 329–352. Jeong, S.M., Nam, J.W., Hwang, S.C., et al., 2013. Numerical prediction of oil amount leaked from a damaged tank using two-dimensional moving particle simulation method. Ocean Eng. 69 (5), 70–78. Khayyer, A., Gotoh, H., 2011. Enhancement of stability and accuracy of the moving particle semi-implicit method. J. Comput. Phys. 230 (230), 3093–3118. Khayyer, A., Gotoh, H., 2010. A higher order Laplacian model for enhancement and stabilization of pressure calculation by the MPS method. Appl. Ocean Res. 32 (1), 124–131. Kim, B., Corradini, M.L., 1988. Modeling of small-scale single droplet fuel/coolant interactions. Nucl. Sci. Eng.; (USA) 98:1 (1), 16–28. Koshizuka, S., Oka, Y., 1996. Moving-particle semi-implicit method for fragmentation of incompressible fluid. Nucl. Sci. Eng. 123 (3), 421–434. Kull, H.J., 1991. Theory of the Rayleigh-Taylor instability. Phys. Rep. 206 (5), 197–325. Li, X.L., 1993. Study of three‐dimensional Rayleigh–Taylor instability in compressible fluids through level set method and parallel computation. Phys. Fluid. Fluid Dynam. 5 (8), 1904–1913. Liu, J., Koshizuka, S., Oka, Y., 2005. A hybrid particle-mesh method for viscous,

144