Flow Turbulence Combust (2007) 79:133–154 DOI 10.1007/s10494-007-9076-4

Analysis and Modeling of the Turbulent Diffusion of Turbulent Kinetic Energy in Natural Convection Laltu Chandra & Günther Grötzbach

Received: 28 August 2006 / Accepted: 2 March 2007 / Published online: 4 April 2007 # Springer Science + Business Media B.V. 2007

Abstract Buoyant flows often contain regions with unstable and stable thermal stratification from which counter gradient turbulent fluxes are resulting, e.g. fluxes of heat or of any turbulence quantity. Basing on investigations in meteorology an improvement in the standard gradient-diffusion model for turbulent diffusion of turbulent kinetic energy is discussed. The two closure terms of the turbulent diffusion, the velocity-fluctuation triple correlation and the velocity-pressure fluctuation correlation, are investigated based on Direct Numerical Simulation (DNS) data for an internally heated fluid layer and for Rayleigh–Bénard convection. As a result it is decided to extend the standard gradient-diffusion model for the turbulent energy diffusion by modeling its closure terms separately. Coupling of two models leads to an extended RANS model for the turbulent energy diffusion. The involved closure term, the turbulent diffusion of heat flux, is studied based on its transport equation. This results in a buoyancy-extended version of the Daly and Harlow model. The models for all closure terms and for the turbulent energy diffusion are validated with the help of DNS data for internally heated fluid layers with Prandtl number Pr=7 and for Rayleigh–Bénard convection with Pr=0.71. It is found that the buoyancy-extended diffusion model which involves also a transport equation for the variance of the vertical velocity fluctuation gives improved turbulent energy diffusion data for the combined case with local stable and unstable stratification and that it allows for the required counter gradient energy flux. Keywords Buoyant flows . Turbulent energy diffusion . Counter gradient flux

G. Grötzbach (*) Forschungszentrum Karlsruhe GmbH, Institut für Kern- und Energietechnik, Postfach 3640, 76021 Karlsruhe, Germany e-mail: [email protected] L. Chandra Equipe MoST/LEGI, B.P. 53, 38041 Grenoble, Cedex 09, France e-mail: [email protected]

134

Flow Turbulence Combust (2007) 79:133–154

1 Introduction Buoyant convection occurs in many technical applications. The basics of buoyant convection have been studied in meteorological and industrial flows, see e.g. Lumley et al. [34], Moeng and Wyngaard [36], and Canuto and Christensen-Dalsgaard [1]. In such flows the thermal stratification plays an important role: e.g. in Rayleigh–Bénard Convection (RBC) the fluid is unstably stratified which means, colder heavier fluid is above warmer lighter fluid and thus the vertical exchange of mass, momentum and heat is amplified. Whereas in an Internally Heated fluid Layer (IHL) and in the Convective Boundary Layer (CBL) unstable situations may occur together with stable stratification; the latter means, warmer lighter fluid is over colder heavier fluid which attenuates the vertical exchange and damps turbulence. This damping cannot be accounted for by isotropic first order ReynoldsAveraged Navier Stokes (RANS) models, like the standard E 0  e0 model, see e.g. Davidson [7]. Here the notation for turbulent kinetic energy E 0 and its dissipation e0 are used instead of k and e. The RANS method is preferred in numerical studies of engineering flow problems. Some tasks involve thermally stratified convection in which the Rayleigh number Ra can attain high values. The quality of the results strongly depends on the reliability of the turbulence models. Dinh and Nourgaliev [10] have demonstrated that the standard low Reynolds number E 0  e0  σt model (where σt is the turbulent Prandtl number for heat) is inadequate to recalculate an experiment of natural convection with an internally heated fluid and cooling from above. The model fails to describe the mean temperature and heat transfer in the regime of interest. One reason for this failure can be found in the inadequacy of the turbulent Prandtl number and gradient-diffusion assumption for the turbulent heat flux, because the IHL shows a counter gradient heat flux over about 50% of the height of the fluid layer [18, 21]. According to Schumann [45] such problems occur in all flows with stable stratification when the sink terms of thermal variances are smaller than the corresponding production terms. In such areas the counter gradient heat flux converts potential energy into kinetic energy. Another reason for the failure may be found in the turbulent momentum transfer model: Standard as well as advanced RANS models, see e.g. Hanjalić et al. [25], Carteciano [2], Carteciano and Grötzbach [3], Otić et al. [39], Hanjalić [26], use the transport equation for the turbulent kinetic energy E 0. In this equation the turbulent diffusion is one of the closure 0 terms. It consists of the turbulent-transport (velocity-fluctuation triple correlation uj E 0 ) and 0 0 the pressure-transport (velocity-pressure fluctuation correlation uj p ). Both are modeled together by means of a gradient-diffusion model, Launder and Spalding [29]. Wörner et al. [49] used Direct Numerical Simulation (DNS) data and showed that this standard model is not adequate in IHL. Areas with counter gradient energy fluxes occur which cannot be reproduced by a gradient-diffusion model. The need for modeling the turbulent diffusion occurs analogously in some Large Eddy Simulations (LES) in the transport equation for the kinetic energy in the sub-grid scales. E.g., such a formally adapted gradient-diffusion model was employed by Deardorff [8] in an LES of the planetary boundary layer. When stable stratification is engaged, the models apply modified length scales to enforce the required damping of turbulence like in Deardorff [9] and Schmidt and Schumann [43]. Applying this approach of modifying locally the length scales in LES to IHL allows according to Seiter [46] for some improvement in LES, but the results are not so convincing that one could expect that a similar approach applied to RANS models would really solve the problem in IHL.

Flow Turbulence Combust (2007) 79:133–154

135

The inadequacy of the gradient model for the turbulent energy diffusion was investigated by Moeng and Wyngaard [36] by using a LES study of the CBL. They observed that the gradient model can indeed not explain the counter gradient transport of E 0 . The importance of the counter gradient energy transport in the upper part of the atmospheric boundary layer was described by Lumley et al. [34]. According to the authors, the anisotropic effect of buoyancy could explain the counter gradient energy flux. Therefore, Moeng and Wyngaard [36] recommended to include the buoyancy effect in the modeling of the turbulent transport of E 0 . 0 0 From analyses of DNS data Wörner and Grötzbach [50] deduced that uj E 0 and uj p0 are of different importance in IHL and RBC: In RBC the pressure correlation is the dominant one, whereas in IHL the triple correlation dominates the turbulent diffusion. This is due to differences in the governing intermittent transport mechanisms. Based on a LES of airflow above forest canopy, Dwyer et al. [16] found that the pressure-transport in the turbulent diffusion plays a 0 0 significant role. These studies indicate that uj E 0 and uj p0 should not be modeled together in IHL and RBC. Moreover, modeling the pressure term needs special attention. In this paper we concentrate on an improved modeling of the turbulent diffusion of kinetic energy for buoyant flows, especially for flows with local stable stratification which may cause counter gradient energy fluxes. The problem and some features of the two different thermally stratified flow types, namely IHL and RBC, will be described based on 0 0 data analyzed from DNS results. Considering the above studies the terms uj E 0 and uj p0 will 0 0 be modeled separately. For the model deduction the transport equation for uj E and the 0 Donaldson [14] approximation for uj p0 will be used as the starting point. The transport 0 equation for one of the closure terms in the model for uj E 0, namely the turbulent convection 0 2 0 of the turbulent heat flux u3 T , will be analyzed. The study will be resulting in a buoyancyextended model for this closure term. Incorporating the model for turbulent convection of turbulent heat fluxes will lead to a buoyancy-extended model for the turbulent diffusion of E 0 . The model and all its contributions will be validated by using DNS data.

