CSF model performance for Rayleigh-Taylor instability

CFD MODEL PERFORMANCE FOR RAYLEIGH–TAYLOR

INSTABILITY

Sergey N. Yakovenko* and K.C. Chang† * Laboratory of Modeling of Turbulent Flows

Institute of Theoretical and Applied Mechanics (ITAM) SB RAS Institutskaya Street, 4/1, Novosibirsk 630090, Russia

e-mail: yakovenk@itam.nsc.ru

Web page: http://www.itam.nsc.ru/~yakovenk/index.htm

†

Department of Aeronautics and Astronautics National Cheng Kung University (NCKU)

Tainan 70101, Taiwan, R.O.C. e-mail: kcchang@mail.ncku.edu.tw

Key words: Rayleigh-Taylor instability, Interface resolution, Continuum surface force model, Kernel function

Abstract. The surface tension effect is introduced according to the continuum surface force model. The mollified volume fraction (color) function varies smoothly across the interface due to convolving the original color function with a smooth kernel. The eight-order polynomial kernel is formulated for plane 2D flows constricted by solid or symmetry boundaries and tested for a Rayleigh-Taylor instability problem.

1 INTRODUCTION

In the present study, the surface tension is parameterized by volume forces due to the continuum surface force (CSF) model1 using the surface tension coefficient σ and the mollified volume fraction (color) function. Smooth variation of this function across the interface is achieved by convolving the volume fraction with an interpolation function K. The properties (compactness, monotonic decreasing, symmetry, differentiation, normality, limit behavior) of K and some examples of its choosing are presented in the paper1.

boundaries and tested for a Rayleigh-Taylor instability (RTI). The RTI flow is a benchmarked problem in free surface hydrodynamics chosen in other theoretical3 and numerical1,4-10 studies too.

2 GOVERNING EQUATIONS

The governing equations for a two-dimensional (2D) incompressible laminar two-fluid flow taking into account surface tension effects are

0 U V x y ∂ ₊∂ ₌ ∂ ∂ (1) 2 1 1 1 S U U U UV p U U F t x y ρ x ρ x η x ρ y η y ⎡ ⎤ ⎡ ⎤ ∂ ₊∂ ₊∂ _{= −} ∂ ₊ ∂ ∂ ₊ ∂ ∂ ₊ ⎢ ⎥ ⎢ ⎥ ∂ ∂ ∂ ∂ ∂ _⎣ ∂ _⎦ ∂ _⎣ ∂ _⎦ (2) 2 1 1 1 S V V UV V p V V g F t x y ρ y ρ x η x ρ y η y ⎡ ⎤ ⎡ ⎤ ∂ ∂ ∂ ∂ ∂ ∂ ∂ ∂ + + = − + _{⎢ ⎥}+ _{⎢ ⎥}− + ∂ ∂ ∂ ∂ ∂ _⎣ ∂ _⎦ ∂ _⎣ ∂ _⎦ (3)

(

1f 2 1 f

)

ρ ρ= +ρ − , η η= ₁f +η₂

(

1− f

)

(4)

( )

( )

0 Uf Vf f t x y ∂ ∂ ∂ _{+ + =} ∂ ∂ ∂ (5)

where U and V are the horizontal and vertical components of the velocity vector, respectively, x and y are the horizontal and vertical coordinates, respectively, g is gravity acceleration, ρ is density, p is pressure, t is time. The values ρ1 and η1 correspond to heavier fluid (liquid) at

, whereas the values

1

f = ρ2 and η2 correspond to lighter fluid (gas) at f =0.

In velocity equations the surface tension effects are introduced as volume forces according to the continuum surface force (CSF) model1 where reformulated volumetric force is expressed (for general 3D case) as

[ ]

,

[ ]

2 1 1 S i i f F f f f f x σκ ρ ∂ = = − = ∂  . Here σ is coefficient of surface tension between fluids 1 and 2, i

