Eur. Phys. J. E (2013) 36: 120 DOI 10.1140/epje/i2013-13120-2

THE EUROPEAN PHYSICAL JOURNAL E

Regular Article

Phenomenological and statistical analyses of turbulence in forced convection with temperature-dependent viscosity under non-Boussinesq condition S.M. Yahya1,a , S.F. Anwer2 , and S. Sanghi1 1 2

Department of Applied Mechanics IIT Delhi, 110016, India Department of Mechanical Engineering, ZHCET, AMU, Aligarh, 202002, India Received 19 March 2013 and Received in final form 18 July 2013 c EDP Sciences / Societ` Published online: 28 October 2013 –  a Italiana di Fisica / Springer-Verlag 2013 Abstract. In this work, Thermal Large Eddy Simulation (TLES) is performed to study the behavior of weakly compressible Newtonian fluids with anisotropic temperature-dependent viscosity in forced convection turbulent flow. A systematic analysis of variable-viscosity effects, isolated from gravity, with relevance to industrial cooling/heating applications is being carried out. A LES of a planar channel flow with significant heat transfer at a low Mach number was performed to study effects of fluid property variation on the near-wall turbulence structure. In this flow configuration the top wall is maintained at a higher temperature (Thot ) than the bottom wall (Tcold ). The temperature ratio (Rθ = Thot /Tcold ) is fixed at 1.01, 2 and 3 to study the effects of property variations at low Mach number. Results indicate that average and turbulent fields undergo significant changes. Compared with isothermal flow with constant viscosity, we observe that turbulence is enhanced in the cold side of the channel, characterized by locally lower viscosity whereas a decrease of turbulent kinetic energy is found at the hot wall. The turbulent structures near the cold wall are very short and densely populated vortices but near the hot wall there seems to be a long streaky structure or large elongated vortices. Spectral study reveals that turbulence is completely suppressed at the hot side of the channel at a large temperature ratio because no inertial zone is obtained (i.e. index of Kolmogorov scaling law is zero) from the spectra in these region.

1 Introduction Prediction of turbulent flows characterized by large temperature gradients and high heat transfer rates is of great importance in engineering. Heat exchangers, combustion chambers, nuclear reactors, solar receiver and cooling systems in electronic devices are the well-known examples of significant temperature variations typically occuring within the flow. The concerned industrial application of the fundamental study of this article is high-temperature solar receivers. In these systems, thermal and viscous properties of fluids are critical design parameters that vary with temperature. For most conventional fluids, density, specific heat and to a lesser extent thermal conductivity are relatively independent of temperature whereas viscosity is a strong function of it. This temperature-dependent variation alters velocity profiles and changes important flow quantities such as the Nusselt number and friction factor. Among analytical and numerical studies dealing with temperature-dependent viscosity problems, a large majority was devoted to the analysis of laminar heat a

e-mail: yahya [email protected]

convection. Shin et al. [1] developed a finite-volume algorithm to study laminar heat transfer in a rectangular duct. The authors focused on the numerical prediction of friction factor and Nusselt number when constant flux boundary conditions are enforced at the top and bottom walls, and imposed temperature conditions are prescribed at the sidewalls. Results revealed an increase in the Nusselt number with respect to the case of constant viscosity. Pinarbasi et al. [2] studied the influence produced by variations in both thermal conductivity and viscosity of a Newtonian fluid for the case of laminar Poiseuille flow in a two-dimensional channel with constant-temperature boundary conditions. Non-uniform viscosity and thermal conductivity affected the temperature field, while very scarce effects were observed on the velocity field. Zonta and Soldati [3], have performed turbulent forced convection with variable fluid properties with small temperature difference between the walls (40 K and 60 K). Serra et al. (see [4]) also performed LES with two Reynolds number 180 and 395 alongwith three temperature ratios discussing only the effects in the inertial range of energy spectrum when the flow is subjected to a high thermal gradient. However we have only taken fluid viscosity as a variable

Page 2 of 14

Fig. 1. Schematic diagram of channel flow.

and discuss the variable viscosity effects without varying thermal conductivity. Serra did not show near-wall turbulent structures in his paper and confined his study only to energy spectrum and model development, however we have elaborated different effects of variable viscosity on flow physics and bulk parameters. In this work Thermal Large Eddy Simulation (TLES) are performed to investigate the behaviour of turbulence when the fluid viscosity varies significantly with temperature while conductivity and specific heats are assumed constant in the absence of buoyancy. Simulations are based on TLES with momentum and energy equations coupled by the explicit dependence of viscosity through Sutherlands law. With reference to the schematics of fig. 1, we consider a weakly compressible and Newtonian turbulent flow of air in a plane channel with differentially heated walls: the hot wall is kept at temperature Th while the cold wall is kept at temperature Tc . The flow is driven by a constant pressure gradient along the streamwise x-direction. The channel walls are normal to the z-direction and are at constant temperature. The no-slip velocity condition is applied for the fluid velocity at walls, the boundaries of domain normal to x- and y-directions are periodic. We propose to investigate numerically the physical mechanisms responsible for the interaction between the thermal and kinetic fields. The typical channel dimensions considered are 4πδ × 4πδ 3 × 2δ with δ = 1. Different values of the friction Reynolds number, Reτ are considered while the Prandtl number is kept constant (P r = 0.71). The range of Reτ was chosen taking into account that the viscous term exhibits a (Reτ )−1 scaling [3] in our governing equation and viscosity effects vanish at very large Reynolds number. Therefore, we selected four values of Reτ in the low/intermediate turbulent range (150, 180, 245, 460) to explore these effects in detail.

2 Governing equations and numerical modeling In low-speed turbulent channel flow applications, the low Mach-number, variable-density approximation of the Navier-Stokes equations is a good basis for simulation, as it supports large density variations while eliminating acoustic waves. This eliminates the need for extremely

Eur. Phys. J. E (2013) 36: 120

small time steps driven by the acoustics. This means that the arising velocities are much smaller than the speed of sound, so that density variations due to pressure variations can be neglected. Firstly the low Mach-number approximation of the Navier-stokes equations is obtained as the low Mach-number asymptotic limit of the compressible Navier-Stokes equations. Details of the derivation of these equations can be found in [5–8]. We Favre-averaged and filtered equations using an implicit filter to obtain the governing equation for channel flow configuration. The filtering operation can be written in terms of convolution integral as  ¯ G(x − x )f (x )dx . (1) f (x) = D

An f turbulent variable is splitted into an f˜ large component and f  sub grid component. Note that ∼ corresponds to the Favre average operator where Favre average is defined as ρf . (2) f˜ = ρ¯ The Favre-filtered low Mach equation for TLES are given below [9] ∂ρ ∂ ρ¯u ˜j + = 0, ∂t ∂xj

(3)

ρu ˜i u −∂ P¯ ˜j ) ∂(¯ ρu ˜i ) ∂(¯ + = ∂t ∂xj ∂xi ∂(˜ τij − σij ) + + F δi1 , (4) ∂xj      ¯  ∂ (λ) ˜ ∂ T˜ ∂xj ρu ˜j T˜) dP0 ∂(¯ ρT˜) ∂(¯ + + =Γ Cp ∂t ∂xj dt ∂xj −Cp

∂ ρ¯(Hj ) , ∂xj

P¯0 = ρ¯T˜, ∂ P¯0 = 0, ∂xi

(5) (6) (7)

where ui is the i-th component of the velocity vector, P and P0 are the fluctuating kinematic and thermodynamic pressure, respectively, T is temperature, while F represents the constant streamwise presssure gradient that γ , where γ is the specific drives the flow. Here, Γ = (γ−1) heat ratio. All the variables have been normalized using the reference state ρref , uref , T ref , Cpref = Cp (T ref ), μref = μ(T ref ), λref = λ(T ref ). Also, Reynolds number and ref ref ref Prandtl number are defined as Re = ρ uμrefL and P r = μref Cpref λref

, respectively. Further, as our system is closed, P0 may vary in time but the total mass (M0 ) of the system remains constant hence P0 =

M0 V

dV T

.

(8)

Eur. Phys. J. E (2013) 36: 120

Page 3 of 14

The balance equation expressing the fluid viscosity is written as (10) μ(T ) = μref + μv (T ), where μref is the viscosity at reference temperature Tref = (Th + Tc )/2 and μv (T ) is the local deviation from μref . Other thermophysical properties like specific heat Cp and thermal conductivity λ were assumed to be constant. This assumption is made primarily because our study focuses on forced-convection heat transfer in zero-gravity conditions and aims at isolating the effect of variable viscosity on the flow field. Taking into account this filtering procedure, non-linear addition terms appear which represent the stresses and the fluxes, here σij and Hj are turbulent shear stress and tubulent heat flux, respectively, which are defined as ˜i u˜j ) = −2νsgs S˜ij , σij = (u

i uj − u   ∂ T˜ ˜j T˜ = κsgs . Hj = u

jT − u ∂xj

4 variation of viscosity of air 3.5 -

Dynamic viscosity x 10 5

Here molecular viscoscity is calculated with the help of Sutherland law as follows:   32   Tref + 110.4 T μ(T ) . (9) = μref Tref T + 110.4

3 2.5 2 1.5 1 0.5 0 0

100

200

300

400

500

600

700

800

900

T

Fig. 2. Variation of viscosity of air with temperature.

use a subgrid-scale diffusivity κsgs given by eq. (12). The subgrid-scale thermal diffusivity can then be expressed as a function of subgrid-scale kinematic viscosity: κsgs =

(11)

νsgs , P rsgs

(13)

Subgrid-scale models must be used to represent these terms and therefore allow the numerical resolution.

where P rsgs is the subgrid-scale Prandtl number. A constant subgrid-scale Prandtl number (P rsgs = 0.71) is used for low-temperature gradient while dynamic Smagoransky formulation is used to evaluate P rsgs at high-temperature gradient as suggested by Nicoud [16].

2.1 Numerical modelling

3 Results and discussion

The employed scheme is an extension to the pressure correction method proposed for variable density flows under low Mach number approximation [9,10]. The algorithm is based on a predictor-corrector time integration scheme that employs a projection method for the momentum equation. A constant-coefficient Poisson equation is solved for the pressure following both the predictor and corrector steps to satisfy the continuity equation at each time step. The proposed algorithm has second-order centrally differenced convective fluxes while diffusive fluxes are viscous implicit fourth order accurate. Spatial discretization is performed on a collocated grid system with Rhie and Chow interpolation [11] to avoid pressure oscillation. Periodicity in x and y is assumed for both velocity and temperature, while no-slip and constant-temperature conditions are imposed at the walls. Once the velocity field is known, the temperature field is computed as a solution of the energy equation and finally viscosity is updated using Sutherlands law because it best matches experimental data of Weast [12] within the range of temperatures examined. For modeling sub-grid stresses, we have used the Wall Adaptive Layer Equation (WALE) model suggested by Nicoud and Ducros [13]. In the WALE model the subgrid scale viscosity accounts for the effect of rotation rate and strain rate of the smallest fluctuations. This model has a correct wall behavior for subgrid stresses [14,15] near walls νsgs ∝ z 3 . For subgrid-scale heat flux Hj we

To the best of the authors knowledge, thermal LES of turbulent channel flow with air as a fluid aimed at the evaluation of variable viscosity effects isolated from gravity at large temperature gradient are not available in the literature. Lessani [17] studied the effect of large temperature gradient on the physics of channel taking air as a fluid but restrict himself to a single Reynolds number Reτ = 180. Serra et al. [18] performed simulation at Reynolds numbers Reτ = 180 and Reτ = 395 to analyse the combined

(12)

Table 1. Simulation parameters. Cases

C-1.1 C-1.2 C-1.3 C-2.1 C-2.2 C-2.3 C-3.1 C-3.2 C-3.3 C-4.1 C-4.3

Reynolds number Reτ 150 150 150 180 180 180 245 245 245 460 460

Mesh spacing in wall unit (Δx+ , Δy + , Δz + ) (19.6, 6.5, 0.15) (23.5, 7.8, 0.18) (32, 8, 0.24) (60.1, 20, 0.46) -

Temperature ratio Rθ = Th /Tc 1.01 2 3 1.01 2 3 1.01 2 3 1.01 3

Page 4 of 14

Eur. Phys. J. E (2013) 36: 120 20

18 Kim et al(1987) Rθ=1.01

15

16 14

2.5*log(x)+ 5.5 /Uτ

/Uτ

12 10

10 8 6

5

4

Nicoud(1998) C-2.2

2 0

0 0.1

1

10

100

0

50

100

150

z+

(a)

250

300

350

(b)

3

4 Kim et al(1987) urms C-2.1 urms Kim et al(1987) vrms C-2.1 vrms Kim et al(1987) wrms C-2.1 wrms

2.5 2 1.5

/Uτ,/Uτ,/Uτ

/Uτ,/Uτ,/Uτ

200 z+

1 0.5 0

Nicoud(1998) urms C-2.2 urms Nicoud(1998) vrms C-2.2 vrms Nicoud(1998) wrms C-2.2 wrms

3.5 3 2.5 2 1.5 1 0.5 0

0

50

100

150

200

0

50

100

z+

150

200

z+

(c)

(d)

1

0.8

0.6 /U2τ

(-Tcold)/(Thot-Tcold)

0.7 0.8

0.6

0.4

0.5 0.4 0.3 0.2

0.2

Nicoud(1998)

0.1

C-2.2

0

0

Kim et al(1987) C-2.1 0

0

0.5

1 z

1.5

20

40

60

80

2

(e)

100

120

140

160

180

z+

(e)

Fig. 3. Mean profile and velocity correlations at Reτ = 180 for validation. (a) Mean velocity profile in log scale at Rθ = 1.01. (b) Mean fluid streamwise velocity at Rθ = 2. (c) RMS of velocity fluctuation at Rθ = 1.01. (d) RMS of velocity fluctuation at Rθ = 2. (e) Mean temperature profile at Rθ = 2. (f) Velocity-velocity correlation at Rθ = 1.01.

effect of thermal conductivity and viscosity of air on the flow physics. While recently Zonta and Soldati [3] took a broad spectrum of Reynolds numbers but at a small temperature gradient with water as a fluid aimed at the evaluation of variable viscosity effects isolated from gravity. In our simulations a broad spectrum of Reynolds numbers is used (see sect. 1) with a large temperature ratio Rθ ; viscosity depends on temperature, which is in turn a

function of the space variable. Local changes in viscosity may produce large fluctuations that are expected to influence significantly the turbulent flow, especially in the near-wall region. A graphical representation of eq. (9) is given in fig. 2, which shows that the ratio of viscosity values between the walls is 1.5 if ΔT = 300 K and 2.1 if ΔT = 600 K. Note that eq. (9) generates an asymmetric velocity profile associated with larger viscosity variations

Eur. Phys. J. E (2013) 36: 120

Page 5 of 14

20

25 C-1.1 C-1.2 C-1.3

18 16 14

20 15 U/Uτ

U/Uτ

12 10 8

C-1.1 cold C-1.1 hot C-1.2 cold C-1.2 hot C-1.3 cold C-1.3 hot

10

6 4

5

2 0

0 0

50

100

150

200

+

z

(a)

0.1

1

10

100

1000

z+

(b)

Fig. 4. (a) Mean fluid streamwise velocity at different temperature ratios for Reτ = 150. (b) Mean velocity prole in log scale at Reτ = 150.

at higher temperatures and vice versa. In the following section we try to characterize this influence from a statistical viewpoint. A summary of the relevant simulation parameters is given in table 1. We performed an extensive code validation against DNS datasets available for the case of uniform and constant fluid viscosity [19] and nonisothermal case with variable viscosity [20] at Reτ = 180. All simulations are carried out at a finer 96×96×96 mesh. Our simulations have been validated by comparisons with DNS in the isothermal and non-isothermal cases. We compared an isothermal LES with the DNS of Kim et al. [19] and a strongly non-isothermal thermal large eddy simulation (Th /Tc = 2) with the DNS carried out by Nicoud [20]. To our knowledge, there do not exist any DNS or LES data for a non-isothermal simulation at Reτ = 245, 460. For all the simulations where DNS data are available, we obtained a fair agreement between our simulations and the reference DNS data for the mean, fluctuation and correlation profiles. For example, in fig. 3, we can see the comparison, for the longitudinal mean velocity profile in log scale for isothermal flow and the longitudinal mean velocity profile for non-isothermal flow, velocity-velocity correlation and the mean temperature profile at Reτ = 180. We have choose a finer mesh in a transverse direction in order to perfectly fit the DNS profiles as suggested by the Chatelain et al. [21] because increasing the number of nodes makes the profile much nearer to the DNS one. However, the corresponding thermal large eddy simulation would not converge within a reasonable time and therefore such a mesh would not allow us to carry out several simulations with a high temperature ratio.

3.1 Influence of variable viscosity on fluid velocity statistics In this section firstly, we analyze the fluid velocity statistics results relative to the simulation at Reτ = 150, since the effects of viscosity we wish to emphasize are qualitatively similar in all cases but trends become more evident at lower Reynolds number. In fig. 4(a) the stream-

wise velocity profile for constant viscosity at Rθ = 1.01 is compared with the profile where there is a large variation in viscosity due to increase in wall temperature ratio. As apparent from figure, the symmetry of the velocity profile is lost due to the occurrence of lower-velocity gradients at the hot wall, favoured by the increase of viscosity with temperature, and higher-velocity gradients at the cold wall. The turbulent profile is somewhat changing towards laminar one at the hot side of the channel and decreases in the near-wall velocity when heat addition increases, while at the cold side the mean velocity profile is almost unaffected but the near-wall velocity decreases. Figure 4(b) shows streamwise velocity in semi-logscale, it is clear from the figure that turbulent region is decreasing on the hot wall as the temperature gradient increases while on the cold wall it increases with temperature ratio. To appreciate further the effect of temperature-dependent viscosity on turbulence intensities, in fig. 5 the root mean square (r.m.s.) of the fluid velocity fluctuations is shown. Again, profiles are not symmetric and deviations from the constant-viscosity situation are observed. Considering the case for constant viscosity i.e. isothermal simulation in fig. 5(a and b) there is an overestimation of the fluctuation level on the hot half of the domain and an underestimation of its level on the cold half. For the longitudinal velocity fluctuation, the profiles related to both simulations differ largely in their peak values at hot wall. Turbulent intensity increases at the cold side while at the hot side it decreases with increase in the variation of temperaturedependent viscosity. For vertical and transverse velocity fluctuation in fig. 5(c and d) the gap between the curves of simulation is more marked. In particular we note that this deviation increases as we go away from the colder half of the channel, and the maximum gap between the two cases is obtained in the vicinity of the hot channel. Figure 6 shows the streamwise and wall normal turbulent heat flux at Reτ = 150. This clearly indicates as that the temperature ratio increases both the flux will increase and is more prominent in the near-wall region of the cold part of the channel, however streamwise flux decreases at the hot wall with increasing thermal gradient. In fig. 7(a) the

Page 6 of 14

Eur. Phys. J. E (2013) 36: 120 3

3 C-1.1 C-1.2 C-1.3

C-1.1 C-1.2 C-1.3

2.5

2

/Uτ

/Uτ

2.5

1.5

2 1.5

1

1

0.5

0.5

0

0 0

20

40

60

80

100

120

140

160

+

180

200

220

z+

z

(a)

(b) 1 C-1.1 C-1.2 C-1.3

1.2

0.8 /Uτ

/Uτ

1

C-1.1 C-1.2 C-1.3

0.8 0.6

0.6 0.4

0.4 0.2

0.2 0

0 0

50

100

150

200

0

50

100

z+

150

200

z+

(c)

(d)

Fig. 5. RMS velocity fluctuations at Reτ = 150 for different cases: (a) urms profile for the cold part of channel; (b) urms profile for the hot part of channel; (c) urms profile for the full channel; (d) urms profile for the full channel.

mean temperature profile shows that the temperature gradient is decreasing near the hot wall for a higher temperature ratio as clearly depicted by the slope while fig. 7(b) shows the Reynolds shear stress, due to the large temperature gradient. The velocity-velocity correlation decreases at the hot side of the wall while there is no remarkable change in the vicinity of the cold wall region. Although we have found that turbulence is promoted in the cold side of the channel, where viscosity is lower and vice versa. The decrease of viscosity near the cold wall, in particular enhances the mean kinetic energy. We have justified these observations with scaling arguments based on local viscosity effects on the inertial and the viscous forces of the flow. To emphasize the Reτ dependence of viscosity variations on turbulence intensities, in fig. 8(a,b,c) we compare the behaviour of RMS velocity at increasing Reτ . Profiles are shown as a function of the wall-normal coordinate z + . With no stratification effects (as occurs in zero-gravity conditions), variable-viscosity effects become less significant as Reτ increases in agreement with the (Reτ )−1 scaling [21]. In fig. 8(d) the kinetic energy profile is shown, the effect of temperature gradient is more pronounced at the hot wall and the kinetic energy decreases with increase in temperature. In table. 2 the Nusselt number N u decreases at the hot side while it changes merely at the cold wall, the coefficient of friction Cf increases as the temperature

ratio increases while the bulk Reynolds number Reb decreases as the temperature increases. It is clear from the table that the friction Reynolds number at the hot side of the wall decreases with increase in heat addition. Hence on the hot side of the wall the profile become almost laminar. In table. 2, the bulk Reynolds number is calculated as follows: Ub h , (14) Reb = νb with Ub =

1 2h



2h

˜ (y)dy. U

(15)

0

and νb is obtained with the bulk temperature Tb  2h 1 Tb = T˜(y)dy. 2h 0

(16)

Reτ h and Reτ c are the turbulent Reynolds number obtained at hot and cold wall, respectively. As can be seen from table. 2 for Reτ = 150, 180 with increasing temperature ratio the turbulent Reynolds number at the hot wall decreases while on the cold wall increases, but for higher Reτ = 245, 460 the flow still remains turbulent at both the walls. The values of friction velocity Uτ at different walls with increasing temperature ratio is given in table. 3. It is clear that at the hot side of the domain value of Uτ h is

Eur. Phys. J. E (2013) 36: 120

Page 7 of 14

8

2 C-1.1 C-1.2 C-1.3

7

-2 /UτTτ

/UτTτ

6

C-1.1 C-1.2 C-1.3

0

5 4 3

-4 -6 -8

2

-10

1

-12

0

-14 0

20

40

60

80

100

120

140

160

z+

(a)

200

220

(b)

1.4

0.8 C-1.1 C-1.2 C-1.3

1.2

C-1.1 C-1.2 C-1.3

0.7 0.6 /UτTτ

1 /UτTτ

180

z+

0.8 0.6 0.4

0.5 0.4 0.3 0.2 0.1

0.2

0

0

-0.1 0

20

40

60

80

100

120

140

160

z+

180

200

220

z+

(c)

(d)

Fig. 6. Turbulent heat flux for different cases: (a) Streamwise turbulent heat flux at Reτ = 150 for cold wall. (b) Streamwise turbulent heat flux at Reτ = 150 for the hot wall. (c) Wall normal turbulent heat flux at Reτ = 150 for the cold wall. (d) Wall normal turbulent heat flux at Reτ = 150 for the hot wall. 1

0.8

0.8

C-2.1 C-2.2 C-2.3

0.6

C-1.2 C-1.3

0.4 /U2τ

(-Tcold)/(Thot-Tcold)

C-1.1

0.6

0.4

0.2 0 -0.2 -0.4

0.2

-0.6 -0.8

0

0 0

0.5

1 z

1.5

2

(a)

50

100

150

200

250

300

350

z+

(b)

Fig. 7. (a) Mean temperature profile at Reτ = 150 along the wall normal direction z. (b) Shear stress profile at Reτ = 180 for different cases varying along z + .

increasing with temperature ratio, however the opposite effect is observed at cold wall with decreasing Uτ c . N u is the Nusselt number at the wall given by

Nu =

δ( ∂T ∂y )w Tw − Tb

,

(17)

where w represent the value at the wall and Tb is the bulk temperature. It is clearly shown in fig. 9 that the rate of change of the mean velocity near the cold wall is higher for the larger temperature ratio between the hot and cold wall while the rate of change of the mean velocity near the hot wall is higher for the lower temperature ratio between the hot and cold wall.

Page 8 of 14

Eur. Phys. J. E (2013) 36: 120 3 C-2.1 urms C-2.3 urms C-2.1 vrms C-2.3 vrms C-2.1 wrms C-2.3 wrms

2.5 2 1.5

/Uτ,/Uτ,/Uτ

/Uτ,/Uτ,/Uτ

3

1 0.5 0

C-3.1 urms C-3.3 urms C-3.1 vrms C-3.3 vrms C-3.1 wrms C-3.3 wrms

2.5 2 1.5 1 0.5 0

0

50

100

150

200

250

300

350

0

50

100 150 200 250 300 350 400 450

z+

z+

(a)

(b) 4 C-4.1 urms C-4.3 urms C-4.1 vrms C-4.3 vrms C-4.1 wrms C-4.3 wrms

3 2.5 2

C-2.1 C-2.2 C-2.3

3.5 3

/Uτ,/Uτ,/Uτ

4 3.5

1.5

2.5 2 1.5

1

1

0.5

0.5

0

0 0

100 200 300 400 500 600 700 800 900

0

50

100

150

z+

200

250

300

350

z+

(c)

(d)

Fig. 8. Root mean square of fluid velocity fluctuations at varying Reynolds number: (a) Reτ = 180. (b) Reτ = 245. (c) Reτ = 460. (d) Turbulent kinetic energy at Reτ = 180.

1.2 C-2.1 C-2.2 C-2.3

1.2

0.8

/Uτ

/Uτ

1.4

C-2.1 C-2.2 C-2.3

1

0.6

1 0.8 0.6

0.4 0.4 0.2

0.2

0 0

1

2

3

4

5

0 350

352

354

z+

(a)

356

358

360

z+

(b)

Fig. 9. Near wall mean velocity profile at Reτ = 180: (a) cold wall, (b) hot wall.

3.2 Influence of variable viscosity on bulk Reynolds number Clearly observing from the previous section, the decrease of turbulence intensity near the hot wall and the increase near the cold wall seem to be a good agreement with bulk Reynolds number. As we see in pipe/channel flow, when viscosity is reduced inertial effects are less damped and

fluctuations are bound to increase. By analogy, one would expect the same reasoning to hold for the two sides of the channel and would intuitively justify the behaviour of turbulence. To discriminate variable-viscosity effects on the boundary layer structure, we can subdivide the flow domain into two distinct sub-channels, one on the hot side and the other on the cold side, separated by the xy plane where the mean shear stress vanishes (located at

Eur. Phys. J. E (2013) 36: 120

Page 9 of 14

Table 2. Effect of variable viscosity on global parameters. Cases

Reb

Reτ h

Reτ c

N uh

N uc

Cf × 103

C-1.1 C-1.2 C-1.3 C-2.1 C-2.2 C-2.3 C-3.1 C-3.2 C-3.3 C-4.1 C-4.3

2120 1671 1556 2787 2120 1789 3239 3000 2415 8760 6142

151 72 51 183 92 63 262 130 86 462 255

154 185 196 187 231 268 277 319 357 468 618

22.23 13.78 10.10 28.44 16.31 7.03 31.19 20.89 11.45 43.75 26.62

23.87 22.65 22.18 29.33 28.88 28.46 32.36 31.06 29.95 44.26 43.33

8.98 9.07 9.97 8.45 8.76 9.02 7.98 8.29 8.64 6.67 7.82

3.3 Influence of variable viscosity on the mean kinetic energy balance

Table 3. Values Uτ at different walls along with Reynolds number. Cases

Reb

Reτ h

Reτ c

Uτ h

Uτ c

C-1.1 C-1.2 C-1.3 C-2.1 C-2.2 C-2.3 C-3.1 C-3.2 C-3.3 C-4.1 C-4.3

2120 1671 1556 2787 2120 1789 3239 3000 2415 8760 6142

151 72 51 183 92 63 262 130 86 462 255

154 185 196 187 231 268 277 319 357 468 618

0.0655 0.0751 0.0955 0.0651 0.0745 0.0799 0.0675 0.0731 0.0788 0.0536 0.0723

0.0647 0.0591 0.0574 0.0648 0.0573 0.0534 0.0667 0.0547 0.0515 0.0533 0.0513

Pk

f1

Tk = − 

  ∂ρku2  1 ∂τi3 ui  ∂ p¯ − u ¯i + , ∂x ∂x Re ∂x3 3   i    f2

f3

f4



 1 ∂u τij i , k = − Re ∂xj   f5

Πk = pressure work,

uhot × hhot × ρhot b b b = 2062. μhot b

(19)

Based on eqs. (18) and (19), we can conclude that, when the pressure gradient that drives the flow is imposed, the mass flow rate in the cold subchannel is higher than the mass flow rate in the hot subchannel. We have examined also the variation of inertial forces, FI , and viscous forces, FV , in the two subchannels. With straightforward meaning of symbols the following scaling holds: (20)

Upon computation of these forces in each subchannel, at Reτ = 150 we obtain F cold FVcold = 0.48 : Ihot = 1.68. hot FV FI

Πk

∂u ¯1 Pk = − ρu1 u3 , ∂x3  

(18)

FI ∼ ρref u2b h2 ∼ μb (ub /h)h2 .

k

Tk

where

ucold × hcold × ρcold b b = b = 2651, cold μb

= Rehot b

In this section, we examine the mean KE balance. We consider the following transport equation of the ensembleaveraged mean kinetic energy of turbulence to analyze the transport processes that take place within the flow Dk = Production +Dissipation+F u1 , (22) +Transport       Dt 

z + = 100). For each sub-channel, we obtain the following bulk Reynolds numbers: Recold b

As expected, in the cold subchannel there is a significant decrease of the viscous forces (−52%) due to the decrease of viscosity. However, there is an even larger increase of the inertial forces (+68%), in agreement with the values bulk Reynolds numbers at both walls. For increasing Re both force ratios should converge to unity and the flow symmetry should be recovered. In our simulations we obtain hcold /hhot ∼ 1.38, 1.33, 1.29 and 1.01 at Reτ = 150 (C-1.3), 180 (C-2.3), 245 (C-3.3) and 460 (C-4.3). This result confirms that the variable-viscosity effects we observe will fade out at higher Reynolds numbers.

(21)

where f 1, f 2, f 3, f 4, f 5, and f 6 correspond to shear production, turbulent transport, mean of the velocity and the pressure gradient, viscous transport, viscous dissipation, and driving pressure gradient work, respectively [22]. The energy for the velocity field is supplied to the fluid by the action of the mean flow through production. This energy can be either dissipated directly by molecular viscosity or converted into turbulent velocity fluctuations. Energy is received entirely by the streamwise component and is transferred to the other components: the mean flow energy derived from the imposed pressure gradient is convected to the wall by turbulence and once in the near-wall region, production increases due to the presence of large mean velocity gradients. We have neglected (f 3) the term which is very small compared to other terms and as expected, the five terms are in balance. In fig. 10 in the immediate vicinity of the walls the viscous transport (f 4) term cancels the viscous dissipation (f 5) term. Slightly away from the walls, the turbulent transport (f 2) term reaches

Page 10 of 14

Eur. Phys. J. E (2013) 36: 120

0.06 Production Viscous Dissipation Viscous transport Turbulent transport Pressure work

0.04

K

0.02 0 -0.02 -0.04 -0.06 0

0.2 0.4 0.6 0.8

1

1.2 1.4 1.6 1.8

2

z

Fig. 10. Mean kinetic energy budget for constant viscosity (C-2.1) case at Reτ = 180.

its maximum and it is balanced by three terms; the production (f 1), the viscous transport (f 4), and the viscous dissipation (f 5). The viscous transport (f 4) term is positive at the walls but it becomes negative slightly away from the walls, and has negligible value at the channel core region. At the channel core region, the balance is between the positive turbulent transport (f 2) term, and the negative viscous dissipation (f 5) term. To analyze possible changes in mean KE redistribution with variable viscosity, we examine fig. 11. Particular emphasis is put on production and dissipation, which are the terms most affected by viscosity variations. The main conclusion that can be drawn from fig. 11 and fig. 12 is that the peak production of mean KE decreases near the hot wall and increases near the cold wall. This is in agreement with the finding that turbulence is damped at the hot wall and promoted at the cold wall (see also fig. 5). In the hot side of the channel, dissipation decreases outside the viscous sublayer but increases right at the wall where it is maximum. This is a direct consequence of the high-velocity fluctuation gradients in this region. In the cold side of the channel, the behaviour is opposite: dissipation increases outside the viscous sublayer and decreases sharply at the wall. As for the other budget terms, only viscous transport shows non-negligible changes, increases at the hot wall whereas a decrease is observed at the cold wall.

3.4 Influence of variable viscosity on instantaneous turbulence structures As mentioned by Lesieur et al. [23], there are different types of methodology available for identifying vortices in a turbulent flow. First route is through pressure, the cores of the coherent structure should be having lowest pressures. The second method is to use high vorticity (H = |ω|), the positive Q and the last is the λ2 criterion. So the isosurfaces of each of these identify the vortices. The choice of the threshold is justified by what visually gives the

best vortices. In this work we use isosurfaces of the positive Q to visualize these vortices. Here, Q is defined as → ω 2 − 2Sij Sij ). Figure 13(a,b,c) shows an isometQ = 14 (− ric view of isosurfaces of the channel with Q = 0.2. It can be seen that, as viscosity varies due to increase in temperature ratio, the vortical activity decreases on the warm side of the channel while it increases at the cold side of the channel. This is evident from the absence of any vortices in fig. 13(c). Fluid velocity streaks are shown in fig. 14 for the cold wall (fig. 14a) and for the hot wall (fig. 14b) for C-2.1 and fig. 14(c,d) for C-2.3 for cold and hot wall, respectively. A greyscale mapping is used for the instantaneous streamwise fluid velocity fluctuations ux in the wall-parallel xy plane located at a distance z + = 12 from each wall. Dark-grey areas indicate regions of negative ux and mark the location of low-speed streaks, lightgrey areas indicate regions of positive ux and mark the location of high-speed streaks. Low-speed streaks appear to be stronger, more pronounced (as shown by the steeper greyscale gradients) and arguably more persistent at the hot wall. This is in an agreement with the behaviour of the strain rate discussed by [3], near the hot wall the strain is increased slightly yet enough to favour streak formation, with an opposite behaviour near the cold wall. The persistence of the streaks, defined by their extent in the wall-normal direction, is connected directly to the local shear rate and to the turbulent dissipation rate.

3.5 Effect of variable viscosity on turbulent kinetic energy spectra In this section, we describe the phenomenon of strong coupling between turbulence and temperature kinetics which disturbed the Kolmogorov −5/3 scaling law when the flow is subjected to a strong temperature gradient, due to which one more time scale arises along with the non-linear time scale called temperature gradient time scale. Here we modify the spectrum model [4] which will take into account the influence of temperature gradient in estimating the time scale for triple decorrelation. The model shows a good agreement with the results of TLES (Thermal Large Eddy Simulation) performed at different friction Reynolds numbers and temperature ratios. The spectrum has been obtained at the middle of the domain in the z + direction. This analysis provides a better way to understand the underlying physics of anisothermal turbulent channel flow. ⎛ E(k) = Ck ⎝

3 Urms h 

⎞2/3 ⎠

Rθ 1+ f1

2/3  1 k −5/3 1+ . f kh

(23)

In the asymptotic limit of the isothermal case, this spectrum is the Kolmogorov spectrum with the −5/3 slope. In the asymptotic limit of the very strong temperature gradient, the slope of the spectrum is equal to −7/3 [4]. For intermediate strength of the temperature gradient, the spectrum varies smoothly between these two limiting forms according to the decrease of the controlling function

Eur. Phys. J. E (2013) 36: 120

0.08

0.08

Production Viscous Dissipation Viscous transport Turbulent transport Pressure work

0.06 0.04 0.02

Production Viscous Dissipation Viscous transport Turbulent transport Pressure work

0.06 0.04 0.02

0

0

K

K

Page 11 of 14

-0.02

-0.02

-0.04

-0.04

-0.06

-0.06

-0.08

-0.08 0

0.1

0.2

0.3

0.4

0.5

1.5

1.6

1.7

z

1.8

1.9

2

z

(a)

(b)

Fig. 11. Mean kinetic energy budget for constant viscosity (C-2.2) case at Reτ = 180: (a) cold wall, (b) hot wall.

0.15

0.15 Production Viscous Dissipation Viscous transport Turbulent transport Pressure work

0.1

0.05

0

0

K

K

0.05

Production Viscous Dissipation Viscous transport Turbulent transport Pressure work

0.1

-0.05

-0.05

-0.1

-0.1

-0.15

-0.15 0

0.1

0.2

0.3

0.4

0.5

1.5

1.6

1.7

z

1.8

1.9

2

z

(a)

(b)

Fig. 12. Mean kinetic energy budget for constant viscosity (C-2.3) case at Reτ = 180: (a) cold wall, (b) hot wall.

f . Figure 15, represents the comparison between turbulent fluctuation profiles obtained by simulation and the spectra calculated using the model (eq. (23)). In fig. 15(a-f) we have ploted model data using square, Kolmogorov −5/3 law using circles and the simulation result using triangles. We can see that the model obtained very good matches in the inertial range of the spectrum with the simulation result. The developed model is fully devoted to consider the effects of temperature gradient on the turbulence spectra. Here two cases follow, isothermal and with large temperature gradient and intermediate cases. Figure 15(a-c) corresponds to the cases for Reτ = 245 while fig. 15(d-f) for Reτ = 460. In fig. 15(b and c) we have seen that for cold side of the wall the spectrum follows the −7/3 slope at Rθ = 2 and Rθ = 3, it shows a good match with simulation results as well as the desired slope. Figure 15(d,e) represent the cold side and hot side of the spectrum, respectively, in which the cold side follows −5/3 law while the hot side deviates from it and decays more rapidly with power −7/3. At a large temperature gradient the spectra at the cold wall follows the −7/3 law while the hot side is

Table 4. Characteristics of the non-isothermal spectra. Reτ = 245 Rθ = 2

Rθ = 3

Rθ = 2

Reτ = 460 Rθ = 3

slope k−σ

slope k−σ

slope k−σ

slope k−σ

Sh

N/A

N/A

7/3

N/A

Sc

7/3

7/3

5/3

7/3

no more turbulent and hence there is no inertial zone from the spectra obtained in these regions. If the value bulk Reynolds number is higher than the critical value (around 2000), then the slope of the spectrum at the hot side is −5/3 while at cold side is −7/3. Sh and Sc are the hot and cold side, respectively. Clearly from table 4 the hot side of the simulations 245-2, 245-3 and 460-3 is no more turbulent, that is why no inertial zone is obtained from the spectra in these regions. While on the cold side of the simulation the slope obtained for the spectrum is k −7/3 .

Page 12 of 14

Eur. Phys. J. E (2013) 36: 120

Fig. 13. Isosurfaces of Q = 0.2 at Reτ = 180 for three cases: (a) C-2.1, (b) C-2.2, (c) C-2.3.

Fig. 14. Visualization of velocity streaks in the near-wall region for the channel flow at Reτ = 180 (C-2.1): (a) cold wall, (b) hot wall and (C-2.3), (c) cold wall, (d) hot wall at Z + = 12.

Eur. Phys. J. E (2013) 36: 120

Page 13 of 14

1

1 spectrum model k(-7/3) simulation 245-2 cold side

0.1

E(k)

E(k)

0.1 0.01

0.01 0.001

spectrum model k(-5/3) simulation 245-iso

0.0001

0.001 0.1

1 kh

10

0.1

1

1 kh

10

1 spectrum model k(-7/3) simulation 245-3 cold side

E(k)

0.1

E(k)

0.1

0.01

0.01 spectrum model k(-5/3) simulation 460-2 cold side

0.001

0.001 0.1

1 kh

10

0.1

0.1

0.1

10

E(k)

1

E(k)

1

1 kh

0.01

0.01 spectrum model

spectrum model

k(-7/3) simulation 460-2 hot side

k(-7/3) simulation 460-3 cold side

0.001

0.001 0.1

1 kh

10

0.1

1 kh

10

Fig. 15. Turbulent fluctuation spectra at (a) Reτ = 245 and Rθ = 1. (b) Reτ = 245 and Rθ = 2 cold side. (c) Reτ = 245 and Rθ = 3 cold side. (d) Reτ = 460 and Rθ = 2 cold side. (e) Reτ = 460 and Rθ = 2 hot side. (f) Reτ = 460 and Rθ = 3 cold side.

4 Conclusion The convective heat transfer in turbulent shear flows is characterized by complex macroscopic phenomena that are rich in physics yet difficult to model and/or to reproduce. Complexity becomes particularly clear in situations where the thermophysical properties of the fluid depend strongly on temperature. This is the case of airbased convective heat transfer equipments, in which the

occurrence of large temperature gradients produces significant variations of density, dynamic viscosity and thermal conductivity. In this paper, we focus on the effects produced by temperature-dependent fluid viscosity, isolated from gravity, on the flow field in turbulent forced convection subjected to large temperature gradient. We performed numerical experiments based on TLES simulations of convective heat transfer in turbulent channel flows with temperature-dependent viscosity. Simulations

Page 14 of 14

were run at Reynolds numbers Reτ = 150, 180, 245 and 460 and Prandtl number P r = 0.71. Having neglected gravity, we have been able to emphasize the non-trivial behaviour due to variable viscosity and the way it is modulated by a change in Reτ . The effects we observe in zerogravity conditions fade out as Re increases while due to the large thermal gradient variable-viscosity effects mingle with high temperature effects and it is sometime difficult to understand both effects simultaneously using LES under non-Boussinesq condition. We find that turbulence is promoted in the cold side of the channel, where viscosity is lower while turbulence is damped in the hot side of the channel, where viscosity is higher. The decrease of viscosity near the cold wall, in particular, enhances the mean KE, while reducing the mean velocity gradient. We have justified these observations with scaling arguments based on local viscosity effects on the inertial and viscous forces of the flow. During investigation of the mean KE budget it is found that moving from the central region of the channel toward the hot wall, both production and dissipation decrease but in the near-wall region production again starts to increase. In the cold side of the channel, the peak of the production decreases in comparison to the constant viscosity case and the dissipation increases reaching a maximum value of a few viscous units from the wall. Finally we have performed an analysis of energy spectra at different temperature ratios and friction Reynolds number and found that the variable-viscosity results due to large temperature gradient changes the spectrum model in the inertial zone and one more time scale arises along with non-linear time scale called temperature gradient time scale. For the isothermal case, this spectrum is the Kolmogorov spectrum with −5/3 slope, while for very strong temperature gradient, the slope of the spectrum is equal to −7/3. Our analysis also reveals that turbulence is completely suppressed at the hot side of the wall at large temperature ratio because no inertial zone is obtained from the spectra in these regions. CFD Lab Applied Mechanics Department (IIT Delhi, India), and High Performance Computing Lab ZHCET (AMU, Aligarh, India) initiative are gratefully acknowledged for generous support and allowance of computer resources.

Eur. Phys. J. E (2013) 36: 120

References 1. S.Y. Shin, Y.I. Cho, W.K. Gringrich, W. Shyy, Intl. J. Heat Mass Transfer 36, 4365 (1993). 2. A. Pinarbasi, C. Ozalp, S. Duman, Phys. Fluids 17, 038109 (2005). 3. F. Zonta, C. Marchioli, A. Soldati, J. Fluid Mech. 697, 150 (2012). 4. S. Serra, A. Toutant, F. Bataille, Y. Zhou, Phys. Lett. A 376, 3177 (2012). 5. A. Majda, J.A. Sethian, Combust. Sci. Technol. 42, 185 (1985). 6. R.G. Rehm, H.R. Baum, N. B. S. J. Res. Tech. 83, 297 (1978). 7. B. Mller, J. Eng. Math. 34, 97 (1998). 8. S. Paolucci, Technical report, Sandia National Laboratories (1982). 9. S.M. Yahya, S.F. Anwer, S. Sanghi, World J. Mech. 5, 253 (2012). 10. S.F. Anwer, H. Naushad Khan, S. Sanghi, A. Ahmad, S.M. Yahya, AIP Conf. Proc. 1440, 683 (2012). 11. C.M. Rhie, W.L. Chow, AIAA j 21, 11 (1983). 12. R.C. Weast, CRC Handbook of Chemistry and Physics (CRC Press, 1988). 13. F. Nicoud, F. Ducros, Flow, Turbulence Combust. 62, 183 (1999). 14. M. Germano, U. Piomelli, P. Moin, W.H. Cabot, Phys. Fluids A. 3, 1760 (1991). 15. P. Moin, K. Squires, W. Cabot, S. Lee, Phys. Fluids. 3, 2746 (1991). 16. F. Nicoud, T. Poinsot, 1st International Symposium on Turbulence and Shear Flow Phenomena, Santa Barbara, Sept. 12-15, 1999. 17. B. Lessani, M.V. Papalexandris, J. Comput. Phys. 212, 218 (2006). 18. S. Serra, A. Toutant, F. Bataille, Heat Transfer Eng. 33, 505 (2011). 19. J. Kim, P. Moin, R. Moser, J. Fluid Mech. 177, 133 (1987). 20. F.C. Nicoud, Annual Research Briefs, Center for Turbulent Research (1998) p. 289. 21. A. Chatelain, F. Ducros, O. Metais, Int. J. Num. Methods Fluids 44, 1017 (2007). 22. B. Lessani, A. Zainali, J. Turbulence 10, (2009). 23. M. Lesieur, O. Metais, P. Comte, Large-Eddy Simulations of Turbulence (Cambridge University Press, Cambridge, 2005) p. 219.