2 Modeling Requirements Based on DNS Data 2.1 DNS specifications and results This section deals with characteristic data of IHL and RBC analyzed from DNS results. In RBC, the fluid layer between two infinite horizontal isothermal walls is heated uniformly from below and cooled from above. The external Rayleigh number is RaE =ggΔTwD3/(3k), where g is the gravity acceleration, g is the volume expansion coefficient, ΔTw is the wall temperature difference, D is the wall distance, ν and k are the diffusivities for momentum, respectively thermal energy. The Prandtl number of the fluid is Pr=ν/κ. In IHL, the fluid is heated by a uniform volumetric energy source qv and cooling is by keeping both walls at a lower temperature than the fluid confined in-between. The internal Rayleigh number is RaI = ggqvD5/(νκl), where l is the thermal conductivity. Hereafter, internal and external Rayleigh numbers are referred to as Rayleigh number Ra. The DNS of RBC and IHL at different Rayleigh numbers and Prandtl numbers are performed with the three-dimensional time-dependent TURBIT code, which is based on a finite volume method, see Schumann [44], Grötzbach [21], Wörner [48]. Results of this code have been intensively used and validated for RBC in different fluids, e.g. in Grötzbach [19, 20, 23], Wörner [48], Wörner and Grötzbach [50], and in Otić et al. [39]. Analyses and validation for the IHL are found e.g. in Grötzbach [17, 18, 21, 22], Wörner et al. [49],

136

Flow Turbulence Combust (2007) 79:133–154

Table 1 DNS datasets used in the analysis Flow type RBC IHL IHL

Rayleigh number Ra 5

6.3×10 5×106, 107, 108 109

Prandtl number Pr

Source of DNS data

0.71 7 7

Wörner [48] Wörner et al. [49] Chandra [4]

Wörner and Grötzbach [50], and Chandra [4]. A review on the knowledge of the convective heat transfer from the less recognized IHL is given by Kulacki and Richards [28]. Table 1 gives the parameters of the simulations which form the basis for the present analysis and model development. The spatial resolution requirements by Grötzbach [17, 20] are fulfilled by all these DNS. This means, all grids are chosen fine enough to resolve everywhere the small scales in the turbulence fields, and the channel size is chosen large enough so that the periodic boundary conditions in the horizontal directions do not hinder the large scale flow structures. The validation of these simulations has been performed in the referred publications. The statistical mean values are calculated by assuming that both flow types are homogeneous in the horizontal plane x1−x2 and thus averaging is performed over these planes and over time; such averages are denoted by an over-bar ðÞ. In all figures this surface-time average is represented by 〈 〉 and the fluctuation with respect to this will be denoted by ( )″. As a result of this averaging the heat transfer in each fluid layer reduces to a one-dimensional problem depending only on the vertical coordinate x3 and as there is no horizontal mean flow, there exists no mean shear. For scaling the wall distance D is used, as well as the wall temperature difference ΔTw, the velocity scale u0 =(ggΔTwD)1/2, and the
2 , ffiffiffiffiffiffiffiffiffiffiffiffi where ρ is the density. The present scaling results in pressure scalepffiffiffiffiffi rffi u0 p Re ¼ u0 D=ν ¼ Gr ¼ Ra=Pr. The such analyzed scaled mean temperature profiles in IHL and RBC at different Ra and Pr in Fig. 1 show, the IHL of water is stably thermally stratified with a very thin unstably stratified upper thermal boundary layer. Here x3 =1 and x3 =2 indicate the vertical positions of the lower and upper walls. In IHL the maximum value of the mean temperature occurs near the edge of the upper boundary layer. This thin upper unstable layer drives the turbulent heat and momentum exchange. At the other vertical positions the fluid is stably stratified which attenuates the turbulent exchange. The thicknesses of the two thermal boundary layers indicate that most of the homogeneously released heat in the fluid is transferred upward. Indeed, even a short estimation from the Nusselt numbers at both walls shows, the heat flux is positive for the upper 70% of the fluid layer, which means the heat is transferred in the upward direction. Thus, a strong counter gradient heat flux occurs at all Ra between roughly x3 =1.3 and the position of the temperature maximum at about x3 =1.9. Detailed profiles are given in Grötzbach [21, 22]. The RBC with air is unstably thermally stratified throughout the height of the channel. The thickness of the thermal boundary layers is at both walls the same. Apart from walls the temperature field shows a widely isothermal core. 1=2 The Root Mean Square (RMS) values of the velocity fluctuations urms ¼ u0 2 ¼ 1=2 ðu  uÞ2 evaluated from the DNS data are shown in Fig. 2. Here ( )′ indicates the fluctuation, which is calculated using ðÞ. The strong anisotropy between the velocity fluctuations is attributed to both the effect of buoyancy and the presence of walls in both flow types. In IHL most of the turbulent kinetic energy is apart from the walls contained in the vertical velocity component. Both the presence of the maximum and the extrusions in the horizontal velocity components indicate a considerable energy transfer from the vertical

Flow Turbulence Combust (2007) 79:133–154

137

1.0

IHL, Pr=7.0

RBC

1.0

0.8 Ra=5E6 Ra=1E7 Ra=1E8 Ra=1E9

0.6 0.4





0.8

0.4 0.2

0.2 0.0 1.0

0.6

1.2

1.4

X3

1.6

1.8

2.0

0.0 1.0

Ra =6.3E5, Pr=0.71

1.2

1.4

1.6 X3

1.8

2.0

Fig. 1 Vertical profiles of time mean temperature T in IHL and RBC

to the horizontal velocity components close to the walls. This is consistent with large values of the pressure strain terms in those areas, Grötzbach [22]. Similarly, in RBC the appearance of maxima in the RMS values of each of the horizontal velocity components close to the walls is caused by energy transfer from the vertical to the horizontal components. The differences between the RMS values of both horizontal components can be explained by the presence of large-scale roll structures which are superimposed to the more random plume structures as could be concluded from these simulation results, see Grötzbach [23]. Recently, the occurrence of such roll systems in RBC at similar Ra numbers was systematically analyzed and classified by Hartlep et al. [27] using DNS from large computational domains and spectral analysis. 0 2 The vertical profiles of the Turbulent Kinetic Energy (TKE) E 0 ¼ 1=2ðui Þ in Fig. 3 0 show an increase of the maximum value of E with Ra in IHL. The maxima for the two highest Ra are about the same. This indicates that the flow is tending with increasing Ra towards the fully turbulent regime. Further on, the vertical profiles of E 0 depict that at these Ra the convection in the IHL is still in-homogeneous in the vertical direction in the core of the fluid layer, whereas in RBC a wide area is observed in which the energy distribution is homogeneous in x3. Data for the terms in the transport equation for E 0 are given for IHL in Grötzbach [21] and Wörner et al. [49] and for RBC in Wörner [48]. IHL, Ra=1E9,Pr=7.0

6

*1E-2 20

RMS of velocity profiles

RMS of velocity profiles

*1E-2 8

U1_RMS U2_RMS U3_RMS

4 2 0 1.0

1.2

1.4

X3

1.6

1.8

2.0

RBC, Ra= 6.3E5, Pr=0.71

15 10 U1_RMS U2_RMS U3_RMS

5 0 1.0

1.2

1.4

X3

Fig. 2 Vertical profiles of RMS values of the velocity fluctuations in IHL and RBC

1.6

1.8

2.0

138

Flow Turbulence Combust (2007) 79:133–154

*1E-2 4

IHL, Pr=7.0

Turbulent kinetic energy

Turbulent kinetic energy

*1E-2 0.5 0.4 0.3 0.2

TKE Ra=5E6 TKE Ra=1E7 TKE Ra=1E8 TKE Ra=1E9

0.1 0.0 1.0

1.2

1.4

X3

1.6

1.8

2.0

RBC

3 2 TKE Ra=6.3E5, Pr=0.71

1 0 1.0

1.2

1.4

X3

1.6

1.8

2.0

Fig. 3 Vertical profiles of TKE evaluated from the DNS data of IHL and RBC

2.2 Modeling requirements This section aims at the investigation of the performance of the standard gradient-diffusion model for the turbulent energy diffusion DE,t in buoyant horizontal fluid layers. According to Launder and Spalding [29] the model is as follows:    @  0 0 @ ν t @E 0 0 ; with j ¼ 1; 2; 3 ð1Þ uj E þ uj p0  DE;t ¼  @xj @xj σk @xj Here, s k is the turbulent Prandtl number of E 0 . In most codes this model is used with s k =1.0. Some problems of the model have been discussed in the introduction. One of the major problems is that it does not account for the counter gradient energy transport which may occur in buoyant flows, see e.g. Lumley et al. [34], Moeng and Wyngaard [36], and for IHL in Wörner et al. [49]. The DNS results for DE,t are given in Fig. 4. The problem of counter gradient energy fluxes occurs in the IHL over wide areas, e.g. in the area between x3 =1.7, which is about the position of the energy maximum, see Fig. 3, and x3 =1.94; in this range the turbulent energy flux is negative. This means, energy is transferred downwards towards the energy maximum, against the energy gradient. A counter gradient upward flux occurs in the area below mid-plane. Consistently, the comparison between the DNS results for DE,t and its gradient-diffusion model as given in (1) depicts the inadequacy of the model in this flow *1E-2 9

