ACTA MECHANICA SINICA, Vol.20, No.4, August 2004 The Chinese Society of Theoretical and Applied Mechanics Chinese Journal of Mechanics Press, Beijing, China Allerton Press, INC., New York, U.S.A.

ISSN 0567-7718

IMPROVED SUBGRID SCALE MODEL FOR DENSE TURBULENT S O L I D - L I Q U I D T W O - P H A S E FLOWS* T A N G Xuelin ( f N ~ )

1,t

QIAN Zhongdong ( ~ , ~ ) 2

WU Yulin ( ~ , $ $ ) 2

1(Department of Hydraulic and Hydropower Engineering, Tsinghua University, Beijin9 100084, China) 2(Department of Thermal Engineering, Tsinghua University, Beijing 100084, China) ABSTRACT: The dense solid-phase governing equations for two-phase flows are obtained by using the kinetic theory of gas molecules. Assuming that the solid-phase velocity distributions obey the Maxwell equations, the collision term for particles under dense two-phase flow conditions is also derived. In comparison with the governing equations of a dilute two-phase flow, the solid-particle's governing equations are developed for a dense turbulent solid-liquid flow by adopting some relevant terms from the dilute two-phase governing equations. Based on Cauchy-Helmholtz theorem and Smagorinsky model, a second-order dynamic sub-grid-scale (SGS) model, in which the sub-grid-scale stress is a function of both the strain-rate tensor and the rotation-rate tensor, is proposed to model the two-phase governing equations by applying dimension analyses. Applying the SIMPLEC algorithm and staggering grid system to the two-phase discretized governing equations and employing the slip boundary conditions on the walls, the velocity and pressure fields, and the volumetric concentration are calculated. The simulation results are in a fairly good agreement with experimental data in two operating cases in a conduit with a rectangular cross-section and these comparisons imply that these models are practical. K E Y W O R D S : kinetic theory, turbulent two-phase flow, dynamic sub-grid-scale model, conduit

