A development of a sharp interface AUSMD scheme for an incompressible preconditiong multi-fluid model

Computers and Fluids 192 (2019) 104269

Contents lists available at ScienceDirect

Computers and Fluids journal homepage: www.elsevier.com/locate/compfluid

Benchmark solutions

A development of a sharp interface AUSMD scheme for an incompressible preconditiong multi-fluid model Yang-Yao Niu∗, Chia-Hung Weng Department of Aerospace Engineering, Tamkang University, New Taipei City, Taiwan, R.O.C.

a r t i c l e

i n f o

Article history: Received 29 August 2018 Revised 4 August 2019 Accepted 19 August 2019 Available online 20 August 2019 Keywords: Preconditiong AUSM Incompressible Flow Interface Capturing MUSCL+THINC

a b s t r a c t In this study, the paper first proposes a sharp-interface cell interface variable reconstruction to solve the preconditioning incompressible multi-fluid flow model using the AUSMD scheme. Here, the preconditioning matrix modified from the Chorin’s artificial compressibility makes the incompressible multi-fluid equations to be hyperbolic. Therefore, a simple AUSMD type flux-splitting with the MUSCL+THINC (Tangent of Hyperbola for Interface Capturing) type cell interface variable reconstruction is used to maintain the preservation of sharp interface evolutions of the multi-fluid simulation. A high-resolution interpolation MUSCL+THINC, the MUSCL method conjunction with the THINC decided under the boundary variation diminishing (BVD) concept, is used to determine the cell interfaces to evaluate the inviscid flux. The interface revolution of the Rayleigh Tayor instability, the dam breaking and the rising bubbles are shown and investigated. © 2019 Elsevier Ltd. All rights reserved.

1. Introduction It is well known that the sharp capturing fluid interface schemes are highly required for simulating immiscible two-phase flows is One of the popular sharp interface methods is called as the Eulerian volume-of-fluid (VOF) methods [7–17] normally use a marker function like the volume fraction to indicate the position of the interface which is expected to resolve interfaces in few mesh cells. Similarly, another widely used front tracking method [7,28] is proposed to simulate unsteady multi-fluid flows in which a sharp interface or a front separates incompressible fluids of different density and viscosity is described. Likewise, level-set methods [18–21] utilize the level set (zero value) of the signed distance function as the marker function. As the behavior of this function is smooth at the interface, such methods generally allow sharper resolution of the curvature of the interface and a more accurate treatment surface-tension effects. However, the level-set does not rigorously satisfy mass conservation, which always induces unnecessary numerical errors. Therefore, a hybrid approach seen in the works of Olsson and Kreiss [20], Sussman and Puckett [21] has been used to combine the VOF model with the level set method to reduce mass-conservation errors and to achieve sharper interfaces. Similarly, an Eulerian type Ghost-fluid method [22] has been shown to provide extremely sharp interface capturing combining ∗

Corresponding author. E-mail address: [email protected] (Y.-Y. Niu).

https://doi.org/10.1016/j.compfluid.2019.104269 0045-7930/© 2019 Elsevier Ltd. All rights reserved.

with the level-set techniques. Recently, the algebraic reconstruction of the VOF function [13–15] using smooth basis functions has been proposed to allow the mesh-scale resolution of the interface for immiscible two-phase flows. However, the reconstruction of the VOF function through artificial algebraic assumptions may not be suitable for compressible two-phase flows due to the possible violation of mass conservation around the interface. Previously [30], we propose the third-order type extrapolations associated with highly compressive limiters using an ENO-type strategy similar to Pilliot and Puckett [11] coupling with Monotonic Upwind Scheme for Conservation Laws (MUSCL) type volume-fraction reconstructions based on extrapolated fluid properties to sharpen the resolution of the volume fraction near the interface. One alternative to deal with multi-scale phenomena in the two-phase flow problems relies on the use of time-derivative preconditioning strategies. Preconditioning schemes originating from Chorin’s artificial compressibility method [23] have been extended to simulate cavitated flow problems using AUSM+ and The AUSMD [24,25] and LDFSS-2001 [26] type all-speed, flux-splitting schemes. Here, the current study considers the extension of the flux splitting scheme to solve an incompressible immiscible two-phase flow model using a sharp-interface Tangent of Hyperbola for Interface Capturing technique (THINC) [15] method for the discontinuity resolution. The currents method develops the AUSMD scheme combining with the MUSCL+THINC formulation under the idea of boundary vale diminishing (BVD) in [33] is used to maintain the preservation of sharp interface evolutions in the immiscible two-phase

2

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

flow simulations. The first test case is a classical Rayleigh–Taylor problem proposed by Daly [33] which considers two fluids of unequal densities subject to a body force field directed normally to their common, planar interface. When the direction of acceleration is from the light to the heavy fluid, the interface is unstable to small perturbations, waves grow exponentially in time, and eventually the interface breaks up into a two-phase mixing zone that grows non-linear algebraically in time. Resolution of phaseinterface capturing will be demonstrated. Secondly, the dam break flow with the consequent wall impact is widely used to benchmark various numerical techniques to simulate the interfacial flows problems. The experimental data in the study of Martin and Moyce [35] is used to validate the propagation of the free surface due to the collapse of the dam. The third test case is the gas bubble rising in liquid which is one of the typical multi-component flows with sharp interfaces widely seen in many engineering problems. During the past decades, the single bubble rising in a viscous liquid attracts many studies. The early theoretical analysis are documented in [1,2]. Experiments and correlations were reviewed by Clift et al. [3]. Later, a more comprehensive study was introduced by Bhaga and Weber [4] who presented experimental data and correlations for the bubble shape and rise velocity. However, due to the insufficient techniques of measurements, the understanding of bubble rising and deformation is strongly restricted. In theoretical works [5–6], small bubble deformations under low Bond number for a wide range of Reynolds numbers were analyzed. In the work of Davies and Taylor [1], the estimated rising velocity of a sphericalcap bubble was derived. Numerical results of the first and second benchmark tests show that the MUSCL + THINC-type preconditioning AUSMD solver is superior to the original MUSCL-type preconditiong AUSMD in resolving interfaces. In the third cases, the revolution of interfaces at different Reynolds number, Bond number and Archimedes number is well simulated. Also, the grid independence study is performed. In our numerical simulations, comparing to experimental investigations has demonstrated the MUSCL+THINC type AUSMD approach also achieve good resolution of the bubble rising flow behavior to a certain range of Reynolds and Bond numbers.