IHL, Ra=1E7, Pr=7

Turbulent diffusion

Turbulent diffusion

*1E-2 0.4

DNS Grad model

0.2

0.0

-0.2 1.0

1.2

1.4

X3

1.6

1.8

2.0

RBC, Ra=6.3E5, Pr=0.71 DNS Grad model

6 3 0 -3 -6 1.0

1.2

1.4

X3

1.6

1.8

Fig. 4 Vertical profiles of DE,t and its model as in (1) analyzed from the DNS data of IHL and RBC

2.0

Flow Turbulence Combust (2007) 79:133–154

*1E-3 0.2



Closure terms

0.0 -0.1 -0.2 -0.3 1.0

*1E-3 1.50

IHL, Ra=1E7, Pr = 7.0

0.1

Closure terms

139

1.2

1.4

X3

1.6

1.8

2.0

RBC, Ra=6.3E5, Pr=0.71

0.75 0.00 -0.75 -1.50 1.0

1.2

1.4

X3

1.6

1.8

2.0

Fig. 5 Vertical profiles of the closure terms of DE,t as in (1) analyzed from the DNS data of IHL and RBC

type; thus its improvement is inevitable. In RBC the deficiency is considerable apart from walls, but it is completely unacceptable in IHL. This may be regarded as one of the reasons that the standard E 0  e0 type models fail completely to describe thermally stratified flow types like IHL as it was experienced by Dinh and Nourgaliev [10]. Therefore, in the present work higher importance is given to the modeling of the turbulent diffusion in IHL, whereas the resulting model will be tested in RBC to ensure that the status quo is at least kept. 0 0 The gradient-diffusion model for DE,t assumes that uj p0 and uj E 0 can be modeled together. On the contrary, Fig. 5 shows the terms have different behavior and importance in IHL and RBC. In IHL u03 E 0 has larger values away from walls, whereas in RBC the values are lower than that of u03 p0 . In DE,t not the size of the terms is of interest but their vertical derivative, but this will behave similarly. Therefore, the contributions of both correlations in the two flows are even qualitatively quite different. Reasons for the differences in the importance of the pressure term are discussed in the light of the acting intermittent transport mechanisms in Wörner and Grötzbach [50], in which the deceleration or acceleration of the upward and downward moving 0 plumes is found to be important. Figure 5 also reveals that u3 p0 and its derivative is significant close to the upper and lower walls in both flows. Based on a DNS study of the effects of shear 0 on turbulent RBC, Domaradzki and Metcalfe [13] had also explained the importance of u3 p0 . They indicated a need of careful modeling of this term. Considering these results, an extended version of the gradient-diffusion model for DE,t will be derived. 3 Modeling of Turbulent Energy Diffusion 3.1 Modeling approach One approximation for linking the two contributions to DE,t was developed by Lumley [33]. 0 He showed that uj p0 ¼ 1=5u0j E0 holds in homogeneous turbulence. Checking this relation against the DNS results in Fig. 5 shows, in the anisotropic RBC there is proportionality between both correlations apart from walls, but the proportionality factor is greater than 5 instead of 1/5, so, in contrast to Lumley’s solution, the pressure term is the dominant one. For IHL the triple correlation is indeed the dominant one but, other than Lumley states, there is no proportionality between the two terms. The DNS values for E 0 in RBC in Fig. 3 depict that apart from the thermal boundary layers they are roughly constant, which means, the energy distribution is homogeneous; consistently the turbulent diffusion is constant over wide areas, Fig. 4. However, in IHL the vertical distribution of E0 is at these Ra strongly

140

Flow Turbulence Combust (2007) 79:133–154

non-constant, and consequently the turbulent energy transfer is hardly locally homogeneous. This may explain why the Lumley approximation fails in IHL. Due to these results 0 0 and those from the literature discussed above, in the present study uj E 0 and uj p0 are modeled separately. Afterwards the models will be coupled. A composed model may have a possibility to predict to some extent counter gradient energy fluxes. 0 0 The following assumptions are used to obtain separate RANS models for uj E 0 and uj p0 : (A1) (A2)

The flow types are horizontally homogeneous and there is no horizontal mean flow so that these are shear free in the current averaging method. The cross-correlations of the velocity fluctuations are smaller than their auto0 0 0 correlations i.e. ui uj  uj 2 for i≠j, with i, j=1, 2, 3. 0

3.2 Derivation of a RANS model for uj E 0 0

The development of the model for uj E 0 starts from a non-dimensional form of the transport 0 equation for uj E 0. Then, some of the closure terms in this equation are approximated according to the available literature. 0 The transport equation for uj E 0 is as follows, see e.g. Hanjalić and Launder [24] or Moeng and Wyngaard [36]: 9 ! ! 0 0 0 0 0 0 > @uj E 0 @uj E 0 @u u > @u u 0 0 j l > > ¼  ul þ ui uj i l þ E 0 > > @t @xl @xl @xl > > > > |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} > > Convection Production due to Rey: stresses and its gradient > > > ! >   > 0 0 > 0 > @u u E @u @u 0 0 0 0 j j l > i 0 >  E ul þ uj ui ul  > > > @xl @xl @xl > > |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} > |fflfflfflfflfflfflffl{zfflfflfflfflfflfflffl} > > Production due to mean shear > Turbulent transport > > >   > > Ra 02 > 0 0 0 > þ 2 E T þ uj T δj3 > = Re Pr |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ð2Þ Buoyancy > 0 1> > > > > > ! > B ( 2 0 0) C 0 > 0 0 0 > B 1 @ uj E C 1 > 0 @ui @ui 0 @ui @uj > C > 2u þB  þ u > j i 2 B Re C Re @x @x @x @x @x > l l l l @ A> l > > |fflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} > > > > Mol: diff Dissipation ðDisÞ > > 8 9 > ! > 0 0 0 0 > < @p0 = 0 > @u u @p u u > j l j l > 0 0 >  E δjl  p δjl : > > : @xj ; @xl @xl > > |fflfflfflfflfflffl{zfflfflfflfflfflffl} > > > |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Pressure transport ; Pressure term

In (2) all terms are closure terms except the convection and the production due to Reynolds stresses and its gradient. Due to assumptionnA1, the.convection and production due to mean o 0 shear vanish. The molecular diffusion 1=Re @ 2 uj E 0 @x2l can be neglected at high Re. The 0 turbulent diffusion of uj E 0 consists of the pressure transport and of the turbulent transport, 0

Turbulent diffusion of uj E 0 ¼ 

 @  0 0 0 0 0 uj ul E þ p0 uj ul δjl : @xl

ð3Þ

Flow Turbulence Combust (2007) 79:133–154

141

Here δjl denotes the Kronecker delta. The fourth-order velocity fluctuation correlations are approximated according to Millionshtchikov [35], i.e. when the triple correlations are small and their distributions do not differ substantially from those of a Gaussian one, the fourthorder correlations may be approximated in terms of second-order correlations. A proof of this simplification is available in Monin and Yaglom [37]. This gives,       0 0 0 0 0 0 0 0 0 0 uj ul ui 2  ul uj ui 2 þ 2 ui uj ui ul : ð4Þ Considering assumption A2 results in, ! 0 0 0 0 @u u 0 0 @ui ul j l  ui uj þ E0 @xl @xl

0

0

@uj ul E0 @xl

!

0

0

 uj 2

0 @uj 2 0 @E δjl  uj 2 δjl : @xl @xl

ð5Þ