i n x κ = −∂ ∂  is surface curvature, i i f n x ∂ = ∂   _{is the unit}

normal to the surface. The mollified volume fraction (color) function can be taken as the first approximation to be equal to the volume fraction function defined numerically due to (possible) smoothing of numerical schemes. However, accurate representation of the mollified color function assumes that it is varies smoothly over a thickness h across the interface by convolving the characteristic function f with an interpolation function (smooth kernel) K.

( )

( )

( ) (

)

d d d1 2 3 K i i i i i f x f x f x x x x x Ω ′ ′ ′ ′ x = ∗K =

∫

⋅K −  _{′ (6)}

A new eight-order polynomial kernel denoted K8 was introduced

2 to be equal to

( )

( )

4 2 8 1 if , 0 if K C r r r r ε ε ε ε ⎧ _{⎡ − ⎤ <} ⎪ _{⎣ ⎦} = ⎨ ⎪ ≥ ⎩ K (7)

with normalization condition ₈

( )

, d 1

K

r ε r

Ω =

∫

K defining the constant A. The parameter ε is

the radius of the region ΩK over which K is nonzero, and

(

) (

2

) (

2

)

1 1 2 2 3 3

r= x′−x + x′−x + x′−x 2

. The simple form of results in a trivial and computationally efficient implementation, which is important since numerical approximation in (6) can require many kernel evaluations.

8

K

For plane 2D flows the convolution (6) can be written as

( )

,

(

,

) (

,

)

d K d f x y f x y x x y y x y Ω ′ ′ ′ ′ ′ =

_∫∫

K − −  _{′ (8)} where

(

)

( )

( )

4 2 8 1 if , , 0 if K C r r x x y y r r ε ε ε ε ⎧ _{⎡ − ⎤ <} ⎪ _{⎣ ⎦} ′ ′ − − = _{= ⎨} ⎪ ≥ ⎩ K K (9)

and r=

(

x′−x

) (

2+ y′−y

)

2 . Normalization condition gives

( )

( )

₂ 4 2 8 0 1 , d d 1 2 d 5 K K K r x y C r r r C ε _π ε ε π Ω ⎡ ⎤ =

∫

K =

∫

_⎣ − _⎦ = ε , and CK =_πε5₂ .

To define smoothed color function f x y

( )

, correctly in the whole (rectangular) computational domain located at x0≤ ≤x xN and y0 ≤ ≤y yN one needs to know values of

(

,

)

f x y′ ′ in the extended region x0− ≤ ≤ε x′ xN +ε and y0− ≤ε y′≤yN +ε . Thus, we should

introduce virtual values of f at x₀− ≤ <ε x′ x₀, x_N < ≤x′ x_N +ε and y₀− ≤ε y′< y₀,

N N

(

,

)

(

2 0 ,

)

f x y′ ′ = f x −x y′ ′ for x0− ≤ <ε x′ x0. On the other hand, for the solid wall boundary located at x x= 0 we will have f x y

(

′ ′,

)

= f x y

(

0, ′

)

for x0− ≤ <ε x′ x0, assuming zero derivative

(

)

0 0

x x

f x ₌

∂ ∂ = normal to the wall (if information on dynamic contact angles is not available).

3 NUMERICAL METHODS AND SCHEMES

The Navier – Stokes solver is described in the paper11 where it was applied for a broken dam problem flow. The time step can be chosen to be sufficiently small and follows the scheme

(

)

* 0 1 0 0 t t ⎡C t t t t ⎤ t ∆ =_⎣ − + _⎦∆ , (C tt0, )0 =(0.1, 0.5) at t t≤ , 0 * t t ∆ = ∆ at t t≥ 0

(

)

2 2 1 * 1 3 1 2 2 2 1/ 2, , 1/ 2 0.01 min , , , , 4 2 j j i j i i r i j i j i j y y x y x x t U U V x y ρ ρ ρ δ _δ πσ η + + + + ⎧ _{− − ∆ ∆ +} ⎫ ⎪ ⎪ ∆ = ⋅ _{⎨ ⎬} ∆ + ∆ ⎪ ⎪ ⎩ ⎭ , δ =min