volume fraction;

nˆ = ∇ c˜ =

∂ α˜ . ∂ xi

(3)

Since ∇ c˜ is nonzero only in the interface transition region, so is the surface volume force. The curvature, κ is calculated form

  κ = − ∇ · nˆ .

(4)

Thus (4) is rewritten as:



nˆ =

∇α |∇α|



or

∂α ∂ xi nˆ =  ∂α ∂α ∂ xk ∂ xk

(5)

The local volume fraction of gas obeys

∂α + uˆ · ∇α = 0 ∂t

(6)

This volume fraction equals 1 and 0 in cells represents pure gas and pure liquid, respectively. Also the volume fraction takes continuous values in cells belonging to the transition region. p is the pressure and the mixture density is defined as ρ = ρv α + ρl (1 − α ).ρ v is the density of the gas phase and ρ l is the density of the liquid phase. The mixture viscosity is defined as μ = μv α + μl (1 − α ).μv is the density of the gas phase and μl is the density of the liquid phase. We non-dimensionalize the equations by introducing dimensionless characteristic variables [38] as follows:

x ∗ x = , uˆ = d0 ∗

uˆ



gd0



∗

,t =

g t, ρ ∗ = d0

ρ ∗ p ,p = , ρl ρl gd0

μ gˆ μ∗ = , κ ∗ = d0 κ , gˆ ∗ = , μl g where d0 is the diameter of a bubble. Thus, we may re-express the momentum equations as

2. Governing equations 2.1. Governing equations The problem of the two-phase flows studied in this paper can be described by a time-dependent incompressible Navier-Stokes model with the assumption of two component immiscible fluids and surface tension effects. The mass conservation for the whole domain under the incompressibility condition may be expressed as:

∇ · uˆ = 0.

σ is the surface tension coefficient, and nˆ is the gradients of the

(1)

The momentum conservation of Navier-Stokes equations takes the form

   ∂ (ρ uˆ ) + ∇ · ρ uˆ 2 = −∇ p + ∇ · μ ∇ uˆ + ∇ T uˆ + σ κ nˆ + ρ gˆ, (2) ∂t where ρ is the fluid density, uˆ is the fluid velocity, p is the pressure, μ is the fluid viscosity, σ is the surface tension coefficient, κ is the interface curvature, nˆ is the unit normal vector to the interface, gˆ is the gravitational acceleration. Note that the surface tension term is a singular term that only comes into effect on the interface between the two fluids. Subsequently, the continuum surface force (CSF) model in [29] is used to model the surface volume force. Here, the ith component of the surface volume force is expressed as a surface integral over a mesh cell. Meanwhile κ is the interface curvature, and

  ∂ ρ ∗ uˆ∗ + ∇ · ρ ∗ uˆ∗ uˆ∗ ∂t∗    1 1 = −∇ p∗ + ∇ · μ∗ ∇ uˆ∗ + ∇ T uˆ∗ + κ ∗ nˆ + ρ ∗ gˆ∗ , (7) Ar

Bo

The non-dimensional Archimedes and Bond numbers used here are thus defined as:

Ar =

ρl gd02 ρl g1/2 d0 3/2 and Bo = . μl σ

And the non-dimensionalized continuity equation and the advection equation of gas volume of fraction is not changed without any extra term here. In addition, The Archimedes number will be replaced by Reynolds number in the calculations of the dam break problem and the Rayleigh-Taylor stability problem asThe non-dimensional Archimedes and Bond numbers used here are thus defined as:

Re =

ρl uˆL μl

Where L is the characteristic length in the considered flow problems. In the numerical formulation of Eqs. (1) and (7), we suggested that the including of extra artificial pressure time derivative terms in momentum equations in addition to mass and volume of fraction convective equations in the current preconditioning matrix,

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

which can modify eigenvalues of the Jacobian matrix to be independent of density and volume fractions. It can be regarded as a modified form of Chorin’s artificial compressibility [24] method for the multi-fluid flow problems. This proposed preconditioning approach reduces the possible eigenvalues disparity from multifluid models to enhance numerical stability and make the single phase predecessor code easier to solve the multi-fluid problems. According to pre-conditioning method proposed by in [24,27,30] to rewrite governing equation and the superscript ∗ is omitted for the following equations as:

⎡

⎤

0 1 ∂V ∂U

κ nx ⎥ ⎢ρ gx + Bo P + + = 0, (8) (E − Ev )k nk Sk − ⎣ 1 ρ gy + Bo κ ny ⎦ ∂τ ∂t 0

formula:

∂ U 1.5U n+1 − 2U n + 0.5U n−1 = ∂t t

⎤

1 ⎢ρ u⎥

⎡

⎤

0 p n ⎥ ⎢ Ek · nk = uˆk ⎣ ⎦ + ⎣ x ⎦ ρv pny α k 0 k

P

∂V V n+1,m+1 − V n+1,m =P . ∂τ τ

1.5U n+1,m − 2U n,m + 0.5U n−1,m V n+1,m + τ t 

n+1,m+1 =− (E − Ev )k nk Sk − H  n+1,m+1 

P

= Res V

(9)

The preconditioning matrix P is following the modification of the Choi-Merkle preconditioner [36] as defined as:

⎡ ⎤

⎡1

1 β 1 ⎢0⎥ ∂U ⎢ 0 P = ⎣ ⎦ [1 , 0 , 0 , 0 ] + =⎣ 0 β 0 ∂V α β

α

0

ρ 0 0

0 0

ρ 0

⎤

0 u(ρv − ρl )⎥ . v ( ρ v − ρl ) ⎦ 1