Index j appears in each term of the lhs only once; thus, the Einstein-summation is not applicable to this index in (5) and following. the approach of Hanjalić and  . According  Launder [24] the pressure strain term p0 @u0j u0l @xl is modeled as in Rotta [41] and the dissipation term as in Zeman and Lumley [51], !) 8 ! 9 0 0 0 < = @u @u u u0 j E 0 1 0 @u @ui 0 @u j j l δjl  C 0 1  þ ui i þ p0 2uj i : : ; Re @xl @xl @xl @xl @xl τ |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} (

0

0

0

ð6Þ

Dissipation

Here, C 0 1 is a coefficient and t is a time scale, t ¼ E 0 e0 , e.g. from Launder et al. [30]. The right hand side of this equation is known as the relaxation term, i.e. triple correlation divided by time scale, see e.g. Lumley et
al. [34].  0 The contribution of buoyancy Ra Re2 Pr E0 T 0 þ uj 2 T 0 δj3 is only along j=3. Figures 2 and 3 indicate that indeed most of the turbulent kinetic energy E 0 is associated with the vertical velocity component in IHL. Thus, the buoyancy term may be simplified to 

 0 0 E0 T 0 þ uj 2 T 0 δj3  2uj 2 T 0 δj3 :

ð7Þ

This assumption may not be sufficient in case of RBC, but in the present case, priority is given to achieve a better modeling for IHL. Considering fully developed flow, that means the steady state, using assumptions A1 and A2, and introducing (5) and (7) in (2) gives @p0 uj @E 0 Ra  0 2 0  δjl þ 2 δjl  2uj T δj3  Re Pr @xl @xl ( !) 8 ! 9 0 02 02 0 0 0 < = 0 @u @u @u 1 @p 0 @u @ui 0 @u 0 j j j 2uj i  uj 2 δjl : þ ui i δjl þ E0  p0 : @xl ; Re @xl @xl @xl @xl @xj @xl 02

0

0  uj 2

0 0 According to Weinstock [47] the terms uj 2 @uj 2 =@xl δjl and E 0 @p0 @xj should be proportional to each other. Indeed, the analysis for RBC in the Appendix shows a proportionality factor of order −1 so that the sum of both terms may be neglected here. Thus, introduction of (6),

142

Flow Turbulence Combust (2007) 79:133–154

rearrangement and calculation of the time scale with t ¼ E 0 e0 results in the required model for the triple correlation, 0 1 0

uj E 0  C1

0 E0 B Ra 0 0 2 @E B u δjl þ 2 2 uj 2 T 0 δj3  B j @xl Re Pr |fflfflffl{zfflfflffl} e0 @ |fflfflfflfflffl{zfflfflfflffl ffl} Buoyancy

Grad: Approx:

0

@p0 uj 2 @xl |fflffl{zfflffl}

C C δjl C: A

ð8Þ

Press: Transport

Here C1  1 C 0 1 is a coefficient with 0.04≤C1 ≤0.17 (Chandra [4]). The range of this coefficient is expected to reduce in IHL at very high Ra as a result of homogeneity away from the walls as in RBC. 0 Equation (8) is the basis for a new RANS model for uj E 0. The first term is a gradientdiffusion approximation including the anisotropy of the velocity fluctuations. The second term is the buoyancy contribution. The last term is the pressure-transport. Thus, this model basis (8) compared to the standard model by Launder [32] has two additional terms, which up to now do not depend on the gradient of the energy, along with an anisotropic form of the common gradient-diffusion approximation. All three still contain closure terms. These will be approximated or discussed in Section 3.4. 0

3.3 Modification of a RANS model for uj p0 0

In order to obtain a modified model for uj p0, this correlation is represented as in Donaldson [14] by the following redistribution approach, ! 0 0 @u u 1=2 0 j l 0 : ð9Þ uj p0  C 2 E 0 l @xl 0

C2 is a coefficient and l is a length scale. Using the Rotta [41] dissipation model for high . 3=2 e0 . Reynolds numbers l is defined as a function of E 0 and e0 , i.e. l ¼ E0 0 Daly and [6] indicate that C2 is a function of the turbulent Reynolds number .
Harlow  2 νe0 . This is usually employed in the near-wall damping functions. It holds Ret ¼ E 0 obviously also in IHL because Wörner et al. [49] found that some model coefficients in the transport equations for E 0 and for the turbulent heat fluxes have to be increased by large Fig. 6 Coefficient C 0 2 at x3 = 1.18

IHL, Pr=7.0

100

'

C2 in '

2 Linear CFit

Ra=5E6 Ra=1E7

'

C2 10 α=0.8 (approx)

Ra=1E8

1 0.01

0.1 Ret

1

Flow Turbulence Combust (2007) 79:133–154

143

factors in this convection type. Figure 6 indeed indicates that C 0 2 ¼ C2 Reat with a =0.8 and 1.5≤C2 ≤7, Chandra [4]. Moreover, using assumption A2 in (9) results in the modified RANS model, ! 02 2 C2 E 0 @uj 0 0 δjl : uj p   α ð10Þ Ret e0 @xl This model does also not involve a gradient of the kinetic energy. In flow types like RBC, 0 where the contribution of uj E0 to DE,t is small, an even better model for the derivative of 0 0 uj p may be needed. A more extended modeling which accounts also for the anisotropy of the velocity field is discussed in Chandra [4] and in Chandra and Grötzbach [5]. This includes the concept of an anisotropic eddy viscosity as in Durbin [15]. 3.4 Derivation of a RANS model for DE,t 0

0

Having deduced the separate RANS models for uj E0 and uj p0 , the model for DE,t is derived by adding the models given in (8) and (10), 0

0

uj

E0

0

þ uj

p0

2 E 0 0 @E 0 E0 Ra 0 2 0 E0 @p0 uj C2  C1 0 uj 2 δjl þ 2C1 0 u T δ C δjl  a j3 1 2 Pr j 0 @xl @x Re Re e e e l t |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflffl}

Buoyancy term

Grad: Approx:

! 02 2 E0 @uj δjl : e0 @xl

Pressure term

ð11Þ This model includes the contribution of buoyancy and a higher-order pressure term. The first and last term on the rhs without their coefficients follow the approximation given in (4). In this equation the buoyancy and pressure terms need to be modeled. Among these terms, the pressure term is problematic. Its importance and behavior in IHL and RBC may be expected to be different due to the different mechanisms in IHL and RBC which lead to pressure and velocity fluctuations. 3.4.1 Modeling of the higher-order pressure term 0

0

*1E-3 0.2

*1E-3 2

IHL, Ra=1E7, Pr=7

Closure terms in RANS model

Closure terms in RANS model

The contributions of the buoyancy and the pressure terms to uj E 0 þ uj p0 analyzed in Fig. 7 indicate their different importance in IHL and RBC. There is only a small contribution of the pressure term in IHL compared to the buoyancy term. In contrast, the pressure term is not negligible in RBC.

Buoyancy_term Pressure_term +

0.1 0.0 -0.1 -0.2 -0.3 1.0

1.2

1.4

X3

1.6

1.8

2.0

RBC, Ra=6.3E5, Pr=0.71

1 0 -1 -2 1.0

Buoyancy term Pressure term +

1.2

1.4

X3

1.6

1.8

2.0

Fig. 7 Vertical profiles of terms in (11) evaluated from the DNS data of IHL with C1 =0.17 and RBC with C1 =0.08

144

Flow Turbulence Combust (2007) 79:133–154

In the present work the primary goal is to obtain an improved model for DE,t in IHL, i.e. for flows with local stable stratification. Since there is no reliable model available for the pressure term and its contribution is small in IHL, it will be neglected in the final form of the model for DE,t which results in, ! 0 0 0 0 0 2 @u 2 E E Ra C E 0 0 0 2 @E 02 j 2 δjl : ð12Þ uj E0 þ uj p0  C1 0 uj δjl þ 2C1 0 uj T 0 δj3  a @xl Ret e e Re2 Pr e0 @xl |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Buoyancy term

As a result of this simplification some deviations between DE,t and its modeled values may be expected in RBC, or some dependencies of the model coefficients on flow parameters may occur. 3.4.2 Modeling of the buoyancy term The buoyancy term from (12) was found in Fig. 7 to be important in the partially stably stratified IHL. Thus, an improved modeling for this term should be achieved. It contains the 0 unknown triple correlation uj 2 T 0 . Using the Daly and Harlow [6] approximation and introducing the assumption A2 results in the following model, ! 0 0 0 @uj 2 E0 02 0 2 @uj T 0 0 0 : ð13Þ uj T  Cθ 0 2uj þ uj T @xj @xj e Here, Cq is a coefficient. In most of the literature the values of this coefficient are assumed to be constant. On the contrary, Dol et al. [11] had shown that it is not a constant in a buoyant flow. Similarly, Wörner et al. [49] found that the coefficient of a turbulent diffusion model in the transport equation for the heat flux had to be considerably increased in IHL. Thus, the model needs improvement. Therefore an extended version of the Daly and Harlow (DH) model is derived. 0 In order to investigate u3 2 T 0 in a horizontal channel in more detail, the terms in its transport equation are analyzed using the DNS data. The non-dimensional form of the transport equation for u03 2 T 0 as given in Dol [12] is as follows, ! 0 0 @u3 2 T 0 @u 2 T 0 ¼  uk 3 þ @t @xk |fflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflffl} 0