{

x_i+1−x y_i, _j+1−y_j

}

. Note, the last term in the expression is due to stability condition which results in the time step constraint imposed by an explicit treatment of surface tension.

*

t

∆

An hierarchy of five numerical schemes for advection terms applied for volume fraction and velocity equations was presented in the paper11. The first-order upwind scheme and third-order QUICK scheme were shown to be unsuitable for free-surface flows because of generation of the thick unrealistic interface due to numerical diffusion of the upwind scheme, and significant distortions due to spurious oscillation features of the QUICK scheme. So, three remaining workable schemes containing the MUSCL scheme with QUICK interpolants and with compressive minmod TVD limiter for U and V equations combined with the same scheme for f (denoted as c4mm) or with slope modification scheme for f (denoted as SMF) or with donor-acceptor upwind-downwind scheme (denoted as VOF) are used in the present work together with appropriate approximations for diffusion terms and for CSF terms.

For the staggered grid used, the latter should be defined at the centers of computational cell surfaces in the equations for the components Ui±1/2,j and Vi,j±1/2 of the mean velocity vector, via the

scalar quantity f determined at the cell center with the coordinate _{i j,}

(

x y : _i, _j

)

( )

( )

,

( )

( )

,

( )

( )

1/ 2, 1/ 2

( )

( )

1/ 2, 1/ 2 , , 1/ 2, 1/ 2 1/ 2, 1/ 2 , , 1/ 2, 1/ 2 1/ 2, 1/ 2 , y , , y x i j i j x i j i j x y x y i j i j i j i j i j i j i j i j n n n n n n n n n n n n + + + + + + + + + + + + = = = =  _G  _G  _G  _G

( )

2

( )

2

( )

2

( )

2 , 1/ 2, 1/ 2 , x i j y , , 1/ 2, 1/ 2 x i j y 1/ 2, 1/ 2 i j i j i j i j n n n n _{+ +} n _{+ +} n + + = + = + G G

( )

_, 1

( )

_{1/ 2,}

( )

_{1/ 2,}

( )

_, 1

( )

_{, 1/ 2}

( )

_{, 1/ 2} , 2 2 x _{i j} x _{i j} x _{i j} y _{i j} y _{i j} y _{i j} n n ₊ n ₋ n n n + − ⎡ ⎤ ⎡ ⎤ = _⎣ + _⎦ = _⎢ + _⎥ ⎣ ⎦

( )

_{1/ 2, 1/ 2} 1

( )

_{1/ 2,}

( )

_{1/ 2, 1}

( )

_{1/ 2, 1/ 2} 1

( )

_{, 1/ 2}

( )

_{1, 1/ 2} , 2 2 x _{i j} x _{i j} x _{i j} y _{i j} y _{i j} y _{i j} n _{+ +} n ₊ n _{+ +} n n n + + + + + ⎡ ⎤ ⎡ ⎤ = _⎣ + _⎦ = _⎢ + _⎥ ⎣ ⎦

( )

1, ,

( )

, 1 , 1/ 2, , 1/ 2 1 1 , i j i j i j i j x _{i j} y _{i j} i i j j f f f n n f x x y + + + + + + − − = = − −     y

4 APPLICATION TO RAYLEIGH-TAYLOR INSTABILITY PROBLEM

The above formulated numerical methods are tested during computations of a Rayleigh-Taylor instability. The initial state of the sudden collapse of rectangular water column is schematically shown in Fig. 1. The main calculation conditions follow the test case4 of the RTI problem: L=0.02, H =1.5L, g=1,

1 2

ρ = , ρ2=1. The coefficients of viscosity and surface tension are varied (see in Table 1). All variables are nondimensionalized using the reference quantities ρr =ρ1, r 1 1

↑ y ↓

g

G

=

(0,

−

g

)