(13)

  m   i, j 3 i, j,k −1 ∂ U I+ P Vi,mj + Pi,−1j m A+ S 1 Vi, j i+ , j τ ∂ V i, j 2 t 2  +   −   −  − A S 1 Vi−1, j − A S 1 Vi, j + A S 1 Vi+1, j i− , j i− , j i+ , j 2   2  +   −  2 + B+ S  V − B S  V − B S 1 1 1 Vi, j i, j i, j−1 i, j+ i, j− i, j− 2 2 2 m   + B− S 1 Vi, j+1 i, j+ 2   3Ui,mj − 4Ui,nj + Ui,nj−1 m m = Pi,−1 Re s −  = m , (14) i, j j i, j i, j 2t

where A = A+ + A− . A+ and A− are computed based on the eigenvalues of the matrix A and the related right-eigenvector matrix T and its inverse T −1 as

(10) Next, the eigenvalues of the Jacobian matrix P −1 ( ∂∂VE · n ) = A are



.

The governing Eq. (11) can be further linearized as



nˆ = nx iˆ + ny jˆ,



(12)

The resulting equations can be obtained by substituting (9) and (10) into Eq. (4)

k

uˆk = u · nk = unx + vny ,

(11)

where n is the index of the physical-time level and m is the index of the artificial-time level. The artificial-time term is discretized by the following implicit Euler finite-difference formula:

where P is the preconditioning matrix,  is the cell area, τ is the pseudo time step and t is the physical time step, nˆ is an outwardpointing normal vector to face k, and Sk is the area of face k, the conservative variable vector is U = [0, ρ u, ρv, α ]T , and the primitive variable vector is V = [ p, u, v, α ]T .u, and v are the Cartesian velocity components. The viscous flux Ev is non-dimensionalized with the Reynolds number for a Newtonian fluid. The inviscid flux vector Ek · n is given as

⎡

3

uˆ, uˆ, uˆ ± (uˆ ) + β , where β is an artificial compressibility parameter, which is suggested as 1 or the scale of (uˆ )2 for the bubble rising flow calculations in the non-dimensionalized governing equations. 2

2.2. Numerical methodology First, to perform the time discretization of the governing equations, the physical-time term is discretized using a second-order, three-point, backward difference formula and a sub-iteration in pseudo-time are used to converge the equations at each physical time step. The artificial-time term is discretized by a first-order implicit Euler finite-difference formula. The evolution of the residual of the primitive variable V at each grid cell (i, j) can be factorized using the LU decomposition method for the pseudo time iteration. Here, the flux extrapolation of Ein+1/2 can split into convective splitting similar to AUSM type splitting and the interface pressure L R flux. Also, the interface primitive variables V1+1 /2 , Vi+1/2 can be determined through an MUSCL+THINC type extrapolation process to achieve the third-order spatial accuracy and sharp capturing of the interface evolutions of bubble rising as shown in the next section. The details are presented in the sub-sections that follow. 2.3. Time discretization In the dual-time stepping formulation, the physical-time term is discretized using a second-order, three-point, backward difference

A± = T ± T −1

(15)

Here, is the diagonal matrix of eigenvalues obtained from (7) as λ1 , λ2 , λ5 = uˆ,λ3 = uˆ + c, λ4 = uˆ − c with the so-called numerical speed of sound like

c=

  2 uˆ

+β

The splitting of the B matrix (associated with the other coordinate direction) is performed similarly. The above Eq. (12) can be factorized using the LU decomposition method for the pseudo time iteration. Res(V n+1,m+1 ) is the unsteady residual vector, Vm is the spatial difference V m i+1, j,k − V m i, j,k . τ is chosen as the local pseudo time step which is determined by the largest eigenvalue of the preconditioning system of governing equations for each grid cell. The diffusion terms are evaluated by the standard central differencing scheme. Furthermore, Eq. (12) is advanced in pseudotime using a LU-type factorization:



m (D − L )D−1 (D + U ) Vi,mj = m , i, j

where



L = Pi,−1 A+ i−1, j Si−1 , j + B+ i, j−1 S j U =

Pi,−1 j

2



A

−



D=

P −1

i+1, j Si+12 , j

∂U ∂V



+ B− i, j+1 S

(16)

 1 i, j− 2 i, j+

1 2



    3  + P −1 A+ i, j − A− i, j SI + P −1 B+ i, j − B− i, j SJ 2 t

 i, j

4

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

Then the primitive variables at the new pseudo-time level m + 1 is then updated. We can sub three steps to solve the Eq. (16) into

 MR −(4) (MR ) =

−0.25αR (MR − 1 ) + (1 − αR )M−(1) (MR ) MR +(1) (MR ) 2

(D − L )Vi,∗j == m i, j Vi,mj

(D + U )

== 

Vi,mj +1

Vi,mj

+ Vi,mj

==

|M R | < 1 otherwise (24)

Vi,∗j (17)

Where

αL =

2.4. Numerical flux

2 ( p L / ρL ) , p L / ρL + p R / ρR

αR =

2 ( p R / ρR ) p L / ρL + p R / ρR

and Next, a modified form of a low-diffusion flux vector splitting schemes, similar to Liou [25], Chang and Liou [26], Neaves and Edwards [27], Niu and Edwards [30] and simplified for incompressible flows, is presented here. Among other properties, AUSM+ exactly captures a grid-aligned stationary contact discontinuity and contains a mechanism for preserving pressure-velocity coupling at low speeds. The baseline compressible formulation is extended to incompressible flows simply by allowing the physical sound speed to approach infinity. Pressure and velocity diffusion mechanisms usually embedded within the preconditioning AUSM+ and LDFSS as shown in [26–27,30] are not necessary in the current formulation of pressure splitting based on artificial compressibility here. An AUSM+ type interface flux is split into convective and pressure contributions can be simplfied as

1  i+ 1/2 = Ei+1/2 · n uˆ1/2 (L + R ) − 2

    uˆ1/2 (R − l ) n i+1/2 , (18)

where

⎡

⎤

L = ⎣

ρv⎦ α L

1 ⎢ρ u⎥



⎤

⎡

1 ⎢ρ u⎥

and R = ⎣

ρv⎦ α R

And

uˆ1/2



1/2

== c1/2 M1/2

(19)

where the L/R notation denotes left and right fluid-property states at the interface. Then we can define the Mach number for left and right states:

ML/R =

uˆL/R c1/2

(20)

with the velocity component normal to the interface given as