Convection

! 0 0 0 @uk T 0 @u u 0 þ 2u3 T 0 3 k @xk @xk |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} 0

u3 2

Production by Reynolds stress and turbulent heat fluxes ðProSÞ

1

! 0 0 B 0 2 0 @T @u3 2 uk T 0 @u3 C Ra  0 0 2  0 0 0 C B u u þ 2u u T þ2 2 u T 3 k @ 3 k @xk Re @xk A @xk Pr 3 |fflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflffl} |fflfflfflfflffl{zfflfflfflfflffl} |fflfflfflfflfflfflfflffl ffl {zfflfflfflfflfflfflfflffl ffl } Buoyancy ðProBÞ ProT |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} Tur: transport ðTurbTÞ Prod: due to mean Temp: and mean shear

0

1

9 > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > > =

C B C B C >: B C > B  0 (  )! 0 0 C > B 0 0 0 0 1 @ @u 1 @T @u3 T @p u3 T > 0 0 C > B >  þ 2u3 T 0 3 þ u3 2 þ 2B p0 Cδk3 > > C > B Re @xk @xk Pr @xk @xk @xk > > |fflfflffl ffl {zfflfflffl ffl } C B |fflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflffl} |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} > C > B > B Press: strain Press: transport C > Molecular terms ðMÞ > A > @ > > > > > ðDputÞ ðPdut Þ > > |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} > > > > Pressure terms ðPÞ > > > ! ( ) ! >    0 2 > 0 > 0 > 2 @u3 1 @u @T 0 > 3 0 > :  T þ 1þ u3 > > Re Pr @xk @xk @xk > > > > |fflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl{zfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflfflffl} ; Dissipative terms ðDÞ

ð14Þ

Flow Turbulence Combust (2007) 79:133–154

*1E-3 2

Terms in transport eqn for

Fig. 8 Vertical profiles of the terms in the transport equation 0 for u3 2 T 0 as in (14) for IHL with Ra=107

145

IHL, Ra = 1E7, Pr = 7.0

ProS -ProT -TurbT ProB Mol (M) Pdut -Dput -Diss (D) Budget

0

-2

-4 1.0

1.2

1.4

1.6

1.8

2.0

X3

pffiffiffiffiffiffi


 Using the present scaling results in Re ¼ Gr and Ra Re2 Pr ¼ 1. In this equation all terms other than the convection and the production due to Reynolds stresses and turbulent heat fluxes are closure terms. In the absence of mean shear the production is only due to the mean temperature gradient (ProT). The analyzed vertical profiles of all terms in this transport equation are shown in Figs. 8 and 9 for IHL and RBC. These also include the budget or out of balance of this equation which is calculated using all terms. It is smaller in comparison to most of the other terms and thus verifies the formal equation and its numerical representation. There is a second form of this transport equation in literature by Sander [42]. The values of the out of balance term calculated from his equation are found to be considerably larger than that of the equation by Dol [12]. Further investigations indicate towards a formal problem in Sanders equation, especially in his molecular and dissipative terms. Figures 8 and 9 show that none of the terms can simply be classified in both flow types, e.g. the production due to the Reynolds stresses and turbulent heat fluxes has both positive as well as negative values in RBC and the dissipation terms are not everywhere positive. The classification of these terms depending on the involved derivatives is only formal and not based on their real action. As the terms can hardly be classified, in the modeling of 0 0 u3 2 T 0 it is difficult to assume local equilibrium. One practical way to model u3 2 T 0 is to identify and use those terms, which are of higher importance.

*1E-3 6

Terms in transport eqn for

Fig. 9 Vertical profiles of the terms in the transport equation 0 for u3 2 T 0 as in (14) for RBC

RBC, Ra = 6.3E5, Pr=0.71

ProS -ProT -TurbT ProB Mol (M) Pdut -Dput -Diss (D) Budget

4 2 0 -2 -4 -6 1.0

1.2

1.4

X3

1.6

1.8

2.0

146

Flow Turbulence Combust (2007) 79:133–154

The production due to Reynolds stresses and turbulent heat fluxes (ProS) and the turbulent transport (TurbT) has higher significance in RBC than in IHL. The production due to the temperature gradient (ProT) is important close to the walls in both flows, see also Fig. 1. The contribution of the buoyancy term (ProB) is comparable to the dissipative terms (D) in IHL. The higher values of the pressure-transport (Dput) and of the pressure-strain term (Pdut) close to the walls in RBC can be accounted to the presence of a local region of high values of the pressure fluctuations as well as of the turbulent heat flux close to the walls. The molecular terms (M) need attention close to the walls due to the presence of strong viscous effects. The appearance of Pr in the denominator of the molecular and dissipative terms depict that their contribution will be enhanced in liquid metals. Therefore, these observations indicate that in addition to the production due to Reynolds stresses (ProS), turbulent transport (TurbT) and dissipative terms (D), the production due to the temperature gradient (ProT), buoyancy contribution (ProB), the molecular terms (M) need also to be included in a model for u03 2 T 0. The DH model for u03 2 T 0 as in (13) can be derived from (14) by balancing the production due to the Reynolds stresses and turbulent heat fluxes (ProS), the turbulent transport (TurbT) and the dissipation terms (D). On the other hand, there are additional terms, which 0 are important in the different flows. In order to obtain an extended model for u3 2 T 0, its transport equation (14) is used in which similar assumptions and approximations are employed as in the velocity fluctuation triple correlations. Further at high Ra with moderate Pr the molecular terms (M) are not considered. Assuming fully developed steady convection, (14) results in the following Daly and Harlow Extended (DHE) model: ! 0 0 0 E0 @u3 2 Ra 0 0 2 0 0 2 0 2 @u3 T 0 0 3 @T 0 0 ð15Þ u3 T  Cθ1 0 2u3 þ u3 T þ u3  2 2 u3 T : @x3 Re Pr @x3 @x3 e Here, C 0 θ1  1=C is a coefficient. . Analogous to the results of Dol et al. [11] and considering Chandra [4], C 0 θ1  Cθ1 Reβt has been used with β=0.52 and Cq1  0:25. The DH model (13) contains only the first two terms on the rhs of (15). This DHE model also includes the production due to the mean temperature gradient and the contribution of 0 0 buoyancy. The last two terms involve the higher-order correlations u3 3 and u3 T 0 2 . The 0 3 0 closure term u3 may be modeled according to Launder [32] and the triple correlation u3 T 0 2 according to an extended model derived by Otić et al. [39]. 3.4.3 Extended RANS model for DE,t 0

In order to obtain the extended model for DE,t the DHE model for u3 2 T 0 (15) is introduced in (12). Taking the partial derivative along j=3 results in, 9 8 E0 0 2 @E0 > > > > C u þ > > 1 0 3 > > > > @x3 e > > > > 1 0 > > 0 0 2 > > 0 > > @u T @u @T 0 2 0 0 3 > > 3 3 0 > > 2u þ u T þ u > > C B 0 0 3 3 3 = < o C E Ra E C θ1 B @x @x @x @ n 0 0 @ 3 3 3 0 0 2C C B 1 u3 E þ u3 p   β 2 0 0 A >: Ret e Re Pr e @ @x3 @x3 > Ra  0 0 2  > > > >  2 u T > > 3 > > Re2 Pr > > > > ! > > > > 0 2 > > 0 @u 2 > > C E 2 > > 3 > >þ ; : 0 Reαt @x e 3

ð16Þ

Flow Turbulence Combust (2007) 79:133–154

147

In this extended RANS model (16) for DE,t the same values of the coefficients and parameters as in (8), (10) and (15) will be used; for details see Appendix B in Chandra [4]. Applications of this model for DE,t require an additional modeled transport equation for u03 2, 0 see e.g. Launder et al. [30], and a closure for u3 T 0 2, see e.g. Otić et al. [39]. In other words, this 0 model extends the E0  e0 model to a E 0  e0  u3 2 model, which is a 3-equation model.