heavier fluid (

ρ ν

₁

,

₁

)

↑ +H ↓ t

=

0

2

,

2

↓ -H

lighter fluid (

ρ ν

)

→ x

0

x

=

x L

=

Figure 1: Scheme of RTI problem (initial state)

η =η =ρ ν,

r

L =L, U_r = gL, t_r = L g producing the following reference Reynolds number

3 Re_r r r r r gL U L ρ η ν = = .

Initial conditions are taken to be the same as in the papers4,5,8-10. A single wavelength perturbation of amplitude is introduced at the fluid interface using the vertical velocity distribution for the whole computational domain in the form

(

, ,0

)

cos x exp y V x y L L α ⎛π ⎞ ⎛ π ⎞ = ⋅ _{⎜ ⎟ ⎜}− _⎟ ⎝ ⎠ _{⎝ ⎠}, 2 LA π δ α =

and the following distributions for the heavier fluid (at y>0):

(

, ,0

)

sin x exp y U x y L L α ⎛π ⎞ ⎛ π ⎞ = ⋅ _{⎜ ⎟ ⎜}− _⎟ ⎝ ⎠ _{⎝ ⎠}, f x y

(

, , 0

)

=1, p x y

(

, , 0

)

= −ρ1gy and for the lighter fluid (at y<0):

(

, ,0

)

sin x exp y U x y L L α ⎛π ⎞ ⎛ π ⎞ = − ⋅ _{⎜ ⎟ ⎜}− _⎟ ⎝ ⎠ _{⎝ ⎠}, f x y

(

, , 0

)

=0, p x y

(

, ,0

)

= −ρ2gy

except for the thin transition region between two fluids (at y=0) where . The initial velocity field is the conservative form

(

, , 0

)

0.5

f x y =

4,5

of the (non-conservative) sinusoidal vertical velocity perturbation of amplitude A and wavelength λ=2L supplied at the interface by specifying V x y

(

, , 0

)

cos x L α ⎛π = ⋅ _⎜ ⎝ ⎠ ⎞

⎟. Note, the modified Reynolds number

3 Re_m r r U_λ g ρ λ λ η ν

= = can be built using the wavelength in reference velocity and length scales. The initial pressure corresponds to a hydrostatic distribution, and negative pressure values should not confuse anyone due to presence of pressure derivatives only in the pressure-velocity computations. The initial volume fraction distribution is assigned to get the density field shown in Fig. 1 from the equality (4).

Left and right vertical boundaries are taken to be symmetry planes, whereas top and bottom horizontal boundaries are considered as solid walls. The pressure boundary conditions are derived8-11 from the normal velocity equation, and the virtual values in the fictitious cells surrounding the mesh are used11. The volume fraction on the walls is extrapolated from the flow interior8-11, whereas the normal velocity is set to zero.

5 RESULTS OF NUMERICAL EXPERIMENTS

Following the previous studies4,5,8-10_{, the computational domain is extended from} 0

x= up to x L= (Fig.1), whereas lower and upper boundaries are set at y= −1.5L and y=1.5L, respectively. The uniform grid 41x121 with δ = H/40 is used for the base computations8-10_. Such a refinement was shown in grid independent studies11

to be sufficient. The main stages of numerical studies performed are listed in the Table 1 as different Cases and discussed below.

coefficient A varied from 1.00 up to 0.02. The simulations show (Fig.2) that the average y_f

of amplitudes of spike and bubble (seen in the right and left sides of the domain in Fig.3) has the clear exponential growth at A<0.1 only. It is realized mainly at − <3 ln

(

y L_f

)

< −2

t t

. This stage corresponds to linear instability3-5 _{which can be defined as or .}

Therefore, one can determine the growth rate

( )

~ exp f y n lny_f ~n (ln _f) d y n dt

= calculated using the difference of

lny_f at the range − ≤3 ln

(

y L_f

)

≤ −2. This rate increases with Reynolds number (Fig.4). However, the non-dimensional growth rate _n*=_n