uˆL/R = nx1/2 uL/R + ny1/2 vL/R ,

c1/2 =

uˆ21/2 + β .

(25)

The pressure component of the interface flux is expressed as:

⎡ ⎤ 0

⎢n ⎥  i+1/2 = ⎣ x ⎦ p1/2 , E p,i+1/2 · n n

(26)

y

0 with the interface pressure in the cell interface (i + 1/2, j) defined as:

p1/2 =

1 [ pL + pR ], 2

(27)

Interface compression For the spatial discretization of Fi+1/2 , the primitive variables on the cell interface are determined through the third-order spatial accurate type MUSCL with compressive limiters [31]. In addition, the MUSCL +THINC modified under the concept of the BVD principle [32] along utilizing the tangent of hyperbola for interface capturing (THINC) technique, to determine the right and left interfacial volume of fraction (VOF) of cell interfaces. The BVD idea is to reconstruct the solution functions so that the jumps at cell boundaries are minimized, which effectively reduces the numerical dissipation in the resulting schemes. The BVD reconstruction for each cell i is to prepare two piece-wisely reconstructed interpolation functions, α i < 1 > and α i < 2 > of the solution function for each cell. Without losing generality, we assume that α i < 1 > is a higher-order polynomial-based interpolation, while α i < 2 > might be of low-order but has better monotonicity. The concept of the BVD principle says that two groups of the functions as α i < 1 > and α i < 2 > of the VOF for each cell can be the MUSCL scheme or the THINC scheme. For widely used, we assume that α i < 1 > are values from higher order interpolation as the MUSCL scheme, while α i < 2 > might be of low order but has better monotonicity as the THINC scheme. After defining α i < p > (x) and α i < q > (x), the p and q can be either 1 or 2, so that the boundary variation (BV) can appear as



(21)

Unless indicated differently above, the “1/2” notation represents an arithmetic average of left and right states:

[ ]1/2 =

1 (M ± |M|) 2







BV (α )i+1/2 = αi

xi+1/2 − αi+1 xi+1/2

The ‘numerical speed of sound’ is given as:



M±(1) (M ) =

1 ( [ ]L + [ ]R ). 2

 

Then the minimum BV can be found after employing the following judging formula to uniquely determine the interfacial values as ⎧      < 1>