4 Validation of Proposed Models In order to validate the new or modified closure models, these will be analyzed from the DNS data of IHL and RBC. The following RANS models will be considered by using 0 model and flow specific coefficients: The model basis for uj E 0 as given in (8), the model for 0 0 uj p0 as in (10), and the Daly and Harlow Extended (DHE) model for uj 2 T 0 as in (15). In the validation of the complete model for DE,t as given in (16) a unique set of coefficients is applied for both flow types. Here, the horizontal plane and time-averaged variables depend only on x3. Therefore, the models are validated only in this direction. There are data for two different Rayleigh numbers used in the validation of the complete turbulent diffusion model for IHL, because the statistical flow data are at these moderate Ra still not yet self-similar, see e.g. Fig. 3, whereas for RBC it is sufficient to consider only one Ra number because Ra is large enough so that the statistical data are in the chosen scaling independent on Ra, except for the change of the boundary layer thickness. 0

4.1 Validation of the RANS model for u3 E 0 0

In the validation of the model basis for u3 E 0 from (8), the coefficient C1 is set to 0.17 for IHL and to 0.08 for RBC. The vertical profiles of u03 E 0 and its ideal model, which contains still closure terms on the rhs, are analyzed from the DNS data and compared in Fig. 10. In 0 IHL the model basis shows an acceptable agreement with u3 E0 . In RBC the RANS model over-predicts the DNS values. However, the modeled values have roughly the required qualitative distribution. 0 0 Figure 5 indicates that the derivative of u3 E 0 is greater than that of u3 p0 in most of the 0 central region ð 1:25  x3  1:8Þ in IHL. Therefore, u3 E 0 is more significant in this flow. 0 Further, the derivative of u3 E 0 is very small close to walls. Consequently, its contribution to IHL, Ra=1E7, Pr=7.0

*1E-3 1.0

DNS_ Mod_

0.0

Triple correlation

Triple correlation

*1E-3 0.1

-0.1 -0.2 -0.3 1.0

1.2

1.4

X3 0

1.6

1.8

2.0

RBC, Ra=6.3E5, Pr=0.71 DNS_ Mod_

0.5 0.0 -0.5 -1.0 1.0

1.2

1.4

X3

1.6

1.8

2.0

Fig. 10 Vertical profiles of u3 E0 and its model basis as in (8) analyzed from the DNS data of IHL with C1 = 0.17 and of RBC with C1 =0.08

Flow Turbulence Combust (2007) 79:133–154

*1E-3 0.2

IHL, Ra=1E7, Pr=7.0

*1E-3 2

DNS_ Mod_

0.1

0.0

-0.1 1.0

1.2

1.4

X3

1.6

1.8

2.0

Vel-press fluctuation correlation

Vel-press fluctuation correlation

148

RBC, Ra=6.3E5, Pr=0.71 DNS_ Mod_

1 0 -1 -2 1.0

1.2

1.4

X3

1.6

1.8

2.0

0

Fig. 11 Vertical profiles of u3 p0 and its model as in (10) analyzed from the DNS data of IHL with C2 =1.5 and of RBC with C2 =3.0, a=0.8

DE,t is not significant in the near-wall region. For RBC the opposite behavior is found, 0 Fig. 5. Indeed, the model also reveals that the derivative of u3 E0 is smaller than that of u3 p0 0 0 in RBC. Therefore, also in the model the contribution of u3 E to the turbulent diffusion is less significant in RBC. 0

0

4.2 Validation of the RANS model for u3 p0 0

In the validation of the modified model for u3 p0 from (10), the coefficient C2 is set to 1.5 for 0 IHL and to 3.0 for RBC and the parameter a =0.8. The vertical profiles of u3 p0 and its model analyzed from the DNS data in Fig. 11 show an acceptable qualitative and 0 quantitative agreement in both flows. Moreover, u3 p0 is very significant close to the walls in 0 both flows. Figure 5 indicates that the derivative of u3 p0 is smaller than the derivative of 0 0 0 u3 E in most of the central region ð 1:25  x3  1:8Þ in IHL. It implies, u3 p0 has a smaller 0 contribution to DE,t in this region for this flow. As a consequence, the model for u3 p0 will not play a significant role in this region, whereas it has already been noted that in RBC the 0 0 0 derivative of u3 p0 is greater than that of u3 E 0 . It implies, that the model for u3 p0 gives the expected higher contribution to DE,t in this flow. *1E-3 0.6

IHL, Ra = 1E7, Pr=7.0 DNS Daly and Harlow (DH) Daly and Harlow Extended (DHE)





0.3 0.0 -0.3 -0.6 -0.9 1.0

*1E-3 2

1.2

1.4 0

X3

1.6

1.8

2.0

RBC, Ra = 6.3E5, Pr=0.71

1 0 DNS Daly and Harlow (DH) Daly and Harlow Extended (DHE)

-1 -2 1.0

1.2

1.4

X3

1.6

1.8

2.0

Fig. 12 Vertical profiles of u3 2 T 0 and its models as in (13) and (15) analyzed from the DNS data of IHL and RBC

Flow Turbulence Combust (2007) 79:133–154

*1E-2 0.4

IHL, Ra=1E7, Pr=7.0 DNS Grad Model Extended Model

0.2

Turbulent diffusion

Turbulent diffusion

*1E-2 0.4

149

0.0

-0.2 1.0

1.2

1.4

X3

1.6

1.8

DNS Grad Model Extended Model

0.2

0.0

-0.2 1.0

2.0

IHL, Ra=1E8, Pr=7.0

1.2

1.4

1.6

1.8

2.0

X3

Fig. 13 Vertical profiles of DE,t and its models as in (1) and (16) analyzed from the DNS data of IHL with Ra=107 and 108

0

4.3 Validation of the DHE model for u3 2 T 0 0

In the validation of the DHE model for u3 2 T 0 from (15), the coefficients Cq1 ¼ 0:25 and 0 β=0.52 are used in IHL and in RBC. The comparisons between u3 2 T 0 and its models in IHL and RBC in Fig. 12 indicate that the DH model from (13) with Cq  0:11 needs improvement in both flow types because it produces not only quantitatively wrong data, but produces even a qualitatively wrong vertical distribution. The Daly and Harlow Extended (DHE) model from (15) shows a considerable qualitative improvement in IHL except close to the lower wall. In RBC the DHE model shows an acceptable agreement close to the walls.

4.4 Validation of the RANS model for DE,t In the validation of the extended model for DE,t from (16), a single set of coefficients and parameters is used: C1 =0.17, Cq 1 ¼ 0:25, C2 =2, a =0.8 and β=0.52. These data belong to

*1E-2 9

Turbulent Diffusion

Fig. 14 Vertical profiles of DE,t and its models as in (1) and (16) analyzed from the DNS data of RBC

RBC, Ra=6.3E5, Pr=0.71 DNS Grad Model Extended Model

6 3 0 -3 -6 1.0

1.2

1.4

1.6 X3

1.8

2.0

150

Flow Turbulence Combust (2007) 79:133–154

their respective regions as explained in the previous sections and subsections. The comparison between DE,t and its models in Fig. 13 reveals significant improvement by the extended model compared to the standard gradient model from (1) in IHL at both Rayleigh numbers. Especially in those areas is a qualitative improvement achieved, in which according to Figs. 3 and 4 a counter gradient energy flux occurs, like below mid-plane and between x3 =1.7 and x3 =1.9. There, the composed model indeed predicts the required counter gradient energy flux. In RBC of air the model for DE,t gives some qualitative improvement, especially as regards to the position of the near-wall change of sign, compared to the standard gradient model as shown in Fig. 14. Nevertheless, in the core of the fluid layer the model should further be improved. Quantitatively, this model gives quite improved results in the near-wall region compared to the standard model.