1 INTRODUCTION So far, most studies concern only dilute solidliquid flows and the continuum theory is applied [1'2] . Assuming that the motions of particles are caused only by the ambient liquid and have no relations with the collisions between particles, the procedures of turbulent models for sing-phase flows may be applied to develop turbulent two-phase models, such as the particle trajectory model [3], the multi-fluid model k-cAp [4'5], the k-c-kp model and the USM model [6], and the LES model [7's]. Despite successful applications of the continuum theory to dilute two-phase flows, it is a fact that the interactions between particles in a dense turbulent solid-liquid two-phase flow are not considered, so other models must be introduced. The dense turbulent two-phase flow is quite different from a dilute one. Except t h a t the particles of a dense flow have the same flow properties as the particles of the di-

lute flow, the collisions between particles can not be neglected; thereby the viscosity and diffusion due to particle collisions should be considered. But based on Boltzmann equation, the kinetic theory can well describe the microscopic interacting properties. Such research techniques on dense solidliquid two-phase flows made no use of the relatively well-developed dilute two-phase flow theory. Ni[9] investigated the sediment-laden flow, and for a dilute two-phase flow. Firstly the particle-phase Boltzmann equation was greatly simplified by ignoring the collisions between particles, and thus the particle-phase velocity distribution function was solved. Finally, the particle-phase flow parameters were obtained; and then for a dense two-phase flow, the collisions between particles were considered and assuming that the collision t e r m is linear to the difference between the particle velocity distribution function of the dense twophase flow and t h a t of the dilute two-phase flow, the

Received 23 October 2002, revised 19 August 2003 * The project supported by the Nationai Natural Science Foundation of China (50176022) t E-maih [email protected]

Vol.20, No.4

Tang XL et al.: Subgrid Scale Model for Turbulent Two-phase Flows

Boltzmann equation was integrated, the dense velocity distributions were determined and the flow parameters were obtained. Ishii and Shih-I Pail 1~ multiplied each phase Boltzmann equation by property parameters and integrating over the velocity space, the continumn and moment equations for the solid-phase two-phase flow were obtained, but the collision concerns a complex integral and can not be integrated analytically, so it can not be used for engineers. Liu Dayou[121 also did the same as Ishii and Shih-I Pal did, and tried to derive the collision terms, but because of this complexity, only obtained the collision term under the dilute two-phase flow condition. Samuelsberg et al. [13], Mathlesen et al. [14'15] and Benyabia et al. [16] studied the flow parameters about the fluidization bed by almost completely employing the gas motion theory, where the relation between gas temperature and velocity is the same as that between particle temperature and velocity, and there is no actual difference between gas molecules and particles. The predicted results are in good agreement with experimental data. Filippova et al. [lz] connected the twophase flow BoRzmann equation at the meshes with the continuum macro-models and calculated the interrelated parameters, and then directly resolved the two-phase flow Boltzmann equation at the meshes, and then obtained each phases flow parameters. So far, based on Boltzmann equation, there are three methods for the solid phase in a two-phase flow: (1) Boltzmann equation velocity distribution function method; (2) Boltzmann equation integrating method; (3) LBM (Lattice Boltzmann Method). At the same time, the solution methods for the governing equations also developed rapidly. There are three types of CFD: direct numerical simulation (DNS)[ ls], Reynolds averaged Navier-Stokes equation model (RANS) and large-eddy simulation (LES). Since the LES requires less computer time than the DNS and uses simpler and more universal models than the RANS, the LES becomes an important and promising method. LES is a turbulence simulation method, which assumes that the turbulent motion could be decomposed into large-scale (grid-scale) eddies and small-scale (subgrid-scale or SGS) eddies. The difference between the two does not have a significant effect on the evolution of the large-eddies, and the small-eddies is independent of the flow geometry and the boundary conditions. The LES solves the large-scale field directly by a set of filtered governing equations, while the SGS is isotropic and mainly serves as a dissipation, so its effect on the grid scale

355

is modeled. LES was proposed in 1963 by Smagorinsky [19] for meteorological applications and was used in the 1970s by Deardorff[ 2~ for industrial applications ( t h a t is, Smagorinsky model (SM)). Though very successful for turbulent flows, SM has some notable drawbacks[ 2q. Cermano et al. and Lilly [2"~'231 proposed a dynamic subgrid scale model (DSM) to solve the problem associated with the Smagorinsky model by computing directly the model parameter with the resolved scales. In other words, the DSM model parameter can be obtained as a function of the flow domain and time. This model can predict correctly the asymptotic behavior near a wall, and permit the energy to backscatter from small scales to large scales. The DSM has been successfully used in some simple flows, but it still needs improvements for complex flows. Yan Zang et al. and Sandip et al.[24,25J analyzed and improved the DSM. The above models are of eddy-viscosity type and based on SM. Furthermore, there are also other noneddy viscosity models, such as the scale-similarity model put forward by Bardina et al. [26]. In the article, by analogy with the gas molecular kinetic theory[ 27,2s], the dense solid-liquid twophase flow is investigated. Multiplying the solid-phase Boltzmann equation by its characteristic parameters and integrating over the velocity space, the governing equations for the solid-phase two-phase flow are obtained. The collision term between particles is also obtained. By adopting some relevant terms in the dilute two-phase governing equations and the dense particle collision term, the governing equations are developed for the dense solid-laden turbulent flow. And then based on Cauchy-Helmholtz theorem and SM, a second-order dynamic model with two parameters is put forward to solve these dense turbulent two-phase governing equations. This model includes both the strain-rate tensor and the rotation-rate tensor. Besides, the slip boundary conditions on the walls as mentioned in Refs.[29,30] are also used here. Finally experimental data from Wang Zhaoyin and Ning Chien [31'32] are used to check these models. 2 SOLID-PHASE

GOVERNING

EQUA-

TIONS

2.1 M i c r o s c o p i c a n d M a c r o V a r i a b l e s o f Particle Phase Firstly assuming that: (1) the particle-phase has only translation, but not rotation; (2) volume unit d R is much less than the macro space, but much larger than every particle volume, that is to say, it includes

356

ACTA MECHANICA SINICA

numerous particles, therefore, the statistical approach can be applied to the parOcle-phase, where R is vector radius in the geometric space. By employing gas molecular kinetic theory to the solid-liquid two-phase flow, the solid-phase particle velocity distribution f can be determined by Boltzmann equation [27]

Of

can be obtained as

~(pp) 4- O@i(flpUpi) = O

Of

where the particle-phase mass density pp = n m and the macro-velocity Upi = (ci}. W h e n 0 = rnci in the transportation equation, the particle-phase m o m e n t u m equation can be derived

(1)

(of)

~

stands

for the

~

collision terms between particles

~(R, c, O, r, t)f(R, c, O, r, t)dcdOdr

//~f //~

dcdOdr

(2)

where n denotes the particle number per unit volume. 2.2 S o l l d - p h a s e G o v e r n i n g E q u a t i o n s

In order to simplify the Boltzmann equation, it is assumed that particles have the same radius r and t e m p e r a t u r e 0, therefore the fourth and fifth terms on the left side of the Boltzmann equation are zero. Then, multiplying the above equation by one certain characteristic p a r a m e t e r 9 of the particle-phase and integrating over the velocity space, the macrovariables' (statistical averaged variables) transportation equation can be finally written as follows

/ ~[(~S-)ppcq- (~-)pfc]d3cj

(pp%

ps) -

f mcj[( ~Of_ ) p p c +

Opp

Ppgpi + Ox5 (Of~

_

ldac j

\ 0~ / pfcJ

(5)

where U p i l t p j ~--- (CiCj} --

Bpij

1

pp = ~(Bpll ~- Bp22 -~-Bp33) pp is the particle-phase partial pressure, Tpij is the particle-phase viscosity stress tensor, Fi = mgpi, gpi is the particle-phase mass force, and Kronecker delta tensor is defined as

5ij =

Of dcdOdr

=

OXki

O

ppc

f )- pfc for those between particles and water and (_O _-~dR de F molecules; ~ - = c, dt - m ' F is the particle mass force; m is the particle mass. In statistical mechanics, the statistical averaged value of one certain property p a r a m e t e r ~ of particlephase can be written as

(0) =

+

Orpij _

where c is particle velocity; 0 is k-phase molecule or particle temperature; r is particle radius;

(4)

as

Df Of + Of d R Of dc Of dO Dt - at OR d~t- + 7c d[ + ~ ~ + Of dr

2004

(3)

When 9 = rn in the transportation equation and assuming that there are no any change on the particle and molecule, the particle-phase continuum equation

{10

i=j i 7s j

2.2.1 Particle-phase Pressure Model For a dense two-phase flow, the particle-phase pressure is assumed to be proportional to its volumetric concentration pp = Cpp, where p is the total pressure of the two-phase flow, Cp is the particle volumetric concentration and related to the liquid-phase volumetric concentration Cf as Cp + Cf = 1

(6)

where Cp = Pp/Pps, Cf = Pf/Pfs, p p s and pf~ are the pure particle density and the pure water density, respectively. 2.2.2 Collision Integral Model The two terms on Boltzmann equation's right side are collision integral terms. There are fundamental differences between gas molecules and particles: (1) Gas molecules have thermal motions, and their kinetic energy depends directly on their temperature, whereas particles do not move by their own temperature and their t e m p e r a t u r e has no relationship to their kinetic energy; (2) gas molecules are all elastic sphere and have no energy loss during collisions between them, whereas particles have collision energy

Tang XL et al.: Subgrid Scale Model for Turbulent Two-phase Flows

Vol.20, No.4

losses between them; (3) every particle volume is often several orders larger t h a n every gas molecule's, so its volume can not be neglected. 2.2.3 Collision Integral Term between Liquid Molecules and Particles The second term on the Boltzmann equation's right side is the collision integral term between liquid molecules and particles, which is assumed to be the interaction force between liquid molecules and particles from the dilute two-phase flow theory, namely: (The interaction force that the liquid-phase acts on the particle-phase Fpfcj = a -w-/mcj(~7)pfcdac j and the force that the particle-phase acts on the liquidphase is equal to Fpfcj in magnitude, but with an opposite sign, namely, Ffpcj = -Fpfcj) fpfcj

= f mcj (~)pfcd3cj= P--e-P

where r~.p = 18~- 1 + ~s

)

eraged relaxation time, dp is the particle radius, pp is the particle density, > is the liquid viscous coefficient, Rep = lUp - u f l d p / ~ is the particle relative Reynold number, zJ is the liquid kinematic viscosity, Up is the particle velocity, uf is the liquid velocity. 2.2.4 Collision Term between Solid Particles Assuming that an area element dS lies at the point /g with unit normal vector n. When the collision between particle 2 and particle 1 on the positive and negative sides of dS occurs, the line of their centers at the collision passes through dS, and the mean positions of the two particle are at the points R - r h and R + r h , respectively, the velocities of particles 1 and 2 are el and c2, respectively, before the collision, C'1 and C2' after the collision, h is the unit vector drawn from the center of particle 1 to t h a t of particle 2 , the velocity of particle 1 relative to particle 2 is c12 = cl - c2. Considering t h a t the particle is not smooth and elastic, the particle elastic recovery coefficient is assumed to be e. With respect to the property ~b before the collision and the corresponding property ~, after the collision, each such collision causes a transfer quantity ( ~ ' - ~) of the property from the negative side to the positive side of dS. Finally, the vector-flux of ~ can be written as follows

,= f f /

f

where c12 9h > 0 assures a collision between particles 1 and 2 and h - n > 0 is obtained from the definition of h and n. /3 is the polar angle between vector h and the relatively velocity c12 as an axis, ~ is the plane angle, fl and f2 are the velocity distribution functions for particles 1 and 2, respectively. When ~ = m e , the particle peculiar velocity c' and the particle velocity distribution function are b o t h Maxwell distribution, f l and f2 are expanded and only terms of zero order remain, and the particlephase velocity maintains approximately a local equilibrium. Finally, the m o m e n t u m variation per unit time and volume, namely, the collision force between particles can be expressed as Fppcj = f m c j ( 0 f ~ d3cj = V . J \ O t / pc

(7)

is the particle av-

h)k(R- h)

357

= 2(i + e)V(GpptC'2})

(9)

For a dilute two-phase flow, some simplifications are made: for the particle-phase, with no collisions between particles, there are no pressure and viscosity inside the particle-phase, namely the effect of collision is very small, therefore, the integral collision term is ignored, thus there are only pressure and viscosity inside the liquid-phase. By applying the molecule kinetic theory, the particle-phase viscous coefficient can be obtained as 2ppr #p -- 3x/~Vp

(10)

g ~c )

Finally, the m o m e n t u m equation for the particlephase can be expressed as 0 u 0 ~-t(PP Pi) ~- ~ p j ( P p % p j ~ p i ) Orpij

Ppgpi q- ~Xpj q- ~ f c i where PP Tr.p

Fppci = fi).

3 LARGE PHASE

2(1 @ r

O(Cpp) Jr0Zp~ 7-

@ /~ppi

(11) Fpfci

=

EDDY SIMULATION FOR TWOGOVERNING EQUATIONS

3.1 S m a g o r i n s k y M o d e l ( S M ) a n d D y n a m i c Smagorinsky Model (DSM) The liquid-phase governing equations can be written as

c12 9 h > 0

n.h>O

( R + r h ) h sin/3d/3d~dct de2

~(pf)-}(8)

~-~i(pfufi)

(12)

ACTA MECHANICA SINICA

358 0

~7 (Pf~fi) + ~

0

Orfij Pfgfi + ~Xfj +

(Pf*tfJ2"tfi) --

O(pfp) pfsOXfi

(la)

Ffpci

where Ffp~ = -Fpf~i. Most SGS models are based on eddy-viscosity assumptions. The most commonly used model is the SM. Based on the assumption about the balance between the energy production and dissipation, applying dimensional analyses and filtering N-S equations with filter function G, Smagorinsky obtained the expression between the sub-grid-scale stress and the eddy viscosity

Tfij =

-- 2Vft

Sfij

(14)

where, ~ft = csA21Sf I is the SGS eddy viscosity coef1 [O~tfi O~tfj ] ficient, Sfij = ~ l~xj + Oxi J is the resolved strain rate tensor, ISfl ---is the magnitude of the strain-rate tensor, and cs is Smagorinsky constant varying from 0.01 to 0.08. The SM is so simple and well designed that it has been applied to m a n y flow fields with a great success. However, the SM has some shortcomings, e.g., the Smagorinsky constant must be optimized for each flow field and the SM requires a damping function to ensure t h a t the SM vanishes at a solid b o u n d a r y and can not account for the energy backscatter, etc. Although the DSM improves the SM, it is less suitable for highly rotational flows. 3.2 Second-order Model

Dynamic

2004

the partic]e as an absolutely rigid b o d y with respect to the point, and the velocity related to the particle deformation. Therefore, based on the modeling principle: a model must be more accurate if including more flow information, the subgrid-scale stress should include above three types of motions to some degrees. Applying dimensional analyses, Smagorinsky model assumes t h a t the subgrid-scale stress is only proportional to the rate of the strain tensor; Joel H. Ferziger [a3] proposed an eddy model, in which the viscosity is related to the rate of the rotation tensor. In most eddy-viscosity models, the turbulent stress is assumed to depend only on the rate of the strain tensor and not on the rotation tensor. This indicates t h a t highly rotational flows would be less suitable to model using the eddy-viscosity concept. It is believed t h a t the vortex stretching is the dominant mechanism by which the turbulence transfers energy from the grid scale to the SGS (forward scatter), and the subgrid scale stress should be a function of the strain rate tensor and the rotational rate tensor. Based on the above analyses, a second-order dynamic model with two parameters is proposed, in which the sub-gridscale stress is a function of both the strain rate tensor and the rotation tensor. In a dynamic eddy-viscosity model, we define two filtering operators 0 and G (the width of the test filter G is assumed to be larger than that of the grid filter G). So the large-scale quantities are expressed by the following formula

Sub-grid-scale

f(xd

=

fo(x

-

xi)f(xi)dx

'

i

'

(15)

f(x )

=

fo(.

-

xi)f(xi)dx

,

,i

(16)

with Two Parameters

Based on the SM, the second-order dynamic subgrid-scale model with two parameters not only includes the strain-rate tensor but also the rotation tensor, which improves DSM and allows these two model parameters to be a function of time and space. By assuming that the smallest resolved scales of the computation provide information about the energy transfer to the subgrid scales and relating the subgrid scale Reynolds stress to two different sizes of filters, the second-order dynamic sub-grid-scale model with two parameters uses the resolved instantaneous velocity field to compute local values for the two dynamic model parameters, so it can provide the correct nearwall asymptotic behavior and energy backscatter and does not require a wall damping function. According to Cauchy-Helmholtz theorem, the velocity of any p o i n t of an infinitesimal particle is composed of the velocity of the translational movement of the particle center, the rotation velocity of

Let ~ = GG, and applying the filter functions and & to the liquid-phase governing equations t(fif) + 0~7(~7~77) = 0

(17)

0

0

OCfij Pfgfi ~- ~ ~- Ffpci

(18)

0 0 ~ (/sf) + ~ x T x / ( ~ ) = 0

(19)

0

--

0

0

0-7 (fl~7) -}- ~TXj( ~ ) +

0rfo. c, x j

+

-

--

Pfs0Xi (P-~) -}(2o)

Vol.20, No.4

Tang XL et al.: Subgrid Scale Model for Turbulent Two-phase Flows

Two filtering operations weighted by the density are used following the approach used by Favre for compressible flows

N~j = - 2 ( A fifRfikRfkj

359

(31)

-- ~'~2pfRfikRfkj)

Applying a least-square approach

f(xi) = pff(xi)/fif

d

(21)

E E (Lij dc--~ i j

The filtered Navier-Stokes equations can be obtained as follows

- - Cl]~/~iJ --

c2Nij)2 = 0

l z 1, 2

(32)

we have the expressions for the two model parameters as

~

0 (/Sf~tfi) -~- ~--~Xj(RfUfjUfi) --

OTfij

pfgf~ + ~

(LijMij)(NktNk~) - (LijXij)(Mk~Nkl)

0(~f~)

+ Pf~r

cl = ( M ~ j M ~ j ) ( N k d V k ~ ) - ( M ~ j ~ j ) ( M k ~ N k ~ )

pfsOX~

Oqfij

(34)

0 (/Of) - ~- L~--'~ (~f~fi) ^ = 0

0 _~

(24)

0

0(~)

~(P~) + ~(~a~a") -

~o~

~ OTfi j Pfgf~ + ~ + FfPr

where, qfij

orfij Oxj

(25) ^

= Uf~fj

-- ~ f i ~ t f j

^

a n d T f i j = ~tfiu~fj -- ~tfi~tfj

are the subgrid-scale stress and the subtest-scale stress, respectively. Therefore, the resolved stress can be defined as

Lij = ~fTfi j -- P~fij = pfU i pfUj/pf -- pfu i pfuj/~f

(26)

Because the model is not directly related to any filter function, the subgrid-scale stress model and subtest-scale stress model can be written, respectively, as

qfij = -2clz~2 ISf ISfij - 2c2A2RfikRfkj

(27)

Tfij

(28)

where

~fij

= --2cl A

=

eijk

(LijNij)(MklMkt) - (LijMij)(MktNkl)

(23)

Oxj

]SflSfij -- 2C2z~ ~fik~fkj 0gfk

OX~-~' eijk is the permutation ten-

sor, cl and c2 are the dynamic model parameters, A is the characteristic filter width associated with G, 21 is the filter width associated with ~. Let us define the following quantities

Lij = clMij + c2Nij

Finally, a second-order dynamic model with two parameters is obtained, where ( ) represents the average taken over a plane parallel to the wail. By applying this procedure for the liquid-phase equations to the solid-phase equations, the filtered solid-phase governing equations can also be obtained. 4

GEOMETRY OF STRAIGHT AND COORDINATES

CONDUIT

In order to validate the accuracy and credibility of these models, the turbulent flow simulation through the duct is compared and analyzed with the experimental results from Wang Zhaoyin and Ning Chien [31'321. Wang Zhaoyin et al. conducted many sets of experiments in a pressured conduit with a rectangular cross-section using the uniform plastic sand (dp = 0.25mm, Pps z 1.067 x 103 kg/m3). The dimension of the conduit used here is identical to that used by Wang Zhaoyin et al. Its crosssection is 100m high and 180mm wide, its length in the streamwise direction used in the simulation is 4 m as shown in Fig.1. The grid is a mesh of I x J x K = 400 x 46 x 26 (I is the X direction length, Y is the Z direction height, K is the Y direction width).

tlow

~

(29) X

M~ : - 2 ( a

(33)

ls~l~s.~- r~l~I~.~)

(30)

Pig.i Geometry of straight conduit and coordinates

360 5

ACTA MECHANICA SINICA COMPARISONS

MENTAL

BETWEEN

DATA AND

2004

10

EXPERI-

SIMULATION

8

9

Inlet conditions: the velocity is assumed to be uniform.

6

Outlet conditions: the velocity gradients are assumed to be zero.

4

.

A

simulation experimental

--

~

2

Solid conditions: the slip conditions are used here. T h e slip velocity on the wall can be w r i t t e n Ou as: u~,~ = K 0 0 ~ , where rt~ is the normal to the

0 1

2

3

4

u/(m.s -1)

wall, u is the velocity in the streamwise direction, (c) X = 2.0m

/%_ T h e box filter function is employed here. There are two experimental cases employed here to compare to the simulation results. Their m e a n velocities and volumetric concentrations are 3.180 and 0.465, 1.875 m / s and 0.457, respectively. 5.1 C o m p a r i s o n s

between

10 8 6

Velocities

simulation A experiment

\

--

\

4

F r o m Figs.1 and 3, it can be shown t h a t the velocity distributions across the Y - Z section are adjusting themselves and reach relatively steady states

2 0 0

1

3

4

u/(m. s-1)

10 8

2

(d) X = 3.0m

simulation A experiment

--

6

A z~

A

4

10 8

2 0

,

1

2

,

,

4

3

zx expel-lment

6

,I

u/(m.s-l)

4

(a) X = 0.0m

2 0

A~ ,

i

1

10

i

o

2

3

i

i

I

4

u/(nl.s -l)

~ t a LX

8

(e) X = 4.0m simulation experiment

--

6

\ ~ )

A

o 4

Fig.2 Case 1 when the mean velocity is 3.180 and the volumetric concentration is 0.465

2 0 1

2

u/(m.s-1) (b) X = 1.0 m

3

4

(namely the fully developed turbulent two-phase flow) at a distance a b o u t X = 2.0 m from the inlet. These simulation results are in fairly g o o d agreement with experimental data.

Tang XL et al.: Subgrid Scale Model for Turbulent Two-phase Flows

Vol.20, No.4

361

10 simulation a experiment

g

--

g

6

-- simulati2