xi+1/2 − αi+1 (xi+1/2 ))(αi−1

xi−1/2 ⎪ ⎨αi (x ), i f (αi  αi

(x ) = . −αi (xi−1/2 )) < 0 ⎪ ⎩ < 2> αi (x ) otherwise

The interface flux M1/2 is defined as

M1/2 = M ( 4 ) ( ML ) + M ( 4 ) ( MR ) , +

With +

ML ( 4 ) ( ML ) =

−



0.25αL (ML + 1 ) + (1 − αL )M+(1) (ML ) ML + ( 1 ) ( ML ) 2

(29) (22)

Then we can calculate the new value of interfacial left side

αiL+1/2 at xi+1/2 and the value of interfacial right side αiR+1/2 at

|ML | < 1 otherwise (23)

and

(28)

xi−1/2 by

    αiL+1/2 = αi

xi+1/2 and αiR−1/2 = αi

xi−1/2 .

(30)

The concept of BVD is expected to reduce the dissipating errors of the residue term of the computed numerical flux for the convective Euler function and the volume-fraction function. For

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

5

smooth solutions, the BVD naturally uses the value at interfaces from the highest order interpolation method because it is better to find an interface value approximated to the continuous exact solution. For discontinuous solutions, purchasing higher order interpolation method does not necessarily lead to reduce the dissipating error at cell interfaces. By minimizing the jumps at cell interface, the concept of BVD chooses the value that is able to represent a jump within the cell where a discontinuity exists. In this idea, the THINC scheme is shown a better choice. We arrange the simple way as in [32] to perform the concept of BVD in following. The values of interface are decided from the MUSCL scheme interpolation under the concept of the BVD principle follows HINC If(αi+1 − αi )(αi − αi−1 ) > 0 and T BVi,Tmin < T BVi,muscl min

αi (x )BV D = αi (x )T HINC , Otherwise

Fig. 1. the patterns of the interfaces distributed across grid points.

αi (x )BV D = αi (x )MUSCL ,

(31)

Where the minimum value of total boundary variation (TBV), for interpolation methods P = MUSCL or THINC is defined as

  L,MUSCL α , − αiR,P | + |αiL,P − αiR,MUSCL i−1/2 −1/2 +1/2 +1/2  L,T HINC  HINC  α − αiR,P | + |αiL,P − αiR,T , i−1/2 −1/2 +1/2 +1/2  L,MUSCL  R,P L,P R,T HINC  α − αi−1/2 | + |αi+1/2 − αi+1/2 , i−1/2  L,T HINC  α − α R,P | + |α L,P − α R,MUSCL  .

T BVi,Pmin = min

i−1/2

i−1/2

i+1/2

i+1/2

(32)

In the THINC method, the tangent hyperbola function is used as a model function for a discontinuous volume fraction within a mesh cell. In principle, however, any model function α (x) that can connect two states αi−1 and αi+1 in a monotone, compact fashion, such that

αi =

1 x

"

xi−1/2

xi−1/2

α (x )dx,

(33)

can be used. The THINC model proposed in [13–15] is considered, the function of the volume fraction is represented as:

αi (x ) =

αmax 2



1 + γ tanh

  x − xi−1/2 η



xi+1/2 − xi−1/2

− x˜i

,

(34)

where η is chosen as 3.5 in order to localize the discontinuity to be approximately in one cell. One can determine the left / right state interpolations by considering the integrated average value of the flux of material crossing a cell boundary, rather than just the end-states. This leads to another representation for left and right states, defined as follows using THINC-M:

"

1 αi (x )dx, ui+1/2 ≥ 0, ui+1/2 t xi+1/2 " xi+1/2 −ui+1/2 t 1 = αi+1 (x )dx, ui+1/2 < 0. ui+1/2 t xi+1/2

αL,i+1/2 = − αR,i+1/2

xi+1/2 −ui+1/2 t

αR,i−1/2

(35)

D+ =



1

η νi+1/2 + ε −

      ln cosh ηνi−1/2 + C sinh ηνi−1/2

η νi−1/2 + ε    (αi − αmin + ε ) B = exp γ η 2 − 1 , ε = 1 × 1012 (αmax + ε ) C = (B/ cosh (η ) − 1 )/ tanh (η )   u i+1/2 · n i+1/2 t  , xxi+1/2 = xc,i+1 − xc,i , vi+1/2 =  (37)   xi+1/2 · n i+1/2  and

αmin = min(αi+1 , αi−1 ) αmax = max(αi+1 , αi−1 ) − αmin γ = sgn(αi+1 − αi−1 ).

(38)

In the above formula, an accurate estimation of void fraction is a key factor to achieve precise evaluation of surface tension effect. Here, it is know that the interface is located over the two neighboring grids by assuming the α of the interface is 0.5. Six modes of the interface distributing over grids can be organized in Fig. 1. As shown in Fig. 2, the volume of fraction between two adjacent grid points usually can be determined by Eqs. (44) and (45) as

 1 αi, j + αi−1, j + αi, j−1 + αi−1, j−1 4  1 = αi, j + αi, j−1 + αi+1, j−1 + αi+1, j 4  1 = αi, j + αi+1, j + αi+1, j+1 + αi, j+1 4  1 = αi, j + αi, j+1 + αi−1, j+1 + αi−1, j , 4

αi−1/2, j−1/2 = αi+1/2, j−1/2

αi−1/2, j+1/2 And

(36)

αi, j+1/2 αi−1/2, j

    ln cosh ηνi+1/2

  tanh (η ) + C sinh ηνi+1/2 1 + C tanh (η )

1

 1 αi, j + αi, j−1 2  1 = αi, j + αi, j+1 2  1 = αi, j + αi−1, j 2  1 = αi, j + αi+1, j . 2

(39)

αi, j−1/2 =

1 = αmin + αmax (1 − γ D+ ) 2 1 = αmin + αmax (1 + γ D− ), 2

with



αi+1/2, j+1/2

Generalizing to an arbitrary cell interface, we have:

αL,i+1/2

D− =

αi+1/2, j 

(40)

Taking diagram (a) in Fig. 1 for example, the curvature of the VOF line as diagram (a) has two possibilities as shown in Fig. 3, one is the curvature opening upward, another is the curvature

6

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

Fig. 2. Computational domain for surface tension.

Fig. 3. Diagrams of the curvature of interface.

opening downward. First, we defined αA = αB = αD = 0.5.We could get A point from αi−1/2, j+1/2 and αi−1/2, j−1/2 by Eq. (41) as

yA = yi−1/2, j+1/2



−

yi−1/2, j+1/2 − yi−1/2, j−1/2





 αi−1/2, j+1/2 − αA  .

αi−1/2, j+1/2 − αi−1/2, j−1/2

(41)

Similarly, B point get from αi+1/2, j+1/2 and αi+1/2, j−1/2 by Eq. (42) as

 yB = yi+1/2, j+1/2 −

yi+1/2, j+1/2 − yi+1/2, j−1/2





 αi+1/2, j+1/2 − αB  .

αi+1/2, j+1/2 − αi+1/2, j−1/2

(42) We used the same way to get D point by Eq. (43)

 yD = yi, j+1/2 −

yi, j+1/2 − yi, j−1/2





 αi, j+1/2 − αD  .

αi, j+1/2 − αi, j−1/2

(43)

The point C is the midpoint of line AB. The position is compared between point C and point D. The point D clould be on the upper position of point C or on the lower position of point C as shown in Fig. 3. For diagram(a)-(b), we defined αA1 = αB1 = αD = 0.5. Then, we calculated yA1 and xB1 by Eqs. (49) and (50) as

 yA1 = yi−1/2, j+1/2 −

yi−1/2, j+1/2 − yi−1/2, j−1/2





 αi−1/2, j+1/2 − αA1  ,

αi−1/2, j+1/2 − αi−1/2, j−1/2

(44)

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

7

Fig. 4. (a) The simulated vector surface tension by J. U. Brackbill et al. [31] (b) The current simulated vector of surface tension (Bo=166).

Fig. 5. Comparison of the interface construction of the evolution of density gradients for a Rayleigh-Taylor problem (t = 5.6 s).

 xB1 = xi−1/2, j+1/2 −

xi+1/2, j−1/2 − xi−1/2, j−1/2





 αi+1/2, j−1/2 − αB1  .

αi+1/2, j−1/2 − αi−1/2, j−1/2

Next, we could calculate αF2 as

αF2 = αF1 −

(48)

Finally, we can obtain the position of the point from

(45) From the above, we could further calculate α E and αF1 from

(y − y )(α − αi+1/2, j−1/2 ) αE = αi+1/2, j+1/2 − i+1/2, j+1/2 A1 i+1/2, j+1/2 , (yi+1/2, j+1/2 − yi+1/2, j−1/2 ) (46)



  αA2 − αD yA2 − yB2   yD = yA2 − , αA2 − αB2 where

αA2 = and

αB2 = αF1 = αi+1/2, j+1/2 −

(yF1 − yF 2 )(αF1 − αB1 ) . (yF1 − yB1 )

(xi+1/2, j+1/2 − xB1 )(αi+1/2, j+1/2 − αi−1/2, j+1/2 ) . (xi+1/2, j+1/2 − xi−1/2, j+1/2 ) (47)

(49)

 1 αA1 + αF2 , 2

(50)

 1 αB1 + αi−1/2, j−1/2 . 2

(51)

The point C is shown to be the midpoint of line A2 B2 . The position between point C and point D is computed to know the curvature of interface whose surface normal vector is toward downward or upward.

8

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

rate of instability can be estimated by measuring the displacement of the interface on the wall. According to the linear theory, the viscosity has an important effect on the growth rate of Rayleigh– Taylor instability. The interface between the different density of fluid is described by the function single wavelength perturbations induced at the interface as y(x)=1.95+0.05cos(2π x). As described in the [33], the interface position is fixed and an initial velocity is assigned to the two fluids as

#

u=

v=

Fig. 6. Dimensionless growth rate versus Reynolds numbers.

Based on the above formula, the computed surface tension in the current approach is shown to be consistent with the results computed by Brackbill et al. [31] shown in Fig. 4. Therefore, the procedure demonstrated from Eqs. (44) to (56) is adopted to achieve accurate curvature of volume of fraction. 3. Results and discussions The first test case considers the gravity-driven Rayleigh-Taylor instability as in Daly [33]. The computational domain is a 1 m × 4 m rectangle with a uniform grid of 100 × 400 used. Following the Daly’s work in which the density of the heavy fluid is 1.225 kg/m3 and the density of the light fluid is 0.1694 kg/m3 . The flow assumption is with the two fluids with equal kinematic viscosities without considering surface tension. The pressure field is set to hydrostatic equilibrium initially. The inviscid wall boundary condition of the computation is imposed on all sides of the boundary. This growth

π A y 2L

−



π A y 2L

π A y 2L



sin πLx exp sin

cos

 πx  L

πx L

 −π |y| 

exp

L

 −π |y|  L

 exp

−π |y| L

as y ≥ 0 as y ≤ 0

(52)

 (53)

where tperturbation of wavelength λ = 2 L is the wavelength of the perturbation, A is the amplitude, and y is the mesh spacing in the vertical direction. For the case under consideration, L = 0.02, A = 0.1. Also in the works of Chandrasekhar [34], the dimensionless growth rate n = n∗ (ν /g2 )1/3 is used for the Rayleigh-Taylor instability as a function of the dimensionless wave numberk = k∗ (ν 2 /g)1/3 , where n∗ is the growth rate and k∗ is the wave number and g is the gravity. For a given wavelength of perturbation and Reynolds number, there exists a single positive exponential growth rate for this type of instability. The computed evolutions of the interfaces performed by the MUSCL and MUSCL+THINC are shown in Fig. 5. The computations simulate that heavy fluid sinks gradually from the center of the interface, then, the interface continues to roll up to a mushroom shape. Due to the nonlinearity of the flow characteristics, the resolution of the interface always demonstrates the typical benchmark test of numerical dissipation effects. The computed results show that the MUSCL+THINC method well resolves interfaces with much less numerical dissipations than the MUSCL method does. Table 1 shows the comparison of the percent mass error for the heavier and lighter fluids of the computations done by the MUSCL and MUSCL+THINC schemes. The maximum error occurs just after initialization for all cases and its relatively high value may be due to the process of adjusting the analytic interface profile to a form consistent with the mesh resolution. It is shown that the mass errors produced by the MUSCL+THINC scheme are slightly than by

Fig. 7. Simulation of free interface of a dam break flow problem.

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

9

Fig. 8. Simulation of the free a broken dam flow problem.

Fig. 10. The regime map of experimental rising bubble shape in liquids [4].

Fig. 9. Effect of mesh sizes on the simulated bubble shapes (Bo=243 Ar=15.23) (50 × 100 (Grid 1), 75 × 150(Grid 2), 100 × 200 (Grid 3), 150 × 300 (Grid 4).

Table 1 Estimated mass error for Rayleigh-Taylor instability. Scheme

MUSCL MUSCL+THINC

the pure MUSCL method. A further comparison of the computed growth rate done by the MUSCL+THINC method with the prediction by linear analysis of Chandrasekhar is given in Fig. 6. Numerical results of both schemes show good agreement with the theoretical curve for Reynolds number 50, 100 and 200 where 3/2 1/2

Mass error (light)

Mass error (heavy)

Average

Max

Average

max

0.152% 0.12%

0.1% 0.1%

0.2% 0.18%

0.12% 0.1%

the Reynolds number is defined as Re = λ νg . Kinematic viscosities and 4.54 × 10−5 m2 /s corresponding to these Reynolds numbers are ν = 2.05 × 10−4 , 1.11 × 10−4 , respectively. The second test case without considering surface tension simulates the dam-break problem. In our variation, the sudden collapse of a rectangular column of water onto a planar surface. The water column, with a base length of 0.06 m and a height of 0.12 m,

10

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

Experiments

Simulations

Test

Observed terminal bubble

Modelling

Predicted terminal

Mass

conditions

shapes

conditions

bubbles shapes

error

Bo=32.2

Bo=32.2

Re=55.3

Re=53.25

Bo=115

Bo=115

Re=94.0

Re=89.36

Bo=243

Bo=243

Re=7.77

Re=7.92

Bo=116

Bo=116

Re=7.16

Re=7.28

Bo=116

Bo=116

Re=13.3

Re=12.22

Bo=116

Bo=116

Re=20.4

Re=19.84

Bo=116

Bo=116

Re=42.2

Re=41.16

Bo=115

Bo=116

Re=94.0

Re=89.36

2.30%

2.70%

2.40%

2.88%

2.81%

2.65%

1.70%

1.67%

∗ Fig. 11. Comparison of the terminal bubble shapes in experiments [1] between the predicted under various Reynolds and Bond numbers (Re = Ar × U∞ ).

is initially supported on the right by a vertical plate drawn up rapidly at time t = 0.0 s. The water falls under the influence of gravity (g = 9.81 m/s2) acting vertically downwards. The density of water ρ =10 0 0 kg/m3. The density of air is taken to be 1 kg/m3.This classical dam break problem, has been frequently employed to validate the code for predicting free surface hydrodynamics. The slip boundary conditions were applied on the bottom and sides of the reservoir. The stress boundary conditions were prescribed on the upper open boundary. They may be changed to fixed pressure and zero normal gradients of the velocities. The gravity causes the water column in the left of the reservoir to seek the lowest possible level of potential energy. Thus, the column will collapse and even-

tually come to rest. The initial stages of the flow are dominated by inertia forces increasing as the water comes to rest. On such a large scale, the effect of surface tension forces is unimportant. Corresponding with those used in the experiment carried out by Martin and. Moyce [35], Fig. 7 shows that both the MUSCL and The MUSCL+THINC scheme simulate volume of fraction plots are shown at times t = 0.281 with the meshes are of size 300 × 80. Both schemes achieve the similar resolution of the free surface. Table 2 shows that the percent mass errors for the gas and liquid phases computed by the MUSCL+THINC are smaller than the MUSCL scheme. Fig. 8 shows the time history of the water column height at the left boundary and also the leading water front posi-

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

11

Fig. 12. The predicted bubble shapes at different Archimedes and Bond numbers. Table 2 Estimated mass errors for dam-break problem. Grid number

MUSCL+THINC MUSCL

Mass error (gas)

Mass error (liquid)

Average

Max

Average

Max

0.004% 0.01%

0.09% 0.1%

0.06% 0.08%

0.1% 0.102%

Table 3 Estimated mass errors for the bubble problem. Grid number

150 × 300 100 × 200 75 × 150 50 × 100

Mass error (gas)

Mass error (liquid)

Average

Max

Average

Max

0.004% 0.006% 0.012% 0.015%

0.08% 0.09% 0.02% 0.02%

0.05% 0.08% 0.1% 0.12%

0.01% 0.01% 0.15% 0.15%

tion at the bottom boundary, where the numerical results obtained using the MUSCL+THINC method are shown to be closer to experimental data than the MUSCL method version of the method. Finally, a grid independence study is executed by simulating the bubble rising the under the condition of (Bo = 243, Ar = 15.23, ρl /ρg = 10 and μl /μg = 100) on different uniformed grids. Four different grid points of 50 × 100, 75 × 150, 100 × 200 and 150 × 300 are used in the calculations. Only the MUSCL+THINC scheme is used here. Table 3 shows that the percent mass errors for the gas and liquid phases computed on different grid points. The effects of grid sizes on the prediction of terminal bubble shapes shown in Fig. 9. It is shown that the numerical errors and the predicted shape of the bubble on 100 × 200 grid points is very close to as those on the finer meshes of 150 × 300. It is noted the computation on 100 × 200 grid points achieves the grid independence. Therefore, the grid points of 100 × 200 are used for the following calculations. Here, the comparison of experimental results [4] with numerical simulated terminal bubble shapes [36–39] for a range of Reynolds and Bond number as shown in Fig. 10. This figure has different locations for the transition lines and introduces additional shape regimes to the diagram introduced by Clift et al. [3]. The eight different shape regimes can be denoted into oblate ellipsoid (OE), Oblate ellipsoidal cap (OEC), oblate ellipsoidal (disk-like and

wobbling) (OED), Spherical (S), Spherical cap with closed steady wake(SCC), spherical cap with open, unsteady wake (SCO), skirted with smooth, steady skirt (SKS), skirted with wavy, unsteady skirt (SKW). Our results are represented by the dot lines as shown in 10. It is shown that the current computations successfully achieve the steady states of the shape of bubbles under the flow conditions with intermediate Reynolds and Bond numbers (10< Re < 100 and 1 < Bo < 100) in which various bubble shapes as oblate ellipsoid, disk-like, oblate ellipsoidal cap, spherical cap with closed steady wake (SCC) and with open, unsteady wake (SCO), skirted spherical cap with smooth, steady skirt and with unsteady wake have been simulated. During this regime, the bubbles rise steadily in a straight path. With a further increase of the Reynolds number (100 < Re < 500) and the Bond number (100 < Bo < 500), the calculations require long iteration time to reach the steadystate in every physical time step in the current computing code. It is noted that the bubble shape may become toroidal due to the unstedy flow regime. As the size of the bubble size increases, the turbulent wakes develop behind the bubble and lead to more unsteady bubble rising motions. At the same time, bubbles oscillating, even breaking up or coalescing are aslo shown. From Figs. 11 and 12, comparisons of the terminal bubble shapes with the data in the works of Bhaga and Weber [4] is depicted. It is shown that the evolution from the flat bubble to the concave shape during the flow regime of the Bond numbers (20 < Bo < 400) and the medium Reynolds number (10< Re < 100). The indentation of the bubble shape is also clearly simulated. The predicted terminal bubble shapes are show to agree well with experimental data [4] for most of test cases. From the results presented here, it can be noted that the current preconditioning incompressible hyperbolic type method is robust enough to reasonably predict the various bubble shapes under the considered flow regimes. Subsequently, the predicted bubble shapes in a wide-ranging of Archimedes and Bond numbers under the flow condition of the ratios of density and viscosity as ρ l /ρ g = 10 and μl /μg = 10 are shown in Fig. 12. While with a small increase in Archimedes number (Ar = 10), the bubble shape is shown to be thinner for the lower Bond number, where the surface tension of the bubble is higher. On the other hand, the bubble shapes become more concave as Bond number increases. However, while the Archimedes number (Ar = 20) is fixed, the bubbles shapes are evolving from flat to concave as Bond number increases (5
12

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269

Fig. 15. The evolution of bubble rising. a) Ar=10, Bo=10;b)Ar=10,Bo=100; c) Ar=50,Bo=100. Fig. 13. The rising velocity vs. time at the different Archimedes number. (Bo=116).