(

ν _g2

)

1/ 3

demonstrates the non-monotonic behaviour in full correspondence with theoretical prediction3 _{and computations of other} authors4,8-10_{. It can also be seen (Fig.3) that viscosity damps the Rayleigh-Taylor instability} development. The effect of different advection schemes is (Fig.6) the same that for dam-break problem11_.

The effect of the surface tension consists in delay of the RTI development (Fig.7) and growth rate (Fig.8,9). This is also in full correspondence with theoretical curve [4] and computations of other authors5_{. Note, the using of the kernel function with appropriate} thickness gives the dramatic improving effect in the computations (Fig.9).

6 CONCLUSIONS

- In the present paper the CSF model with an eight-order polynomial kernel is applied to the Rayleigh-Taylor instability problem.

- The damping effects of both viscosity and surface tension are described in full correspondence with theoretical predictions.

- The used approach is to be applied to other laminar and turbulent flows with noticeable surface tension effects to clarify limitations and capabilities of the CSF model.

7 ACKNOLEDGEMENTS

The authors gratefully acknowledge the grant support from the National Science Council of the Republic of China for this work under the contract NSC-92-2212-E006-102.

REFERENCES

[1] J.U.Brackbill, D.B.Kothe and C.Zemach, “A continuum method for modeling surface tension”, J. Comp. Phys., 100, 335-354 (1992).

[2] M.W.Williams, D.B.Kothe and E.G.Puckett, “Accuracy and convergence of continuum surface-tension method”, Fluid Dynamics at Interface, Eds. Wei Shyy, Ranga Narayanan, Cambridge Univ. Press, Cambridge, UK, 294-305 (1999).

[3] S.Chandrasekar, Hydrodynamic and hydromagnetic stability, Clarendon Press, Oxford, § 94, 441-453 (1961).

10, 297-307 (1967).

[5] B.J.Daly, “Numerical study of the effect of surface tension on interface instability,” Phys. Fluids, 12, 1340-1354 (1969).

[6] C.W.Hirt and B.D.Nichols, “Volume of fluid (VOF) method for the dynamics of free boundaries”, J. Comp. Phys., 39, 201-225 (1981).

[7] M.Rudman, “Volume-tracking method for interfacial flow calculations”, Int. J. for Numerical Methods in Fluids, 24, 671-691 (1997).

[8] F.J.Kelecy and R.H.Pletcher, “The development of a free surface capturing approach for multidimensional free surface flows in closed containers,” J. Comp. Phys., 138, 939-980 (1997).

[9] C.H.Chang, “The development of a multi-fluid surface-capturing method for the computation of incompressible free surface flows,” Ph.D.Thesis, NCKU, Tainan (1998).

[10] D.Pan and C.H.Chang, “A free surface capturing method for incompressible multi-fluid flows,” Proc. of the 3rd ASME/JSME Joint Fluids Engng Conf., FEDSM99-7105, 1-6 (1999).