5 Conclusions Buoyant flows are a challenge in turbulence modeling, especially when areas with stable stratification are involved. Peculiarities are the consequences of the large anisotropy of the turbulence field on the turbulent fluxes and the tendency to develop areas with counter gradient fluxes of several turbulence quantities. Studies in meteorology e.g. by Moeng and Wyngaard [36] indicate the need of improving the standard model for the turbulent diffusion of turbulent kinetic energy. The inadequacy of this gradient-diffusion model has already been observed by Dinh and Nourgaliev [10] in numerically investigating the cooling of a horizontal fluid layer with a combination of mostly stable and a thin unstable stratified layer. Analyses of DNS data of the two closure terms, the velocity-pressure fluctuation correlation and the velocity-fluctuation triple correlation, reveal their different importance in an Internally Heated Layer (IHL) and in Rayleigh–Bénard Convection (RBC). Some detailed underlying mechanisms of these terms are explained by Wörner and Grötzbach [50]. Basing on these studies separate RANS models for these two closure terms are derived by starting from existing approximations. For the velocity-pressure fluctuation correlation a simple model by Donaldson [14] is adopted, in which a coefficient is introduced which depends on the local turbulent Reynolds number. The modeling for the velocity-fluctuation triple correlation starts from its transport 0 equation. Therein, a buoyancy term containing the triple correlation uj 2 T 0 appears as one of the important closure terms. In a horizontal fluid layer it is non zero along the vertical direction j=3. The Daly and Harlow [6] model for this term has already been reported to be inadequate in buoyant flows e.g. by Dol et al. [11] and Wörner et al. [49]. All terms in the transport equation for this triple correlation are analyzed from DNS data. From this analysis a buoyancy-extended version of the Daly and Harlow model is deduced and incorporated in the model for the velocity-fluctuation triple correlation. Coupling both models for the closure terms in the turbulent energy diffusion results in an extended RANS model, in which three additional closure terms occur: The triple 0 0 correlations u3 3 and u3 T 0 2 can be calculated using models by Launder [32] and by Otić et 0 al. [39], respectively. The third term, u3 2 , should be calculated by an additional modeled transport equation e.g. as in Launder et al. [30]. Thus, this diffusion model extends the 0 E 0  e0 model to a E0  e0  u3 2 3-equation model. It allows for a better representation of the consequences of the large anisotropy which may occur in buoyant flows. In addition, the extended turbulent diffusion model is not exclusively based on gradient-diffusion assumptions; therefore, it also allows predicting counter gradient energy fluxes.

Flow Turbulence Combust (2007) 79:133–154

151

The models for all closure terms and the coupled model for the turbulent diffusion of energy are validated using DNS data. Each of the models was found to represent in a roughly acceptable manner its contributions especially in those areas in which the corresponding terms are relevant in both flow types. Also the shift of the relative importance between the pressure and triple correlation term in the energy diffusion between IHL and RBC is well reproduced. This shows the acceptability of the models in both flow types. A significant improvement in the prediction of the turbulent diffusion of energy in IHL is achieved, as the extended model well accounts for the strong counter gradient energy fluxes. In case of RBC the diffusion model shows a slight improvement compared to the existing standard gradient-diffusion model. Another less general version of this diffusion model in which a modified form of the Daly and Harlow approximation has been used, also shows an acceptable predictive capability in IHL and RBC, see Chandra [4]. It has also been validated in RBC with liquid metals (Pr=0.025). So far one could expect that incorporating one of these extended models in current CFD codes could result in more accurate predictions of such buoyant flow types. Unfortunately, partially stable stratified fluid layers like the IHL tend to form also counter gradient turbulent heat fluxes which cannot be reproduced by the thermal gradient-diffusion assumption, which is nowadays implemented in most codes. The simplest models which are suited to reproduce counter gradient heat fluxes are algebraic heat flux models as in Launder [31] or in Otić and Grötzbach [38, 40]. Such models require for buoyant flows an additional modeled transport equation for the temperature variance T 0 2 . For more complicated flows like in liquid metal heat transfer specially suited models for the 0 dissipation of the temperature variances eT are required as developed in Otić and Grötzbach 0 [40] or an additional modeled transport equation for eT is required. For both transport equations variants are available for low Prandtl number fluids which were already applied in a CFD code to mixed convection, see Carteciano and Grötzbach [3], or which were recently developed for buoyant convection, see Otić et al. [39]. Thus, a suitable concept to calculate the turbulent heat transfer in buoyant fluid layers with partially stable stratification should at least be based on an algebraic heat flux model and an E 0  e0 model involving the here discussed extensions of the turbulent diffusion, so that the resulting model will be 0 0 0 based on the E0  e0  u3 2  T 0 2 or on the E 0  e0  u3 2  T 0 2  eT equations. Further improvements could be achieved by applying the Algebraic Stress Model concept, e.g. as in Durbin [15], also to the shear stress modeling so that the anisotropy occurring in buoyant flows is better accounted for. Acknowledgements The authors are thankful to Prof. Dr. T. Schulenberg, head IKET at FZK, Prof. Dr. H. Oertel Jr., head ISL at Universität Karlsruhe (TH), Prof. Dr. C. Günther, Prof. Dr. U. Müller and other colleagues for their valuable suggestions and personal support during the Ph.D. work of the first author. Special thanks go to Dr. I. Otić and Dr. M. Wörner for providing valuable DNS data for Rayleigh–Bénard convection and for the fruitful discussions with them.

Appendix .

0 0 Weinstock [47] proposed, the terms uj 2 @uj 2 @xl δjl and E0 @p0 @xj should be proportional to each other with a proportionality coefficient of magnitude greater than or equal to 1. Analysis of the corresponding vertical profiles for IHL and RBC for j=3 is given in Fig. 15. Indeed, for RBC it is found that both terms have an opposite behavior in the homogeneous core of the fluid layer where the turbulent kinetic energy is roughly constant, see Fig. 3. So, the proportionality factor is of order −1. As a result, the sum of both terms may be neglected in this flow type. Unfortunately, this proportionality does not hold in the IHL at

152

Flow Turbulence Combust (2007) 79:133–154

*1E-3 3

IHL, Ra = 1E7, Pr=7.0

0.08

Weinstock terms

Weinstock terms

*1E-3 0.12

0.04 0.00 -0.04 -0.08 -0.12 1.0

RBC, Ra=6.3E5, Pr=0.71

2 1 0 -1 -2

1.2

1.4

1.6 X3

1.8

2.0

-3 1.0

1.2

1.4

X3

1.6

1.8

2.0

Fig. 15 Vertical profiles of the terms in the Weinstock approximation in IHL and RBC

these low Rayleigh numbers. Nevertheless, a similar behavior may be expected also in IHL at high Ra in which the vertical distribution of the turbulent kinetic energy should be homogeneous in the core of the fluid layer as in RBC.

References 1. Canuto, V.M., Christensen-Dalsgaard, J.: Turbulence in astrophysics. Annu. Rev. Fluid Mech. 30, 167– 198 (1998) 2. Carteciano, L.N.: Entwicklung eines Turbulenzmodells für Auftriebsströmungen. Ph.D. thesis, University of Karlsruhe, Forschungszentrum Karlsruhe KfK 5775 (1996) 3. Carteciano, L.N., Grötzbach, G.: Validation of turbulence models in the computer code FLUTAN for a free hot sodium jet in different buoyancy flow regimes. Forschungszentrum Karlsruhe FZKA 6600 (2003) 4. Chandra, L.: A model for the turbulent diffusion of turbulent kinetic energy in natural convection. Ph.D. thesis, University of Karlsruhe, Forschungszentrum Karlsruhe FZKA 7158 (2005) 5. Chandra, L., Grötzbach, G.: Modeling the derivative of the pressure-velocity fluctuation correlation in turbulent natural convection. In: Proceedings of 9th AeSI CFD Symposium, 11–12th August, NAL, Bangalore, Paper No: CP19, (2006) 6. Daly, B.J., Harlow, F.H.: Transport equations in turbulence. Phys. Fluids 18, 2634–2649 (1970) 7. Davidson, L.: Second-order corrections of the k−ɛ model to account for non-isotropic effects due to buoyancy. Int. J. Heat Mass Transfer 33, 2599–2608 (1990) 8. Deardorff, J.W.: Three-dimensional numerical study of the height and mean structure of a heated planetary boundary layer. Boundary – Layer Meteorol. 7, 81–106 (1974) 9. Deardorff, J.W.: Stratocumulus-capped mixed layers derived from a three-dimensional model. Boundary – Layer Meteorol. 18, 495–527 (1980) 10. Dinh, T.N., Nourgaliev, R.R.: Turbulence modelling for large volumetrically heated liquid pools. Nucl. Eng. Des. 169, 131–150 (1997) 11. Dol, H.S., Hanjalić, K., Kenjereš, S.: A comparative assessment of the second-moment differential and algebraic models in turbulent natural convection. Int. J. Heat Fluid Flow 18, 4–14 (1997) 12. Dol, H.S.: Turbulence models for natural convection in side-heated enclosures. Ph.D. thesis, Technische Universiteit Delft (1998) 13. Domaradzki, J.A., Metcalfe, W.: Direct numerical simulations of the effect of shear on turbulent Rayleigh–Bénard convection. J. Fluid Mech. 193, 499–531 (1988) 14. Donaldson, C. duP.: A computer study of an analytical model of boundary layer transition. AIAA J. 7, 272–278 (1969)