a further increase in the Bond number, the terminal bubble velocity would become slower. Finally, the evolutions of bubble rising demonstrated in different conditions are shown in Fig. 15. It is shown that the shape of bubble gets curved as the Bond number increases, then the bubble evolves to be closer to be horseshoeshaped as the Archimedes number is getting larger. 4. Conclusions A pre-conditioning hyperbolic incompressible multi-fluid flow model based on MUSCL+THINC sharp interface technique is first developed to simulate the Raleigh-Taylor instability, dam breaking and bubble rising flow problems. The grid independence study is performed here. The simulated results were compared with the validated data and satisfactory comparisons. To sum up numerical investigations, we can conclude that

Fig. 14. The rising velocity vs. time at different Bond numbers. (Ar=20).

number (50
1. The current pre-conditioning incompressible multi-fluid flow model combined with the MUSCL+THINC interface compression well resolves the interfaces with less dissipative effects than the original MUSCL-type preconditiong approach does. 2. For the bubble rising flow problems, the MUSCL+THINC type preconditiong AUSMD approach has achieved satisfactory simulations under the flow regime with the intermediate Reynolds and Bond numbers (10< Re < 100 and 1 < Bo < 100). Acknowledgment The author wishes to acknowledge the finical support sponsored by the Ministry of Science and Technology of Taiwan, R.O.C. under contract MOST 106-2221-E-032 -034 -MY3. References [1] Davies RM, Taylor GI. The mechanism of large bubbles rising through extended liquids and through liquids in tubes. In: Proceedings of the Royal Society of London. Series A, Mathematical and Physical Sciences, 200; 1950. p. 375–90. [2] Hartunian RA, Sears WR. On the instability of small gas bubbles moving uniformly in various liquids. J Fluid Mech 1957;3:27–47. [3] Clift R, Grace JR, Weber ME. Bubbles, drops, and particles. New York: Academic Press; 1978. [4] Bhaga D, Weber ME. Bubbles in viscous liquids: shapes, wakes and velocities. J Fluid Mech 1981;105:61–85. [5] Moore DW. The rise of a gas bubble in viscous liquid. J Fluid Mech 1959;6: 113–130.

Y.-Y. Niu and C.-H. Weng / Computers and Fluids 192 (2019) 104269 [6] Taylor TD, Acrivos A. On the deformation and drag of a falling viscous drop at low Reynolds number. J Fluid Mech 1964;18:466–76. [7] Stewart HB, Wendroff B. Two-phase flow: models and methods. J Comput Phys 1984;56:363–409. [8] Hirt CW, Nichols BD. Volume of fluid (VOF) method for the dynamics of free boundaries. J Comput Phys 1981;39:201–25. [9] Ubbink O, Issa RI. A method for capturing sharp fluid interfaces on arbitrary meshes. J Comput Phys 1999;153:26–50. [10] Puckett EG, Almgren AS, Bell JB, Marcus DL, Rider WJ. A high-order projection method for tracking fluid interfaces in variable density incompressible flows. J Comput Phys 1997;130:269–82. [11] Pilliot JE, Puckett EG. Second order accurate volume-of-fluid algorithms for tracking material interfaces. J Comput Phys 2004;199:465–502. [12] Pan D, Chang CH. The capturing of free surfaces in incompressible multi-fluid flows. Int J Numer Methods Fluids 20 0 0;33:203–22. [13] Xiao F, Honma Y, Kono T. A simple algebraic interface capturing scheme using hyperbolic tangent function. International Journal of Numerical Methods in Fluids 2005;48:1023–40. [14] Yokoi K. Efficient implementation of THINC scheme: a simple and practical smoothed VOF algorithm. J Comput Phys 2007;226:1985–2002. [15] Cassidy DA, Edwards JR, Tian M. An investigation of interface sharpening schemes for multiphase mixture flows. J Comput Phys 2009;228:5628–49. [16] Niu YY. Simple conservative flux splitting for multi-component flow calculations. Numer. Heat Transfer, Part B 20 0 0;38:203–22. [17] LeVeque RJ, Shyue KM. Two-dimensional front tracking based on high resolution wave propagation methods. J Comput Phys 1996;123:354–68. [18] Sussman M, Smereka P, Osher S. A level set approach for computing solutions of incompressible two-phase flows. J Comput Phys 1994;114:146–59. [19] Sethian JA, Smereka P. Level set methods for fluid interfaces. Annu Rev Fluid Mech 2003;35:341–72. [20] Olsson E, Kreiss G. A conservative level set method for two phase flow. J Comput Phys 2005;210:225–46. [21] Sussman M, Puckett EG. A coupled level set and volume of fluid method for computing 3D and axisymmetric incompressible two-phase flows. J Comput Phys 20 0 0;162:301–37. [22] Son G. Efficient implementation of a coupled level-set and volume-of-fluid method for three-dimensional incompressible two-phase flows, numerical heat transfer, part B. Fundamentals 2003;43:549–65.

13

[23] Fedkiw RP, Aslam T, Merriman B, Osher S. A non-oscillatory Eulerian approach to interfaces in multimaterial flows (the ghost fluid method). J Comput Phys 1999;152:457–92. [24] Chorin AJ. A numerical method for solving incompressible viscous flow problems. J Comput Phys 1967;2:12–26. [25] Liou MS. A sequel to AUSM, part II, AUSM+-UP for all speeds. J Comput Phys 2006;214:137–70. [26] Chang CH, Liou MS. A new approach to the simulation of compressible multifluid flows with AUSM+Up scheme. J Comput Phys 2008;227:840–73. [27] Neaves MD, Edwards JR. All-speed time-accurate underwater projectile calculations using a preconditioning algorithm. J Fluids Eng 2006;128:284–96. [28] Unverdi SO, Tryggvason G. A front-tracking method for viscous, incompressible, multi-fluid flows. J Comput Phys 1992;100:25–37. [29] Brackbill JU, Kothe DB, Zemach C. A continuum method for modeling surface tension. J Comput Phys 1992;100:335–54. [30] Niu YY, Edwards JR. A simple incompressible flux splitting for sharp free surface capturing. Int J Numer Methods Fluids 2012;69:1661–78. [31] Leer BV. Towards the ultimate conservative difference, V: a second order sequel to Godunov’s method. J Comput Phys 1979;32:179–86. [32] Sun Z, Inaba S, Xiao F. Boundary Variation Diminishing (BVD) reconstruction: a new approach to improve Godunov schemes. J Comput Phys 2016;322:309–25. [33] Daly BJ. Numerical study of two fluid Rayleigh–Taylor instability. Phys Fluid 1967;10:297–307. [34] Chandrasekhar S. Hydrodynamic and hydromagnetic stability. Oxford University Press; 1961. [35] Martin JC, Moyce WJ. An experimental study of the collapse of liquid columns on a rigid horizontal plane. Philos Trans A 1952;244:312–24. [36] Turkel E. Preconditioning techniques in computational fluid dynamics. Annu Rev Fluid Mech 1999;31:385–416. [37] Hysing S, Turek S, Kuzmin D, Parolini N, Burman E, Ganesan S, Tobiska L. Quantitative benchmark computations of two-dimensional bubble dynamics. Int J Numer Methods Fluids 2009;60:1259–88. [38] Hua J, Lou J. Numerical simulation of bubble rising in viscous liquid. J Comput Phys 2007;222:769–95. [39] Weiss JM, Smith WA. Preconditioning applied to variable and constant density time-accurate flows on unstructured meshes. AIAA 1994:1994–2209.