Cases A Re_m σ Kernel Advection scheme 1 1.00 800 0 no K c4mm 2 0.50 800 0 no K c4mm 3 0.20 800 0 no K c4mm 4 0.10 800 0 no K c4mm 5 0.05 800 0 no K c4mm 6 0.02 800 0 no K c4mm 7 0.05 176 0 no K c4mm 8 0.05 72 0 no K c4mm 9 0.05 39 0 no K c4mm 10 0.05 20 0 no K c4mm 11,12,13 0.05 400 0 no K c4mm, SMF, VOF 14,15,16 0.05 400 0.00001 no K c4mm, SMF, VOF 17,18,19 0.05 400 0.00002 no K c4mm, SMF, VOF 20,21,22 0.05 400 0.00003 no K c4mm, SMF, VOF 23,24,25 0.05 400 0.00004 no K c4mm, SMF, VOF 26,27,28 0.05 400 0.00005 no K c4mm, SMF, VOF 29,30,31 0.05 400 0.00006 no K c4mm, SMF, VOF 32,33,34 0.05 400 0.00007 no K c4mm, SMF, VOF 35,36,37,38 0.05 400 0.000010 K8,ε =0.1, 0.2, 0.3, 0.4 VOF 39,40,41,42 0.05 400 0.000020 K8,ε =0.1, 0.2, 0.3, 0.4 VOF 43,44,45,46 0.05 400 0.000030 K8,ε =0.1, 0.2, 0.3, 0.4 VOF 47,48,49,50 0.05 400 0.000035 K8,ε =0.1, 0.2, 0.3, 0.4 VOF 51,52,53,54 0.05 400 0.000040 K8,ε =0.1, 0.2, 0.3, 0.4 VOF 55,56,57,58 0.05 400 0.0000410, 415, 418, 420 K8,ε =0.1 VOF 59,60 0.05 400 0.0000430, 0.0000445 K8,ε=0.2 VOF 61,62,63 0.05 400 0.0000430, 450, 460 K8,ε=0.3 VOF 64,65 0.05 400 0.0000450, 0.0000500 K8,ε=0.4 VOF

(a) t/t_r y f /L 0 2 4 6 8 0 1 c4mm, A=1.00 c4mm, A=0.50 c4mm, A=0.20 c4mm, A=0.10 c4mm, A=0.05 c4mm, A=0.02 (b) t/t_r ln (y f /L ) 0 2 4 6 8 10 -5 -4 -3 -2 -1 0 c4mm, A=1.00 c4mm, A=0.50 c4mm, A=0.20 c4mm, A=0.10 c4mm, A=0.05 c4mm, A=0.02

x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99

(a) t/t_r y f /L 0 2 4 6 8 10 0 0.5 1 1.5 c4mm, Rem= 800 c4mm, Rem= 400 c4mm, Rem= 176 c4mm, Re_m= 72 c4mm, Re_m= 39 c4mm, Rem= 20 (b) t/t_r ln (y f /L ) 0 2 4 6 8 10 -5 -4 -3 -2 -1 0 c4mm, Rem= 800 c4mm, Re_m= 400 c4mm, Rem= 176 c4mm, Rem= 72 c4mm, Rem= 39 c4mm, Rem= 20

Re

_m

n*

0 250 500 750 0 0.05 0.1 0.15 0.2 Theory c4mm

Figure 5: Growth rate of Rayleigh-Taylor instability: theoretical curve3 (solid line) and calculated values plotted as points (symbols) for different Cases versus the Reynolds number

(where n d(lnyf)

( )

dt n

x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/L 0 1 -1 0 1 2 0.01 0.50 0.99

x/W y/H 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/ L 0 1 -1 0 1 2 0.01 0.50 0.99 x/W y/H 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/ L 0 1 -1 0 1 2 0.01 0.50 0.99 x/L y/ L 0 1 -1 0 1 2 0.01 0.50 0.99

(a) t/t_r y f /L 0 5 10 0 0.5 1 1.5 T = 0.0000200 (+K8) T = 0.0000350 (+K8) T = 0.0000400 (+K8) T = 0.0000418 (+K₈) T = 0.0000420 (+K₈) T = 0 T = 0.00004 T = 0.00005 T = 0.00006 T = 0.00007 (b) t/t_r ln (y f /L ) 0 5 10 -5 -4 -3 -2 -1 0

Φ

n/

n

0

0 0.5 1 1.5 0 0.4 0.8 1.2 theory VOF without K c4mm without K SMF without K

VOF with K₈, eps=0.1 VOF with K₈, eps=0.2 VOF with K₈, eps=0.3 VOF with K₈, eps=0.4

Figure 9: Growth rate of Rayleigh-Taylor instability: theoretical curve3 (solid line) and calculated values plotted as points (symbols) for different Cases versus the stability parameter

(

)

(

)

2 2 1 L g σ π ρ ρ Φ = − (where

( )

(ln f) d y n dt