- Email: [email protected]
Two-fluid large-eddy simulation approach for particle-laden turbulent flows
International Journal of Heat and Mass Transfer 45 (2002) 4753–4759 www.elsevier.com/locate/ijhmt
Two-fluid large-eddy simulation approach for particle-laden turbulent flows R.V.R. Pandya, F. Mashayek
*
Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 842 West Taylor Street, Chicago, IL 60607, USA Received 1 April 2002; received in revised form 18 May 2002
Abstract In this paper, effects of particles on the subgrid scales of turbulence are properly accounted for during the modeling of subgrid scale stresses in the large-eddy simulation (LES) of fluid phase. In doing so, we propose closed filtered kinetic equations for phase space density of the particle. The various moments of these equations give the ‘fluid’ equations which can be considered as the LES equations for the particle phase. The influence of subgrid scales motion on the particles is included in these ‘fluid’ equations. Ó 2002 Elsevier Science Ltd. All rights reserved.
1. Introduction In recent years, large-eddy simulation (LES) has emerged as an important predictive technique for fluid turbulence with growing number of models for statistical properties of small scales or subgrid scales which influence the instantaneous flow properties of the large scales or eddies. The pioneering model for subgrid scale stress tensor is due to Smagorinsky [1], followed by other dynamic and similarity models which are reviewed recently by Meneveau and Katz [2]. LES is also emerging in the field of two-phase turbulent flows for accurate prediction of particle or droplet laden flows [3,4]. In these flows, each particle moves under the influence of fluid forces (e.g. fluid drag force, history force), and the instantaneous fluid flow velocity Ui in the vicinity of the particle governs its trajectory. Here, we consider only the fluid drag force and assume the particle to be spherical with time constant sp ¼ qp d 2 =18l, where l is the fluid viscosity, and qp and d are the particle mass density and diameter, respectively. In this case, the governing Lagrangian equations for the spherical particle’s position, Xi , and velocity, Vi , can be written as *
Corresponding author. Tel.: +1-312-996-5317; fax: +1-312413-0447. E-mail address: [email protected] (F. Mashayek).
dXi ¼ Vi ; dt
ð1:1Þ
dVi ¼ bv ðUi Vi Þ; dt
bv ¼
1 þ 0:15Re0:687 p ; sp
ð1:2Þ
where particle Reynolds number Rep ¼ dqjV Uj=l and q is the fluid density. In this paper, we consider small particles with Rep < 1 and take bv ¼ 1=sp . The presence of particles modifies the flow field which is accounted for through the source term Si in the governing equations for fluid phase velocity ui , written for incompressible fluid as ouj ¼ 0; oxj
oui o op o2 ui þ ½ui uj ¼ þm þ Si : oxi ot oxj oxj oxj
ð1:3Þ
Here, p represents the fluid pressure divided by q at location xi and time t, m is the fluid kinematic viscosity, and Si accounts for the effects of particles on the fluid flow and is given later in this paper. In the case of one-way coupling, Si ¼ 0 and the fluid phase is considered independent of the particle phase. Si has finite values in the case of two-way coupling. Following the work of Germano et al. [5], we denote the grid filtering operation with filter width D on any function f as Z e ðx; x0 Þdx0 ; f~ðxÞ ¼ f ðx0 Þ G ð1:4Þ
0017-9310/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 7 - 9 3 1 0 ( 0 2 ) 0 0 1 9 1 - 6
4754
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Nomenclature C C
d D ~, G G Gjk G0jk G00jk k, K mp np p Rep Si Sij S~ij ; S~ij t Tbi bj ui
subgrid stress model coefficient model coefficient in Eq. (3.9) particle diameter model coefficient in Eq. (3.8) filters Green’s function, Eq. (2.10) Green’s function, Eq. (2.17) Green’s function, Eq. (2.32) subgrid scale kinetic energies per unit mass of fluid mass of the particle total number of particles fluid pressure divided by fluid density particle Reynolds number source term in Eq. (1.3) strain rate tensor filtered strain rate tensors time Lagrangian integral time scale ith component of fluid velocity u
and the test filtering operation as Z f ðxÞ ¼ f ðx0 ÞGðx; x0 Þdx0 ;
ð1:5Þ
e filter has D as filter width. e ¼ GG and G e to the Navier–Stokes Now, application of the filter G equations gives o~ ui o o~ p o2 u~i o þ ½~ ui u~j ¼ þm sðui ; uj Þ þ e Si; oxi ot oxj oxj oxj oxj ð1:6Þ where the unknown subgrid stress tensor ~i u~j sðui ; uj Þ ¼ ug i uj u
ð1:7Þ
and e S i pose closure problems and need to be modelled. While writing (1.6), and also later in this paper, commutation between filtering and derivative operators is used. In case of a single-phase flow, various models for subgrid stress tensor are suggested and are recently reviewed by Meneveau and Katz [2]. In this paper, we use the framework of dynamic localization model, proposed by Ghosal et al. [6], and derive the model for sðui ; uj Þ in case of two-phase flow by incorporating the effects of particles on the subgrid scales. Also, a closed expression for e S i is derived from the kinetic equation for phase space density of the particle phase. In the early studies on the LES of two-phase flows with one-way coupling [7,8], particles are tracked by using the LES velocity field, in the Lagrangian Eqs. (1.1) and (1.2), instead of the instantaneous velocity field ui of
Ui Vi vi W xi Xi
ith component of fluid velocity U along the particle path ith component of particle velocity V phase space variable corresponding to Vi phase space density physical space coordinates ith component of particle position vector X
Greek symbols l fluid viscosity m fluid kinematic viscosity q fluid density qp particle mass density sp particle time constant D, D filter widths Subscript p particle Symbol h i ensemble average
the fluid phase. The effect of subgrid scale velocity field u0i ¼ ui u~i on the particle are not taken into account in these studies. Armenio et al. [9] studied these effects in detail by performing and comparing the results from independent direct numerical simulation and LES. The two-way coupling effects are incorporated in LES by Boivin et al. [10] through the source term Si . Only very recently, the two-way coupling is treated completely by also taking into account the effects of particles in the model for subgrid stress tensor [11]. Yuu et al. [11] modelled subgrid-scale turbulent mass flux of particle by gradient transport using the analogy of the molecular transport. In all of the existing studies on LES, the particle phase is simulated in Lagrangian framework (see the recent review [3] and references cited therein). Rigorous modeling of the effects of particle phase on the subgrid scales and vice versa, still, remains as an important unfinished task. In this paper, an attempt is made for such modeling and, in doing so, ‘fluid’ or Eulerian equations are derived for particles which can be considered as LES equations for the dispersed phase.
2. ‘Fluid’ equations for particle phase Using the phase space density W ðx; v; tÞ, the source term Si can be written as Z mp Si ¼ ðui vi ÞW ðx; v; tÞdv; ð2:1Þ qsp
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
where mp is mass of the particle, x and v are phase space variables corresponding to X and V, and W is governed by [3,12,13] o o o W ðx; v; tÞ þ ½vi W þ ½b ðui vi ÞW ¼ 0 ot oxi ovi v
ð2:2Þ
Application of the grid filtering process to Eq. (2.1) gives Z mp e e Si ¼ ð ug ð2:4Þ i W vi W Þdv; qsp e for which equations are needed for ug i W and W . Applying the grid filter to (2.2), gives o e o e þ o ½bv ð~ e W ðx; v; tÞ þ ½vi W ui vi Þ W ot oxi ovi o e Þ; ½b ð ug ~i W ¼ iW u ovi v
is defined by the following functional derivative (denoted by dðÞ=dðÞ) of Xk [14]: Gjk ðx1 ; t1 ; X; tÞ ¼
and is related to the total number of particles np by Z ð2:3Þ W ðx; v; tÞdxdv ¼ np :
dXk ðtÞ : bv dfuj ðx1 ; t1 Þ huj ðx1 ; t1 Þigdt1 ð2:10Þ
The usual procedure in LES of fluid phase is to assume the model for subgrid scale stress having a form which is similar to a model for Reynolds stress. Adopting this procedure for particle phase also and using (2.6), we e , ~i W propose an analogous expression for bv ½ ug iW u written as o 0 o 0 0 e e g bv ½ ui W u~i W ¼ ð2:11Þ k þ l ci W ; oxk ki ovk ki where the tensors k0ki , l0ki , and k0i are Z t dt1 s ui ðx; tÞ; uj ðx; v; tjt1 Þ G0jk ðx1 ; t1 ; x; tÞ; k0ki ¼ b2v 0
ð2:12Þ ð2:5Þ
Z
d 0 G ðx1 ; t1 ; x; tÞ; dt1 s ui ðx; tÞ; uj ðx; v; tjt1 Þ dt jk ð2:13Þ
Z t oui ðx; tÞ dt1 s ; uj ðx; v; tjt1 Þ G0jk ðx1 ; t1 ; x; tÞ: c0i ¼ b2v oxk 0 ð2:14Þ
l0ki ¼ b2v
t
0
which poses closure problem due to the unknown term ug appearing on its right-hand side. The term iW 0 0 e Þ is of the order of ug ~i W ~i ð ug iW u i W where ui ¼ ui u represents the subgrid scale velocity fluctuations. A similar type of closure problem arising in the ensemble averaged (denoted by hi) of (2.2) due to hui W i hui ihW i was solved by Reeks [12], using Kraichnan’s Lagrangian history direct interaction (LHDI). The expression derived by Reeks is o o bv ðhui W i hui ihW iÞ ¼ kki þ lki ci hW i; oxk ovk ð2:6Þ
Here for any variables A and B f A eB e sðA; BÞ ¼ AB
ð2:15Þ
and the correlations of the form sðbi ðx; tÞ; bj ðx; v; tjt1 ÞÞ, which appear in (2.12)–(2.14), can be approximated by usual exponential function with Lagrangian integral time scale Tbi bj , written as s bi ðx; tÞ; bj ðx; v; tjt1 Þ ¼ s bi ðx; tÞ; bj ðx; tÞ eðt1 tÞ=Tbi bj :
where the tensors kki , lki , and ci are given by kki ¼ b2v
Z
ð2:16Þ The Green’s function G0jk ðx1 ; t1 ; X; tÞ is defined as
t
dt1 ½hui ðx; tÞuj ðx; v; tjt1 Þi hui ðx; tÞi 0
huj ðx; v; tjt1 ÞiGjk ðx1 ; t1 ; x; tÞ; Z t lki ¼ b2v dt1 ½hui ðx; tÞuj ðx; v; tjt1 Þi hui ðx; tÞi
ð2:7Þ
d Gjk ðx1 ; t1 ; x; tÞ; dt Z t oui ðx; tÞ dt1 uj ðx; v; tjt1 Þ ci ¼ b2v oxk 0 ohui ðx; tÞi huj ðx; v; tjt1 Þi Gjk ðx1 ; t1 ; x; tÞ: oxk
G0jk ðx1 ; t1 ; X; tÞ ¼
dXk ðtÞ bv dfuj ðx1 ; t1 Þ u~j ðx1 ; t1 Þgdt1
ð2:17Þ
and is governed by
0
huj ðx; v; tjt1 Þi
4755
ð2:8Þ
ð2:9Þ
Here uj ðx; v; tjt1 Þ is velocity of fluid in the vicinity of the particle at time t1 which had (or will have) a velocity v at position x at time t. The Green’s function Gjk ðx1 ; t1 ; X; tÞ
d2 0 d o~ uk Gjk ðx1 ; t1 ; x; tÞ þ bv G0jk bv G0ji ¼ djk dðt t1 Þ: 2 dt dt oxi ð2:18Þ W from (2.11) into (2.4) and after Substituting for ug i integration we obtain mp e mp o 0 e 0e e Si ¼ ðk N Þ ci N N ð~ ui vei Þ þ ð2:19Þ qsp q oxk ki where Z e ¼ e ðx; v; tÞdv; N W
vei ¼
1 e N
Z
e ðx; v; tÞdv; vi W
ð2:20Þ
4756
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Z
for which equations can be obtained by taking moments of the filtered Eq. (2.5) after substituting for ug i W from (2.11). These equations are e oN o e ¼ 0; þ ½ vei N ð2:21Þ ot oxi
d 00 G ðx1 ; t1 ; x; tÞ; dt1 T ui ðx; tÞ; uj ðx; v; tjt1 Þ dt jk 0 ð2:29Þ
Z t oui ðx; tÞ c00i ¼ b2v dt1 T ; uj ðx; v; tjt1 Þ G00jk ðx1 ; t1 ; x; tÞ: oxk 0
and o~ vj o~ vj 1 o e þ bv ð~ þ v~i ¼ ½ðg vi vj v~i v~j Þ N uj v~j Þ e ot oxi N oxi 0 e okkj o ln N þ c0j ; ð2:22Þ k0kj oxk oxk
ð2:30Þ
where 1 g vi vj ¼ e N
Z
e ðx; v; tÞdv vi vj W
t
Here for any variables A and B f A eB e: T ðA; BÞ ¼ AB
ð2:31Þ
Green’s function G00jk ðx1 ; t1 ; X; tÞ ¼
ð2:23Þ
dXk ðtÞ ; bv dfuj ðx1 ; t1 Þ u~j ðx1 ; t1 Þgdt1
ð2:32Þ
is governed by
and is governed by Z o e o e dv ½N g vi vj þ vi vj vn W ot oxn e ð~ e ðl0 þ l0 þ c0 v~j þ c0 v~i Þ vi u~j þ v~j u~i 2~ vi v~j Þ þ N ¼ bv N ij ji i j
l00ki ¼ b2v
o 0 e o 0 e ½k v~i N ½k v~j N : oxk kj oxk ki
ð2:24Þ
R e dv, it would If we now write an equation for R vi vj vn W e dv posing contain the higher order term vi vj vn vm W closure problem and so on. We close the set of equations (2.21)–(2.24) by considering third-order correlations of particle velocity fluctuation over v~i , to be approximately equal to zero, that is Z e ðx; v; tÞdv ’ 0 ðvi v~i Þðvj v~j Þðvn v~n Þ W ð2:25Þ and thus Z e dv ’ N e ½~ ~j vg ~n g vi vj 2~ vi vg vi v~j v~n : vi vj vn W j vn þ v i vn þ v ð2:26Þ Eqs. (2.21), (2.22), and (2.24) along with (2.26) are in the Eulerian framework and can be considered as the LES ‘fluid’ equations for the particle phase. In the framework of dynamic localization model, we e also need expression for e S i . Now, application of filter G e to (2.2) produces equation for W with unknown terms ug i W . Following a procedure similar to that described e ~i , we can above to obtain ug i W and equations for N and v write h i e ¼ o k00ki þ o l00ki c00i W e; ~i W bv ug ð2:27Þ iW u oxk ovk where the tensors k00ki , l00ki , and k00i are Z t dt1 T ui ðx; tÞ; uj ðx; v; tjt1 Þ G00jk ðx1 ; t1 ; x; tÞ; k00ki ¼ b2v 0
ð2:28Þ
uk d2 00 d o~ Gjk ðx1 ; t1 ; x; tÞ þ bv G00jk bv G00ji ¼ djk dðt t1 Þ: 2 dt dt oxi ð2:33Þ e i is The expression for S mp e mp o 00 e 00 e e ui vei Þ þ N ð~ ð2:34Þ k N ci N Si ¼ qsp q oxk ki and equations for Z Z e ¼ e ðx; v; tÞdv; vei ¼ 1 e ðx; v; tÞdv; vi W N W e N Z 1 e ðx; v; tÞdv g vi vj W vi vj ¼ e N and Z e ðx; v; tÞdv; vi vj vn W
ð2:35Þ
ð2:36Þ
have the forms similar to the Eqs. (2.21), (2.22), (2.24), and (2.26) and can be obtained by changing f, k0ij , l0ij , and c0i to f, k00ij , l00ij , and c00i , respectively, in these equations. 3. Dynamic subgrid scale model for fluid In this section, we use the framework of dynamic localization model [6] and derive the model for sðui ; uj Þ which accounts for the effects of particles on the subgrid scale turbulent motion. The model for sðui ; uj Þ is written as sðui ; uj Þ
dij sðul ; ul Þ ¼ 2CDk 1=2 e S ij ; 3
ð3:1Þ
o~ u o~ ui þ oxji and k ¼ 12 sðul ; ul Þ is governed by where e S ij ¼ 12 ½ox j ok o o 1 ok þ sðui ; ui ; uj Þ þ sðp; uj Þ m ½~ uj k ¼ ot oxj oxj 2 oxj
oui oui þ sðui ; Si Þ: S ij ms ; sðui ; uj Þ e oxj oxj
ð3:2Þ
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Here, sðf ; g; hÞ for any variables f, g, and h is defined as [15]
which is modelled as [6]
sðf ; g; hÞ ¼ g fgh f~sðg; hÞ g~sðf ; hÞ h~sðf ; gÞ f~g~h~:
T ðui ; uj Þ
dij T ðul ; ul Þ ¼ 2CDK 1=2 e S ij : 3
4757
ð3:12Þ
ð3:3Þ Also, sðui ; Si Þ is source/sink term arising due to the presence of particles and is given by ~i e sðui ; Si Þ ¼ ug Si; i Si u where mp ug i Si ¼ qsp
Z
ð3:4Þ
2 g ½ ug i W ðx; v; tÞ vi ui W ðx; v; tÞdv;
u2i ¼ ui ui
ð3:5Þ
is written using the expression for Si from (2.1). Using 2 the definition of sðf ; g; hÞ from (3.3), ug i W can be written as 2 e þ uei uei W e þ sðui ; ui ; W Þ: ug ei þ sðui ; ui Þ W i W ¼ 2sðui ; W Þ u
ð3:6Þ The unknown term sðui ; ui ; W Þ is of the order of 02 0 with W 0 ¼ W W e . The third-order correlation ug i W correlation sðui ; ui ; W Þ is small in comparison to other second-order correlations in (3.6) and is not taken into account now onwards. Thus, from (2.11), (3.5), and (3.6) i eh mp N ug sðui ; ui Þ þ u~i u~i u~i v~i i Si ¼ qsp 2mp e u~i o k0ki N e u~i c0i N q oxk mp o 0 e e c0i v~i N e : ð3:7Þ kki v~i N l0ii N q oxk
Here, for any variables f, g, and h fgh f~T ðg; hÞ g~T ðf ; hÞ T ðf ; g; hÞ ¼ g h~T ðf ; gÞ f~ g~ h~;
ð3:8Þ
~i e T ðui ; Si Þ ¼ ug S i; i Si u
ð3:15Þ
and ug i Si ¼
mp qsp
Z
2 g ug W ðx; v; tÞ v W ðx; v; tÞ dv; u i i i
u2i ¼ ui ui :
ð3:16Þ
Introducing the models for different terms,
t ¼ mT ð3:9Þ
the modelled equation for k is ok o o ok o ok þ þ ½~ uj k ¼ DDk 1=2 m ot oxj oxj oxj oxj oxj
ð3:17Þ
and
and
oui oui k 3=2 ; ¼ C
; g ¼ ms oxj oxj D
ð3:14Þ
1 oK T ðui ; ui ; uj Þ þ T ðp; uj Þ ¼ DDK 1=2 2 oxj
Introducing the models for different terms, 1 ok sðui ; ui ; uj Þ þ sðp; uj Þ ¼ DDk 1=2 2 oxj
The coefficient C is assumed to be independent of filter width while writing (3.12). A recent work by Meneveau and Lund [16] suggests that the coefficient appearing in the dynamic Smagorinsky model is scale-dependent, which could also be the case for model coefficients C, D, and C . Here we do not consider the scale-dependency and take all the coefficients independent of filter width. Similar to (3.2), equation for K ¼ T ðul ; ul Þ=2 is written as oK o o 1 oK uj K ¼ ½~ þ T ðui ; ui ; uj Þ þ T ðp; uj Þ m ot oxj oxj 2 oxj
oui oui S ij mT T ðui ; uj Þ e ; þ T ðui ; Si Þ: oxj oxj ð3:13Þ
oui oui ; oxj oxj
¼ C
K 3=2 ; D
ð3:18Þ
the modelled equation for K can be written in the form oK o o oK o oK þ uj K ¼ þ ½~ DDK 1=2 m ot oxj oxj oxj oxj oxj 3=2
3=2
k sðui ; uj Þ e þ sðui ; Si Þ: ð3:10Þ S ij C
D
K þ T ðui ; Si Þ: T ðui ; uj Þ e S ij C
D ð3:19Þ
Now we discuss the method to find the model coefficients C, D, and C appearing in (3.1), (3.8) and (3.9). e to (1.3) produces subgrid scale Application of filter G stress
3.1. Determination of coefficient C
~i u~j ; T ðui ; uj Þ ¼ ug i uj u
Lij ¼ T ðui ; uj Þ sðui ; uj Þ ¼ u~i u~j u~i u~j ;
ð3:11Þ
Eqs. (1.7) and (3.11) give the Germano identity ð3:20Þ
4758
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
where sðui ; uj Þ is obtained by applying the filter G to s. Now, from (3.1) and (3.12) Lij
dij Lll ¼ aij C Cbij ; 3
ð3:21Þ
u~i sðui ; uj Þ þ u~i ½Lij þ sðui ; uj Þ Lj :
where S ij ; aij ¼ 2DK 1=2 e
S ij : bij ¼ 2Dk 1=2 e
ð3:22Þ
Ghosal et al. [6] have proposed a variational formulation which obtains C such that it minimizes integral F½C, over the entire flow domain x, of the squares of the error Eij ðxÞ ¼ Lij that is dF ¼ 2
oK ok DDk 1=2 oxj oxj ! ! u~i u~i u~i u~i u~j p~ þ k þ ¼ u~j p~ þ k þ 2 2
DDK 1=2
Z
dij Lll aij C þ Cbij ; 3
Eij ðxÞdEij ðxÞdx ¼ 0;
ð3:23Þ
ð3:24Þ
under the conditions when Lij , aij , and bij are known. In general, Lij , aij , and bij are not known prior to the calculations of u~i and u~i , and their values can be computed by using LES equations. In view of this, any change in C would cause a change in Lij , aij , and bij as the governing equations for the velocity field u~i and u~i depend on C. This dependency makes it difficult to solve the variational problem posed by (3.24) when u~i and u~i are not known. Thus the solution given by Ghosal et al. [6] can be considered as a first approximation as it does not include the terms accounting for the dependency. Then the approximate equation obtained from (3.24) is [6] Z aij Eij dC þ Eij dCbij dx ¼ 0 ð3:25Þ and for the condition C P 0 it gives Z CðxÞ ¼ H 1 aij ðxÞLij ðxÞ bij ðxÞ Lij ðyÞGðy; xÞdy Z þ aij ðxÞ bij ðyÞCðyÞGðx; yÞdy Z þ bij ðxÞ aij ðyÞCðyÞGðy; xÞdy Z Z bij ðxÞ dy bij ðyÞCðyÞ dzGðz; yÞGðz; xÞ ; ð3:26Þ with H ¼ aij ðxÞaij ðxÞ and CðxÞ ¼ 0 when (3.26) gives negative values for C. While writing (3.26), condition of incompressibility is used. 3.2. Determination of coefficient D Applying filter G to (3.8) and then subtracting from (3.17), we obtain
To obtain D, the integral Z oK ok þ DDk 1=2 D½D ¼ Lj DDK 1=2 oxj oxj oK ok dy;
Lj DDK 1=2 þ DDk 1=2 oxj oxj
ð3:27Þ
ð3:28Þ
is to be minimized. A dD change in D would also cause changes in K, k, and Lj . Neglecting these changes, an approximate solution for D with the constraint D P 0 can be written as Z 1 DðxÞ ¼ J aj ðxÞLj ðxÞ bj ðxÞ Lj ðyÞGðy; xÞdy Z þ aj ðxÞ bj ðyÞDðyÞGðx; yÞdy Z þ bj ðxÞ aj ðyÞDðyÞGðy; xÞdy Z Z bj ðxÞ dy bj ðyÞDðyÞ dzGðz; yÞGðz; xÞ ð3:29Þ and D ¼ 0 when (3.29) gives negative value for D. Here, J ¼ al ðxÞal ðxÞ, aj ¼ DK 1=2 oK=oxj , bj ¼ Dk 1=2 ok=oxj , and the result given by (3.29) is identical to that given by Ghosal et al. [6,17] with some difference in notations. 3.3. Determination of coefficient C
To determine C following the procedure given by Ghosal et al. [6], we obtain C K 3=2 C k 3=2 ¼ v; D D where 1 oLii o 1 o oLii uj Lii Þ þ þ ð~ m v¼ 2 ot oxj 2 oxj oxj
ð3:30Þ
S ij þ sðui ; uj Þ e S ij S ij sðui ; uj Þ e Lij e 2 ! ! u~i u~i u~i u~i o 4 u~j p~ þ u~j p~ þ þ 2 2 oxj 3
h i þ u~i ðLij þ sðui ; uj ÞÞ u~i sðui ; uj Þ5 þ u~i e S i u~i e Si : ð3:31Þ
It should be noted that the last term in the square bracket on the right-hand side of (3.31) accounts for the effect of particles and appears explicitly in the
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
4759
determination of coefficient C . Whereas, the form of the expressions for C and D are identical to that obtained by Ghosal et al. [6,17]. And the effect of particles appears in coefficients C and D only through the modification of LES flow fields due to the presence of the source term Si . The variational problem of minimizing the integral #2 Z " C K 3=2 C k 3=2 þ v dy; ð3:32Þ H½C ¼ D D
Acknowledgements
gives an approximate solution for C when the changes in K, k, and v due to the change dC in C are neglected. The approximate solution with constraint C P 0 can be written in the form Z 1 aðxÞvðxÞ bðxÞ vðyÞGðy; xÞdy C ðxÞ ¼ aðxÞaðxÞ Z þ aðxÞ bðyÞC ðyÞGðx; yÞdy Z þ bðxÞ aðyÞC ðyÞGðy; xÞdy Z Z bðxÞ dy bðyÞC ðyÞ dzGðz; yÞGðz; xÞ
[1] J. Smagorinsky, General circulation experiments with the primitive equations. I. The basic experiment, Mon. Weather Rev. 91 (3) (1963) 99–164. [2] C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Ann. Rev. Fluid Mech. 32 (2000) 1–32. [3] F. Mashayek, R.V.R. Pandya, Analytical description of two-phase turbulent flows, Progress in Energy and Combustion Science, submitted for publication. [4] J. Bellan, Perspectives on large-eddy simulations for sprays: issues and solutions, Atomization Spray 10 (2000) 409–425. [5] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A 3 (7) (1991) 1760–1765. [6] S. Ghosal, T.S. Lund, P. Moin, K. Akselvoll, A dynamic localization model for large-eddy simulation of turbulent flows, J. Fluid Mech. 286 (1995) 229–255. [7] F. Yeh, U. Lei, On the motion of small particles in a homogeneous isotropic turbulent flow, Phys. Fluids 3 (11) (1991) 2571–2586. [8] F. Yeh, U. Lei, On the motion of small particles in a homogeneous turbulent shear flow, Phys. Fluids 3 (11) (1991) 2758–2776. [9] V. Armenio, U. Piomelli, V. Fiorotto, Effect of the subgrid scales on particles motion, Phys. Fluids 11 (1999) 3030– 3042. [10] M. Boivin, O. Simonin, K.D. Squires, On the prediction of gas-solid flows with two-way coupling using large eddy simulation, Phys. Fluids 12 (2000) 2080–2090. [11] S. Yuu, T. Ueno, T. Umekage, Numerical simulation of the high Reynolds number slit nozzle gas-particle jet using subgrid-scale coupling large eddy simulation, Chem. Eng. Sci. 56 (2001) 4293–4307. [12] M.W. Reeks, On the continuum equations for dispersed particles in nonuniform flows, Phys. Fluids 4 (6) (1992) 1290–1303. [13] L.I. Zaichik, A statistical model of particle transport and heat transfer in turbulent shear flows, Phys. Fluids 11 (6) (1999) 1521–1534. [14] K.E. Hyland, S. McKee, M.W. Reeks, Derivation of a PDF kinetic equation for the transport of particles in turbulent flows, J. Phys. A: Math. Gen. 32 (1999) 6169– 6190. [15] M. Germano, Turbulence: the filtering approach, J. Fluid Mech. 238 (1992) 325–336. [16] C. Meneveau, T.S. Lund, The dynamic Smagorinsky model and scale-dependent coefficients in the viscous range of turbulence, Phys. Fluids 9 (1997) 3932–3934. [17] S. Ghosal, T.S. Lund, P. Moin, K. Akselvoll, Corrigendum, J. Fluid Mech. 297 (1995) 402.
ð3:33Þ and C ¼ 0 when (3.33) gives negative value for C . Here a ¼ K 3=2 =D and b ¼ k 3=2 =D.
4. Concluding remarks The LES ‘fluid’ equations for particles have been derived which govern the large scale structure of particle phase concentration and velocity field arising due to the LES flow and subgrid scale field of fluid phase. The subgrid scales contributions in these equations appear through statistical properties of the form sðui ðx; tÞ; uj ðx; v; tjt1 ÞÞ, sðoui ðx; tÞ=oxk ; uj ðx; v; tjt1 ÞÞ, T ðui ðx; tÞ; uj ðx; v; tjt1 ÞÞ, and T ðoui ðx; tÞ=oxk ; uj ðx; v; tjt1 ÞÞ. The presence of particles influence the subgrid scale motions and their influence has been properly accounted for during the subgrid-scale stress modeling in the dynamic localization model framework. The variational method, due to Ghosal et al. [6], has been used to obtain model coefficients C, D, and C and it has been pointed out that the method gives approximate expressions due to the neglect of the dependency of various LES flow fields on these coefficients. The main objective has been to propose a new approach of LES for two-phase flows where the particle phase can also be described by the ‘fluid’ equations. Further simplifications and reduction in the number of ‘fluid’ equations by deriving constitutive relations for ðg vi vj v~i v~j Þ are possible. This remains as a part of future work on the detailed assessments of the present model equations.
This work was supported in part by the Office of Naval Research and the National Science Foundation.
References
Two-fluid large-eddy simulation approach for particle-laden turbulent flows R.V.R. Pandya, F. Mashayek
*
Department of Mechanical and Industrial Engineering, University of Illinois at Chicago, 842 West Taylor Street, Chicago, IL 60607, USA Received 1 April 2002; received in revised form 18 May 2002
Abstract In this paper, effects of particles on the subgrid scales of turbulence are properly accounted for during the modeling of subgrid scale stresses in the large-eddy simulation (LES) of fluid phase. In doing so, we propose closed filtered kinetic equations for phase space density of the particle. The various moments of these equations give the ‘fluid’ equations which can be considered as the LES equations for the particle phase. The influence of subgrid scales motion on the particles is included in these ‘fluid’ equations. Ó 2002 Elsevier Science Ltd. All rights reserved.
1. Introduction In recent years, large-eddy simulation (LES) has emerged as an important predictive technique for fluid turbulence with growing number of models for statistical properties of small scales or subgrid scales which influence the instantaneous flow properties of the large scales or eddies. The pioneering model for subgrid scale stress tensor is due to Smagorinsky [1], followed by other dynamic and similarity models which are reviewed recently by Meneveau and Katz [2]. LES is also emerging in the field of two-phase turbulent flows for accurate prediction of particle or droplet laden flows [3,4]. In these flows, each particle moves under the influence of fluid forces (e.g. fluid drag force, history force), and the instantaneous fluid flow velocity Ui in the vicinity of the particle governs its trajectory. Here, we consider only the fluid drag force and assume the particle to be spherical with time constant sp ¼ qp d 2 =18l, where l is the fluid viscosity, and qp and d are the particle mass density and diameter, respectively. In this case, the governing Lagrangian equations for the spherical particle’s position, Xi , and velocity, Vi , can be written as *
Corresponding author. Tel.: +1-312-996-5317; fax: +1-312413-0447. E-mail address: [email protected] (F. Mashayek).
dXi ¼ Vi ; dt
ð1:1Þ
dVi ¼ bv ðUi Vi Þ; dt
bv ¼
1 þ 0:15Re0:687 p ; sp
ð1:2Þ
where particle Reynolds number Rep ¼ dqjV Uj=l and q is the fluid density. In this paper, we consider small particles with Rep < 1 and take bv ¼ 1=sp . The presence of particles modifies the flow field which is accounted for through the source term Si in the governing equations for fluid phase velocity ui , written for incompressible fluid as ouj ¼ 0; oxj
oui o op o2 ui þ ½ui uj ¼ þm þ Si : oxi ot oxj oxj oxj
ð1:3Þ
Here, p represents the fluid pressure divided by q at location xi and time t, m is the fluid kinematic viscosity, and Si accounts for the effects of particles on the fluid flow and is given later in this paper. In the case of one-way coupling, Si ¼ 0 and the fluid phase is considered independent of the particle phase. Si has finite values in the case of two-way coupling. Following the work of Germano et al. [5], we denote the grid filtering operation with filter width D on any function f as Z e ðx; x0 Þdx0 ; f~ðxÞ ¼ f ðx0 Þ G ð1:4Þ
0017-9310/02/$ - see front matter Ó 2002 Elsevier Science Ltd. All rights reserved. PII: S 0 0 1 7 - 9 3 1 0 ( 0 2 ) 0 0 1 9 1 - 6
4754
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Nomenclature C C
d D ~, G G Gjk G0jk G00jk k, K mp np p Rep Si Sij S~ij ; S~ij t Tbi bj ui
subgrid stress model coefficient model coefficient in Eq. (3.9) particle diameter model coefficient in Eq. (3.8) filters Green’s function, Eq. (2.10) Green’s function, Eq. (2.17) Green’s function, Eq. (2.32) subgrid scale kinetic energies per unit mass of fluid mass of the particle total number of particles fluid pressure divided by fluid density particle Reynolds number source term in Eq. (1.3) strain rate tensor filtered strain rate tensors time Lagrangian integral time scale ith component of fluid velocity u
and the test filtering operation as Z f ðxÞ ¼ f ðx0 ÞGðx; x0 Þdx0 ;
ð1:5Þ
e filter has D as filter width. e ¼ GG and G e to the Navier–Stokes Now, application of the filter G equations gives o~ ui o o~ p o2 u~i o þ ½~ ui u~j ¼ þm sðui ; uj Þ þ e Si; oxi ot oxj oxj oxj oxj ð1:6Þ where the unknown subgrid stress tensor ~i u~j sðui ; uj Þ ¼ ug i uj u
ð1:7Þ
and e S i pose closure problems and need to be modelled. While writing (1.6), and also later in this paper, commutation between filtering and derivative operators is used. In case of a single-phase flow, various models for subgrid stress tensor are suggested and are recently reviewed by Meneveau and Katz [2]. In this paper, we use the framework of dynamic localization model, proposed by Ghosal et al. [6], and derive the model for sðui ; uj Þ in case of two-phase flow by incorporating the effects of particles on the subgrid scales. Also, a closed expression for e S i is derived from the kinetic equation for phase space density of the particle phase. In the early studies on the LES of two-phase flows with one-way coupling [7,8], particles are tracked by using the LES velocity field, in the Lagrangian Eqs. (1.1) and (1.2), instead of the instantaneous velocity field ui of
Ui Vi vi W xi Xi
ith component of fluid velocity U along the particle path ith component of particle velocity V phase space variable corresponding to Vi phase space density physical space coordinates ith component of particle position vector X
Greek symbols l fluid viscosity m fluid kinematic viscosity q fluid density qp particle mass density sp particle time constant D, D filter widths Subscript p particle Symbol h i ensemble average
the fluid phase. The effect of subgrid scale velocity field u0i ¼ ui u~i on the particle are not taken into account in these studies. Armenio et al. [9] studied these effects in detail by performing and comparing the results from independent direct numerical simulation and LES. The two-way coupling effects are incorporated in LES by Boivin et al. [10] through the source term Si . Only very recently, the two-way coupling is treated completely by also taking into account the effects of particles in the model for subgrid stress tensor [11]. Yuu et al. [11] modelled subgrid-scale turbulent mass flux of particle by gradient transport using the analogy of the molecular transport. In all of the existing studies on LES, the particle phase is simulated in Lagrangian framework (see the recent review [3] and references cited therein). Rigorous modeling of the effects of particle phase on the subgrid scales and vice versa, still, remains as an important unfinished task. In this paper, an attempt is made for such modeling and, in doing so, ‘fluid’ or Eulerian equations are derived for particles which can be considered as LES equations for the dispersed phase.
2. ‘Fluid’ equations for particle phase Using the phase space density W ðx; v; tÞ, the source term Si can be written as Z mp Si ¼ ðui vi ÞW ðx; v; tÞdv; ð2:1Þ qsp
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
where mp is mass of the particle, x and v are phase space variables corresponding to X and V, and W is governed by [3,12,13] o o o W ðx; v; tÞ þ ½vi W þ ½b ðui vi ÞW ¼ 0 ot oxi ovi v
ð2:2Þ
Application of the grid filtering process to Eq. (2.1) gives Z mp e e Si ¼ ð ug ð2:4Þ i W vi W Þdv; qsp e for which equations are needed for ug i W and W . Applying the grid filter to (2.2), gives o e o e þ o ½bv ð~ e W ðx; v; tÞ þ ½vi W ui vi Þ W ot oxi ovi o e Þ; ½b ð ug ~i W ¼ iW u ovi v
is defined by the following functional derivative (denoted by dðÞ=dðÞ) of Xk [14]: Gjk ðx1 ; t1 ; X; tÞ ¼
and is related to the total number of particles np by Z ð2:3Þ W ðx; v; tÞdxdv ¼ np :
dXk ðtÞ : bv dfuj ðx1 ; t1 Þ huj ðx1 ; t1 Þigdt1 ð2:10Þ
The usual procedure in LES of fluid phase is to assume the model for subgrid scale stress having a form which is similar to a model for Reynolds stress. Adopting this procedure for particle phase also and using (2.6), we e , ~i W propose an analogous expression for bv ½ ug iW u written as o 0 o 0 0 e e g bv ½ ui W u~i W ¼ ð2:11Þ k þ l ci W ; oxk ki ovk ki where the tensors k0ki , l0ki , and k0i are Z t dt1 s ui ðx; tÞ; uj ðx; v; tjt1 Þ G0jk ðx1 ; t1 ; x; tÞ; k0ki ¼ b2v 0
ð2:12Þ ð2:5Þ
Z
d 0 G ðx1 ; t1 ; x; tÞ; dt1 s ui ðx; tÞ; uj ðx; v; tjt1 Þ dt jk ð2:13Þ
Z t oui ðx; tÞ dt1 s ; uj ðx; v; tjt1 Þ G0jk ðx1 ; t1 ; x; tÞ: c0i ¼ b2v oxk 0 ð2:14Þ
l0ki ¼ b2v
t
0
which poses closure problem due to the unknown term ug appearing on its right-hand side. The term iW 0 0 e Þ is of the order of ug ~i W ~i ð ug iW u i W where ui ¼ ui u represents the subgrid scale velocity fluctuations. A similar type of closure problem arising in the ensemble averaged (denoted by hi) of (2.2) due to hui W i hui ihW i was solved by Reeks [12], using Kraichnan’s Lagrangian history direct interaction (LHDI). The expression derived by Reeks is o o bv ðhui W i hui ihW iÞ ¼ kki þ lki ci hW i; oxk ovk ð2:6Þ
Here for any variables A and B f A eB e sðA; BÞ ¼ AB
ð2:15Þ
and the correlations of the form sðbi ðx; tÞ; bj ðx; v; tjt1 ÞÞ, which appear in (2.12)–(2.14), can be approximated by usual exponential function with Lagrangian integral time scale Tbi bj , written as s bi ðx; tÞ; bj ðx; v; tjt1 Þ ¼ s bi ðx; tÞ; bj ðx; tÞ eðt1 tÞ=Tbi bj :
where the tensors kki , lki , and ci are given by kki ¼ b2v
Z
ð2:16Þ The Green’s function G0jk ðx1 ; t1 ; X; tÞ is defined as
t
dt1 ½hui ðx; tÞuj ðx; v; tjt1 Þi hui ðx; tÞi 0
huj ðx; v; tjt1 ÞiGjk ðx1 ; t1 ; x; tÞ; Z t lki ¼ b2v dt1 ½hui ðx; tÞuj ðx; v; tjt1 Þi hui ðx; tÞi
ð2:7Þ
d Gjk ðx1 ; t1 ; x; tÞ; dt Z t oui ðx; tÞ dt1 uj ðx; v; tjt1 Þ ci ¼ b2v oxk 0 ohui ðx; tÞi huj ðx; v; tjt1 Þi Gjk ðx1 ; t1 ; x; tÞ: oxk
G0jk ðx1 ; t1 ; X; tÞ ¼
dXk ðtÞ bv dfuj ðx1 ; t1 Þ u~j ðx1 ; t1 Þgdt1
ð2:17Þ
and is governed by
0
huj ðx; v; tjt1 Þi
4755
ð2:8Þ
ð2:9Þ
Here uj ðx; v; tjt1 Þ is velocity of fluid in the vicinity of the particle at time t1 which had (or will have) a velocity v at position x at time t. The Green’s function Gjk ðx1 ; t1 ; X; tÞ
d2 0 d o~ uk Gjk ðx1 ; t1 ; x; tÞ þ bv G0jk bv G0ji ¼ djk dðt t1 Þ: 2 dt dt oxi ð2:18Þ W from (2.11) into (2.4) and after Substituting for ug i integration we obtain mp e mp o 0 e 0e e Si ¼ ðk N Þ ci N N ð~ ui vei Þ þ ð2:19Þ qsp q oxk ki where Z e ¼ e ðx; v; tÞdv; N W
vei ¼
1 e N
Z
e ðx; v; tÞdv; vi W
ð2:20Þ
4756
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Z
for which equations can be obtained by taking moments of the filtered Eq. (2.5) after substituting for ug i W from (2.11). These equations are e oN o e ¼ 0; þ ½ vei N ð2:21Þ ot oxi
d 00 G ðx1 ; t1 ; x; tÞ; dt1 T ui ðx; tÞ; uj ðx; v; tjt1 Þ dt jk 0 ð2:29Þ
Z t oui ðx; tÞ c00i ¼ b2v dt1 T ; uj ðx; v; tjt1 Þ G00jk ðx1 ; t1 ; x; tÞ: oxk 0
and o~ vj o~ vj 1 o e þ bv ð~ þ v~i ¼ ½ðg vi vj v~i v~j Þ N uj v~j Þ e ot oxi N oxi 0 e okkj o ln N þ c0j ; ð2:22Þ k0kj oxk oxk
ð2:30Þ
where 1 g vi vj ¼ e N
Z
e ðx; v; tÞdv vi vj W
t
Here for any variables A and B f A eB e: T ðA; BÞ ¼ AB
ð2:31Þ
Green’s function G00jk ðx1 ; t1 ; X; tÞ ¼
ð2:23Þ
dXk ðtÞ ; bv dfuj ðx1 ; t1 Þ u~j ðx1 ; t1 Þgdt1
ð2:32Þ
is governed by
and is governed by Z o e o e dv ½N g vi vj þ vi vj vn W ot oxn e ð~ e ðl0 þ l0 þ c0 v~j þ c0 v~i Þ vi u~j þ v~j u~i 2~ vi v~j Þ þ N ¼ bv N ij ji i j
l00ki ¼ b2v
o 0 e o 0 e ½k v~i N ½k v~j N : oxk kj oxk ki
ð2:24Þ
R e dv, it would If we now write an equation for R vi vj vn W e dv posing contain the higher order term vi vj vn vm W closure problem and so on. We close the set of equations (2.21)–(2.24) by considering third-order correlations of particle velocity fluctuation over v~i , to be approximately equal to zero, that is Z e ðx; v; tÞdv ’ 0 ðvi v~i Þðvj v~j Þðvn v~n Þ W ð2:25Þ and thus Z e dv ’ N e ½~ ~j vg ~n g vi vj 2~ vi vg vi v~j v~n : vi vj vn W j vn þ v i vn þ v ð2:26Þ Eqs. (2.21), (2.22), and (2.24) along with (2.26) are in the Eulerian framework and can be considered as the LES ‘fluid’ equations for the particle phase. In the framework of dynamic localization model, we e also need expression for e S i . Now, application of filter G e to (2.2) produces equation for W with unknown terms ug i W . Following a procedure similar to that described e ~i , we can above to obtain ug i W and equations for N and v write h i e ¼ o k00ki þ o l00ki c00i W e; ~i W bv ug ð2:27Þ iW u oxk ovk where the tensors k00ki , l00ki , and k00i are Z t dt1 T ui ðx; tÞ; uj ðx; v; tjt1 Þ G00jk ðx1 ; t1 ; x; tÞ; k00ki ¼ b2v 0
ð2:28Þ
uk d2 00 d o~ Gjk ðx1 ; t1 ; x; tÞ þ bv G00jk bv G00ji ¼ djk dðt t1 Þ: 2 dt dt oxi ð2:33Þ e i is The expression for S mp e mp o 00 e 00 e e ui vei Þ þ N ð~ ð2:34Þ k N ci N Si ¼ qsp q oxk ki and equations for Z Z e ¼ e ðx; v; tÞdv; vei ¼ 1 e ðx; v; tÞdv; vi W N W e N Z 1 e ðx; v; tÞdv g vi vj W vi vj ¼ e N and Z e ðx; v; tÞdv; vi vj vn W
ð2:35Þ
ð2:36Þ
have the forms similar to the Eqs. (2.21), (2.22), (2.24), and (2.26) and can be obtained by changing f, k0ij , l0ij , and c0i to f, k00ij , l00ij , and c00i , respectively, in these equations. 3. Dynamic subgrid scale model for fluid In this section, we use the framework of dynamic localization model [6] and derive the model for sðui ; uj Þ which accounts for the effects of particles on the subgrid scale turbulent motion. The model for sðui ; uj Þ is written as sðui ; uj Þ
dij sðul ; ul Þ ¼ 2CDk 1=2 e S ij ; 3
ð3:1Þ
o~ u o~ ui þ oxji and k ¼ 12 sðul ; ul Þ is governed by where e S ij ¼ 12 ½ox j ok o o 1 ok þ sðui ; ui ; uj Þ þ sðp; uj Þ m ½~ uj k ¼ ot oxj oxj 2 oxj
oui oui þ sðui ; Si Þ: S ij ms ; sðui ; uj Þ e oxj oxj
ð3:2Þ
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
Here, sðf ; g; hÞ for any variables f, g, and h is defined as [15]
which is modelled as [6]
sðf ; g; hÞ ¼ g fgh f~sðg; hÞ g~sðf ; hÞ h~sðf ; gÞ f~g~h~:
T ðui ; uj Þ
dij T ðul ; ul Þ ¼ 2CDK 1=2 e S ij : 3
4757
ð3:12Þ
ð3:3Þ Also, sðui ; Si Þ is source/sink term arising due to the presence of particles and is given by ~i e sðui ; Si Þ ¼ ug Si; i Si u where mp ug i Si ¼ qsp
Z
ð3:4Þ
2 g ½ ug i W ðx; v; tÞ vi ui W ðx; v; tÞdv;
u2i ¼ ui ui
ð3:5Þ
is written using the expression for Si from (2.1). Using 2 the definition of sðf ; g; hÞ from (3.3), ug i W can be written as 2 e þ uei uei W e þ sðui ; ui ; W Þ: ug ei þ sðui ; ui Þ W i W ¼ 2sðui ; W Þ u
ð3:6Þ The unknown term sðui ; ui ; W Þ is of the order of 02 0 with W 0 ¼ W W e . The third-order correlation ug i W correlation sðui ; ui ; W Þ is small in comparison to other second-order correlations in (3.6) and is not taken into account now onwards. Thus, from (2.11), (3.5), and (3.6) i eh mp N ug sðui ; ui Þ þ u~i u~i u~i v~i i Si ¼ qsp 2mp e u~i o k0ki N e u~i c0i N q oxk mp o 0 e e c0i v~i N e : ð3:7Þ kki v~i N l0ii N q oxk
Here, for any variables f, g, and h fgh f~T ðg; hÞ g~T ðf ; hÞ T ðf ; g; hÞ ¼ g h~T ðf ; gÞ f~ g~ h~;
ð3:8Þ
~i e T ðui ; Si Þ ¼ ug S i; i Si u
ð3:15Þ
and ug i Si ¼
mp qsp
Z
2 g ug W ðx; v; tÞ v W ðx; v; tÞ dv; u i i i
u2i ¼ ui ui :
ð3:16Þ
Introducing the models for different terms,
t ¼ mT ð3:9Þ
the modelled equation for k is ok o o ok o ok þ þ ½~ uj k ¼ DDk 1=2 m ot oxj oxj oxj oxj oxj
ð3:17Þ
and
and
oui oui k 3=2 ; ¼ C
; g ¼ ms oxj oxj D
ð3:14Þ
1 oK T ðui ; ui ; uj Þ þ T ðp; uj Þ ¼ DDK 1=2 2 oxj
Introducing the models for different terms, 1 ok sðui ; ui ; uj Þ þ sðp; uj Þ ¼ DDk 1=2 2 oxj
The coefficient C is assumed to be independent of filter width while writing (3.12). A recent work by Meneveau and Lund [16] suggests that the coefficient appearing in the dynamic Smagorinsky model is scale-dependent, which could also be the case for model coefficients C, D, and C . Here we do not consider the scale-dependency and take all the coefficients independent of filter width. Similar to (3.2), equation for K ¼ T ðul ; ul Þ=2 is written as oK o o 1 oK uj K ¼ ½~ þ T ðui ; ui ; uj Þ þ T ðp; uj Þ m ot oxj oxj 2 oxj
oui oui S ij mT T ðui ; uj Þ e ; þ T ðui ; Si Þ: oxj oxj ð3:13Þ
oui oui ; oxj oxj
¼ C
K 3=2 ; D
ð3:18Þ
the modelled equation for K can be written in the form oK o o oK o oK þ uj K ¼ þ ½~ DDK 1=2 m ot oxj oxj oxj oxj oxj 3=2
3=2
k sðui ; uj Þ e þ sðui ; Si Þ: ð3:10Þ S ij C
D
K þ T ðui ; Si Þ: T ðui ; uj Þ e S ij C
D ð3:19Þ
Now we discuss the method to find the model coefficients C, D, and C appearing in (3.1), (3.8) and (3.9). e to (1.3) produces subgrid scale Application of filter G stress
3.1. Determination of coefficient C
~i u~j ; T ðui ; uj Þ ¼ ug i uj u
Lij ¼ T ðui ; uj Þ sðui ; uj Þ ¼ u~i u~j u~i u~j ;
ð3:11Þ
Eqs. (1.7) and (3.11) give the Germano identity ð3:20Þ
4758
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
where sðui ; uj Þ is obtained by applying the filter G to s. Now, from (3.1) and (3.12) Lij
dij Lll ¼ aij C Cbij ; 3
ð3:21Þ
u~i sðui ; uj Þ þ u~i ½Lij þ sðui ; uj Þ Lj :
where S ij ; aij ¼ 2DK 1=2 e
S ij : bij ¼ 2Dk 1=2 e
ð3:22Þ
Ghosal et al. [6] have proposed a variational formulation which obtains C such that it minimizes integral F½C, over the entire flow domain x, of the squares of the error Eij ðxÞ ¼ Lij that is dF ¼ 2
oK ok DDk 1=2 oxj oxj ! ! u~i u~i u~i u~i u~j p~ þ k þ ¼ u~j p~ þ k þ 2 2
DDK 1=2
Z
dij Lll aij C þ Cbij ; 3
Eij ðxÞdEij ðxÞdx ¼ 0;
ð3:23Þ
ð3:24Þ
under the conditions when Lij , aij , and bij are known. In general, Lij , aij , and bij are not known prior to the calculations of u~i and u~i , and their values can be computed by using LES equations. In view of this, any change in C would cause a change in Lij , aij , and bij as the governing equations for the velocity field u~i and u~i depend on C. This dependency makes it difficult to solve the variational problem posed by (3.24) when u~i and u~i are not known. Thus the solution given by Ghosal et al. [6] can be considered as a first approximation as it does not include the terms accounting for the dependency. Then the approximate equation obtained from (3.24) is [6] Z aij Eij dC þ Eij dCbij dx ¼ 0 ð3:25Þ and for the condition C P 0 it gives Z CðxÞ ¼ H 1 aij ðxÞLij ðxÞ bij ðxÞ Lij ðyÞGðy; xÞdy Z þ aij ðxÞ bij ðyÞCðyÞGðx; yÞdy Z þ bij ðxÞ aij ðyÞCðyÞGðy; xÞdy Z Z bij ðxÞ dy bij ðyÞCðyÞ dzGðz; yÞGðz; xÞ ; ð3:26Þ with H ¼ aij ðxÞaij ðxÞ and CðxÞ ¼ 0 when (3.26) gives negative values for C. While writing (3.26), condition of incompressibility is used. 3.2. Determination of coefficient D Applying filter G to (3.8) and then subtracting from (3.17), we obtain
To obtain D, the integral Z oK ok þ DDk 1=2 D½D ¼ Lj DDK 1=2 oxj oxj oK ok dy;
Lj DDK 1=2 þ DDk 1=2 oxj oxj
ð3:27Þ
ð3:28Þ
is to be minimized. A dD change in D would also cause changes in K, k, and Lj . Neglecting these changes, an approximate solution for D with the constraint D P 0 can be written as Z 1 DðxÞ ¼ J aj ðxÞLj ðxÞ bj ðxÞ Lj ðyÞGðy; xÞdy Z þ aj ðxÞ bj ðyÞDðyÞGðx; yÞdy Z þ bj ðxÞ aj ðyÞDðyÞGðy; xÞdy Z Z bj ðxÞ dy bj ðyÞDðyÞ dzGðz; yÞGðz; xÞ ð3:29Þ and D ¼ 0 when (3.29) gives negative value for D. Here, J ¼ al ðxÞal ðxÞ, aj ¼ DK 1=2 oK=oxj , bj ¼ Dk 1=2 ok=oxj , and the result given by (3.29) is identical to that given by Ghosal et al. [6,17] with some difference in notations. 3.3. Determination of coefficient C
To determine C following the procedure given by Ghosal et al. [6], we obtain C K 3=2 C k 3=2 ¼ v; D D where 1 oLii o 1 o oLii uj Lii Þ þ þ ð~ m v¼ 2 ot oxj 2 oxj oxj
ð3:30Þ
S ij þ sðui ; uj Þ e S ij S ij sðui ; uj Þ e Lij e 2 ! ! u~i u~i u~i u~i o 4 u~j p~ þ u~j p~ þ þ 2 2 oxj 3
h i þ u~i ðLij þ sðui ; uj ÞÞ u~i sðui ; uj Þ5 þ u~i e S i u~i e Si : ð3:31Þ
It should be noted that the last term in the square bracket on the right-hand side of (3.31) accounts for the effect of particles and appears explicitly in the
R.V.R. Pandya, F. Mashayek / International Journal of Heat and Mass Transfer 45 (2002) 4753–4759
4759
determination of coefficient C . Whereas, the form of the expressions for C and D are identical to that obtained by Ghosal et al. [6,17]. And the effect of particles appears in coefficients C and D only through the modification of LES flow fields due to the presence of the source term Si . The variational problem of minimizing the integral #2 Z " C K 3=2 C k 3=2 þ v dy; ð3:32Þ H½C ¼ D D
Acknowledgements
gives an approximate solution for C when the changes in K, k, and v due to the change dC in C are neglected. The approximate solution with constraint C P 0 can be written in the form Z 1 aðxÞvðxÞ bðxÞ vðyÞGðy; xÞdy C ðxÞ ¼ aðxÞaðxÞ Z þ aðxÞ bðyÞC ðyÞGðx; yÞdy Z þ bðxÞ aðyÞC ðyÞGðy; xÞdy Z Z bðxÞ dy bðyÞC ðyÞ dzGðz; yÞGðz; xÞ
[1] J. Smagorinsky, General circulation experiments with the primitive equations. I. The basic experiment, Mon. Weather Rev. 91 (3) (1963) 99–164. [2] C. Meneveau, J. Katz, Scale-invariance and turbulence models for large-eddy simulation, Ann. Rev. Fluid Mech. 32 (2000) 1–32. [3] F. Mashayek, R.V.R. Pandya, Analytical description of two-phase turbulent flows, Progress in Energy and Combustion Science, submitted for publication. [4] J. Bellan, Perspectives on large-eddy simulations for sprays: issues and solutions, Atomization Spray 10 (2000) 409–425. [5] M. Germano, U. Piomelli, P. Moin, W.H. Cabot, A dynamic subgrid-scale eddy viscosity model, Phys. Fluids A 3 (7) (1991) 1760–1765. [6] S. Ghosal, T.S. Lund, P. Moin, K. Akselvoll, A dynamic localization model for large-eddy simulation of turbulent flows, J. Fluid Mech. 286 (1995) 229–255. [7] F. Yeh, U. Lei, On the motion of small particles in a homogeneous isotropic turbulent flow, Phys. Fluids 3 (11) (1991) 2571–2586. [8] F. Yeh, U. Lei, On the motion of small particles in a homogeneous turbulent shear flow, Phys. Fluids 3 (11) (1991) 2758–2776. [9] V. Armenio, U. Piomelli, V. Fiorotto, Effect of the subgrid scales on particles motion, Phys. Fluids 11 (1999) 3030– 3042. [10] M. Boivin, O. Simonin, K.D. Squires, On the prediction of gas-solid flows with two-way coupling using large eddy simulation, Phys. Fluids 12 (2000) 2080–2090. [11] S. Yuu, T. Ueno, T. Umekage, Numerical simulation of the high Reynolds number slit nozzle gas-particle jet using subgrid-scale coupling large eddy simulation, Chem. Eng. Sci. 56 (2001) 4293–4307. [12] M.W. Reeks, On the continuum equations for dispersed particles in nonuniform flows, Phys. Fluids 4 (6) (1992) 1290–1303. [13] L.I. Zaichik, A statistical model of particle transport and heat transfer in turbulent shear flows, Phys. Fluids 11 (6) (1999) 1521–1534. [14] K.E. Hyland, S. McKee, M.W. Reeks, Derivation of a PDF kinetic equation for the transport of particles in turbulent flows, J. Phys. A: Math. Gen. 32 (1999) 6169– 6190. [15] M. Germano, Turbulence: the filtering approach, J. Fluid Mech. 238 (1992) 325–336. [16] C. Meneveau, T.S. Lund, The dynamic Smagorinsky model and scale-dependent coefficients in the viscous range of turbulence, Phys. Fluids 9 (1997) 3932–3934. [17] S. Ghosal, T.S. Lund, P. Moin, K. Akselvoll, Corrigendum, J. Fluid Mech. 297 (1995) 402.
ð3:33Þ and C ¼ 0 when (3.33) gives negative value for C . Here a ¼ K 3=2 =D and b ¼ k 3=2 =D.
4. Concluding remarks The LES ‘fluid’ equations for particles have been derived which govern the large scale structure of particle phase concentration and velocity field arising due to the LES flow and subgrid scale field of fluid phase. The subgrid scales contributions in these equations appear through statistical properties of the form sðui ðx; tÞ; uj ðx; v; tjt1 ÞÞ, sðoui ðx; tÞ=oxk ; uj ðx; v; tjt1 ÞÞ, T ðui ðx; tÞ; uj ðx; v; tjt1 ÞÞ, and T ðoui ðx; tÞ=oxk ; uj ðx; v; tjt1 ÞÞ. The presence of particles influence the subgrid scale motions and their influence has been properly accounted for during the subgrid-scale stress modeling in the dynamic localization model framework. The variational method, due to Ghosal et al. [6], has been used to obtain model coefficients C, D, and C and it has been pointed out that the method gives approximate expressions due to the neglect of the dependency of various LES flow fields on these coefficients. The main objective has been to propose a new approach of LES for two-phase flows where the particle phase can also be described by the ‘fluid’ equations. Further simplifications and reduction in the number of ‘fluid’ equations by deriving constitutive relations for ðg vi vj v~i v~j Þ are possible. This remains as a part of future work on the detailed assessments of the present model equations.
This work was supported in part by the Office of Naval Research and the National Science Foundation.
References