Flow Turbulence Combust (2007) 79:133–154

153

15. Durbin, P.A.: Near wall turbulence modeling without damping functions. Theor. Comput. Fluid Dyn. 3, 3–13 (1991) 16. Dwyer, M.J., Patton, G.E., Shaw, R.H.: Turbulent kinetic energy budgets from a large eddy simulation of airflow above and within a forest canopy. Boundary – Layer Meteorol. 84, 23–43 (1997) 17. Grötzbach, G.: Spatial resolution requirements for numerical simulation of internally heated fluid layers. In: Taylor, C., Schrefler, B.A. (eds.) Numerical Methods in Laminar and Turbulent Flow, pp. 593–604. Pineridge Press, Wales, UK (1981) 18. Grötzbach, G.: Direct numerical simulation of the turbulent momentum and heat transfer in an internally heated fluid layer. In: Grigull, U., Hahne, E., Stephan, K., Straub, J. (eds.) Heat Transfer 1982, pp. 141–146. Hemisphere, Washington (1982) 19. Grötzbach, G.: Direct numerical simulation of laminar and turbulent Bénard convection. J. Fluid Mech. 119, 27–53 (1982) 20. Grötzbach, G.: Spatial resolution requirements for direct numerical simulation of the Rayleigh–Bénard convection. J. Comput. Phys. 49, 241–264 (1983) 21. Grötzbach, G.: Direct numerical and Large Eddy simulation of turbulent channel flows. In: Cheremisinoff, N.P. (ed.) Encyclopedia of Fluid Mechanics, Complex Flow Phenomena and Modeling, vol. 6, pp. 1337–1391. Gulf, Houston, TX (1987) 22. Grötzbach, G.: Turbulent heat transfer in an internally heated fluid layer. In: Iwasa, Y., Tamai, N., Wada, A. (eds.) Refined Flow Modelling and Turbulence Measurements, pp. 267–275. Universal Academy Press, Tokyo, Japan (1989) 23. Grötzbach, G.: Simulation of turbulent flow and heat transfer for selected problems of nuclear thermal hydraulics. In: Japan Atomic Energy Research Institute (ed.), The 1st International Conference on Supercomputing in Nuclear Applications (SNA’90), pp. 29–35, March 12–16, Nuclear Energy Data Center, Japan (1990) 24. Hanjalić, K., Launder, B.E.: A Reynolds stress model of turbulence and its application to thin shear flows. J. Fluid Mech. 52, 609–638 (1972) 25. Hanjalić, K., Kenjereš, S., Durst, F.: Natural convection in partitioned two-dimensional enclosures at higher Rayleigh numbers. Int. J. Heat Mass Transfer 39, 1407–1427 (1996) 26. Hanjalić, K.: Second-moment turbulence closures for CFD: needs and prospects. Int. J. Comput. Fluid Dyn. 12, 67–97 (1999) 27. Hartlep, T., Tilgner, A., Busse, F.H.: Large scale structures in Rayleigh–Bénard convection at high Rayleigh numbers. Phys. Rev. Lett. 91, 064501–064504 (2003) 28. Kulacki, F.A., Richards, D.E.: Natural convection in plane layers and cavities with volumetric energy sources. In: Kakac, S., Aung, W., Viskanta, R. (eds.) Natural Convection – Fundamentals and Applications, pp. 179–255. Hemisphere, Washington (1985) 29. Launder, B.E., Spalding, D.B.: Lectures in Mathematical Models of Turbulence. Academic, London (1972) 30. Launder, B.E., Reece, G.J., Rodi, W.: Progress in the development of a Reynolds-stress turbulence closure. J. Fluid Mech. 68, 537–566 (1975) 31. Launder, B.E.: On the computation of convective heat transfer in complex turbulent flows. ASME J. Heat Transfer 110, 1112–1128 (1988) 32. Launder, B.E.: Second-moment closure: present... and future. Int. J. Heat Fluid Flow 10, 282–300 (1989) 33. Lumley, J.L.: Computational modeling of turbulent flows. Adv. Appl. Mech. 18, 123–176 (1978) 34. Lumley, J.L., Zeman, O., Siess, J.: The influence of buoyancy on turbulent transport. J. Fluid Mech. 84, 581–597 (1978) 35. Millionshtchikov, M.D.: On the theory of homogeneous isotropic turbulence. C. R. Acad. Sci. S.S.S.R. 32, 615–619 (1941) 36. Moeng, C.-H., Wyngaard, J.C.: Evaluation of turbulent transport and dissipation closures in second-order modeling. J. Atmos. Sci. 46, 2311–2330 (1989) 37. Monin, A.S., Yaglom, A.M.: Statistical Fluid Mechanics: Mechanics of Turbulence. MIT Press, Cambridge, MA (1971) 38. Otić, I., Grötzbach, G.: Direct numerical simulation and statistical analysis of turbulent convection in lead– bismuth. In: International Conference on Supercomputing in Nuclear Applications, Paris Sept. 22–24, 2003, Proceedings on CD-ROM, paper Nr. S10/2, SNA 2003, CEA, 91191 Gif sur Yvette, Cedex, F. (2003) 39. Otić, I., Grötzbach, G., Wörner, M.: Analysis and modelling of the temperature variance equation in turbulent natural convection for low-Prandtl-number fluids. J. Fluid Mech. 525, 237–261 (2005) 40. Otić, I., Grötzbach, G.: Turbulent heat flux and temperature variance dissipation rate in natural convection in lead–bismuth. Nucl. Sci. Eng. 155, 489–496 (2007) 41. Rotta, J.: Statistische Theorie Nichthomogener Turbulenz, 1. Mitteilung. Z. Phys. 129, 547–572 (1951) 42. Sander, J.: Dynamical equations and turbulent closure in geophysics. Contin. Mech. Thermodyn. 10, 1–28 (1998), Springer

154

Flow Turbulence Combust (2007) 79:133–154

43. Schmidt, H., Schumann, U.: Coherent structure of the convective boundary layer derived from largeeddy simulations. J. Fluid Mech. 200, 511–562 (1989) 44. Schumann, U.: Subgrid scale model for finite difference simulations of turbulent flows in plane channels and annuli. J. Comput. Phys. 18, 376–404 (1975) 45. Schumann, U.: The counter gradient heat flux in turbulent stratified flows. Nucl. Eng. Des. 100, 255–262 (1987) 46. Seiter, C.: Numerische Simulation turbulenter Auftriebsströmungen in horizontalen Kanälen. Ph.D. thesis, Univ. of Karlsruhe, Forschungszentrum Karlsruhe FZKA 5505 (1995) 47. Weinstock, J.: A theory of turbulent transport. J. Fluid Mech. 202, 319–338 (1989) 48. Wörner, M.: Direkte Simulation turbulenter Rayleigh–Bénard–Konvektion in flüssigem Natrium. Ph.D. thesis, Univ. of Karlsruhe, Forschungszentrum Karlsruhe KfK 5228 (1994) 49. Wörner, M. Schmidt, M., Grötzbach, G.: Direct numerical simulation of turbulence in an internally heated convective fluid layer and implications for statistical modeling. J. Hydraul. Res. 35, 773–797 (1997) 50. Wörner, M., Grötzbach, G.: Pressure transport in direct numerical simulations of turbulent natural convection in horizontal fluid layers. Int. J. Heat Fluid Flow 19, 150–158 (1998) 51. Zeman, O., Lumley, J.L.: Modeling buoyancy driven mixed layers. J. Atmos. Sci. 33, 1974–1988 (1976)