- Email: [email protected]
Effects of Rayleigh-Bénard convection on spectra of viscoplastic fluids
International Journal of Heat and Mass Transfer 147 (2020) 118947
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Effects of Rayleigh-Bénard convection on spectra of viscoplastic fluids Sahin Yigit a, Josef Hasslberger a, Nilanjan Chakraborty b,⇑, Markus Klein a a b
Instute of Applied Mathematics and Scientific Computing, Bundeswehr University Munich, Germany School of Engineering, Newcastle University, Newcastle upon Tyne, United Kingdom
a r t i c l e
i n f o
Article history: Received 19 June 2019 Received in revised form 21 October 2019 Accepted 22 October 2019
Keywords: Direct Numerical Simulation Rayleigh-Bénard convection Non-Newtonian Yield stress fluid High nominal Rayleigh number
a b s t r a c t Direct Numerical Simulations (DNS) have been conducted to analyse velocity and temperature spectra in unsteady three-dimensional Rayleigh-Bénard convection (differentially heated horizontal walls, i.e. heated from below and cooled from above whereas the other walls are adiabatic) at high nominal Rayleigh number of viscoplastic fluids obeying a Bingham model in cubic domains and the corresponding Newtonian cases with same nominal Rayleigh and Prandtl numbers for reference purposes. The simulations have been carried out for nominal Rayleigh numbers of 107 and 108 for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution). The effects of yield stress on velocity and temperature spectra and heat transfer characteristics of the flow fields have been investigated in terms of isotherms, path lines, apparently unyielded regions and mean Nusselt number, along with the power spectral densities of temperature and velocity fluctuations. It has been found that thermal convection strengthens (weakens) with increasing nominal Rayleigh number Ra (Bingham number Bn). Consistent with this observation, the range of flow structures becomes larger with higher (smaller) values of Ra (Bn). Additionally, it has been observed that the decay of the power spectral densities of temperature fluctuations in the viscous-convective range is close to the 1 inverse of frequency f (i.e. it scales as f ) for Pr > 1 in both Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases with moderate values of Bn. The fluid motion (including its unsteady features) weakens and the mean Nusselt number drops progressively with an increase in Bingham number. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Non-Newtonian fluids are immensely important in many industries such as polymer and chemical processing, biochemical synthesis. Yield stress fluids belong to a special type of non-Newtonian fluids, which flow once a critical stress (i.e. yield stress) is surpassed but act like a rigid solid below this threshold stress. Rayleigh-Bénard convection (i.e. the buoyancy-driven flow of a fluid heated from below and cooled from above) of yield stress fluids in enclosed spaces has a wide range of engineering applications, such as preservation of canned food, chemical processing, solar and nuclear power systems – to name only a few. Aqueous solutions of Carbopol are often used as model yield stress fluids in laboratory experiments [1,2]. Some magneto and electro-rheological fluids exhibit yield stress behaviour, where it is possible to modify the yield stress by applying electrical or magnetic fields [3,4]. This can be useful in reducing heat losses or designing new adaptive thermal systems to control convective thermal transfer in enclosed spaces. ⇑ Corresponding author. E-mail address: [email protected] (N. Chakraborty). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118947 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
Several recent studies [1,2,5–14] have analysed laminar Rayleigh-Bénard convection of yield stress fluids obeying a Bingham model in two-dimensional rectangular enclosures for small and moderate values of nominal Rayleigh number. These analyses revealed that the flow progressively weakens with increasing Bingham number (i.e. non-dimensional yield stress) due to additional flow resistance arising from the yield stress. This leads to a reduction in convective transfer rate, and thermal transport takes place principally due to conduction for large values of Bingham number. This behaviour arises because the fluid flow practically stops for large values of Bingham number and it does not influence the thermal transport under such a condition. A few studies [1,2,5–8] focused on critical conditions for the onset of Rayleigh-Bénard convection of Bingham fluids, whereas others [9–14] concentrated on heat transfer characteristics under conditions which are far beyond the critical value for the onset of fluid motion. It is worth noting that all the aforementioned studies in the existing literature [1,2,5–14] are restricted to laminar flow (i.e. nominal Rayleigh number from 103 to 105). However, Rayleigh-Bénard convection of yield stress fluids at high nominal values of Rayleigh number for which a steady laminar solution may not exist, has not yet been analysed in the existing literature
2
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
and consequently, little is known about the energy cascade in this kind of flows. This gap has been addressed here by conducting Direct Numerical Simulations (DNS) of unsteady threedimensional Rayleigh-Bénard convection of yield stress fluids obeying a Bingham model in cubic domains at high values of nominal Rayleigh number for which the Newtonian fluids (and also Bingham fluids of small values of Bingham number) are expected to exhibit turbulent natural convection. The simulations have been 7
8
carried out for nominal Rayleigh numbers of 10 and 10 for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution [2]). Accordingly, the main objectives of this study are: i) To demonstrate the effects of Rayleigh and Bingham numbers on Rayleigh-Bénard convection of yield stress fluids in cubic domains for high nominal Rayleigh numbers. ii) To provide physical explanations for the observed Rayleigh and Bingham number dependence of mean Nusselt number in the aforementioned configuration. iii) To analyse the effects of Rayleigh, Prandtl and Bingham numbers on velocity and temperature spectra at the centre of the flow field. The structure of the paper is as follows: Section 2 will introduce the governing equations and the physical background. The Numerical scheme is detailed in Section 3, followed by the discussion of the results in Section 4. The concluding remarks are presented at the end of the paper.
2.1. Modelling of the yield stress fluid The bi-viscosity regularisation proposed by O’Donovan and Tanner [15] has been used to mimic the stress-shear rate characteristics for Bingham fluids:
s ¼ sy c_ =c_ þ lc_ for c_ > sy =lyield
beyond
which
the
linear
stress-strain
rate
relation
(i.e.
s ¼ sy ðc_ =c_ Þ þ lc_ ) according to the Bingham model is invoked. Therefore, the precise shapes and sizes of Apparently Unyielded Regions (AURs) do not affect the mean Nusselt number in the present study. O’Donovan and Tanner [15] indicated that a value of lyield equal to 1000 l mimics the true Bingham model in a satisfactory manner but here lyield =l ¼ 104 is kept fixed to ensure higher fidelity of the simulations. This can also be confirmed by the mean Nusselt number comparison (i.e. Nu) for Ra ¼ 107 , Pr ¼ 320 and Bn ¼ 0:43, as Nu is obtained to be 11.093 and 11.095 for the values of
lyield =l ¼ 103 and 104 respectively.
Additionally, several studies analysed natural convection of Bingham fluids in a square cavity with differentially heated sidewalls obeying the Bingham model without using any regularisation [17,18] with numerical simulations. They compared the results obtained using the bi-viscosity regularisation [5] with the ones without regularisation [17,18]. It has been reported by Refs. [17,18] that the mean Nusselt number values obtained for low Ra (i.e. Ra ¼ 104 ) remain in good agreement, but small differences have been reported for high values of Ra (i.e. Ra ¼ 105 ) between the Bingham model without regularisation and bi-viscosity regularisation (i.e. a maximum of 5% difference in the mean Nusselt number has been reported for Bn ¼ 3, Ra ¼ 105 at Pr ¼ 1 between Refs. [18] and [19] in the Table 2 of Ref. [18]). However, this difference remains within the typical computational uncertainty and may also arise because of the differences in the mesh and numerical schemes.
2. Mathematical background
s ¼ lyield c_ for c_ sy =lyield
Above the critical shear rate (where jsj > sy ), qualitative and quantitative distributions of flow patterns remain independent of the value of lyield =l. This is because, as indicated by Eq. (2), the value of lyield determines the critical shear rate (i.e. c_ c ¼ sy =lyield )
ð1Þ ð2Þ
where c_ ij ¼ ð@ui =@xj þ @uj =@xi Þ are the components of the rate of strain tensor c_ , s is the stress tensor,sy is the yield stress, l is the plastic viscosity and lyield is the yield viscosity. Here, s and c_ are evaluated based on the second invariants of the stress and the rate of strain tensors in a pure shear flow respectively 0 " #1=2 " #1=2 ! @i:e: s ¼ 1=2 s : s _ _ _ , c ¼ 1=2c : c .
2.2. Computational domain The simulation configuration is schematically shown in Fig. 1, which demonstrates that the differentially heated horizontal walls are subjected to constant wall temperature boundary conditions. The bottom wall is taken to be at higher temperature than the top wall (i.e. T bottom > T top ) and all the other walls are considered to be adiabatic (i.e. the temperature gradient in the wall normal direction vanishes at the wall). No-slip and impermeability conditions are specified for all walls. The enclosure is taken to be cubic (i.e. H ¼ L ¼ W).
O’Donovan and Tanner [15] indicated that a value of lyield equal to 1000l mimics the true Bingham model in a satisfactory manner. Eqs. (1) and (2) transform ‘‘unyielded” region to a zone of high viscosity here in order to avoid singularities arising from a pure conduction solution (i.e. where c_ ¼ 0). As noticeable from Eqs. (1) and (2) the value of lyield represents the sensitivity of critical shear rate (i.e. c_ sy =lyield or c_ > sy =lyield ) here in order to mimic satisfactorily the Bingham model. This value affects only the size of the unyielded regions (regions where jsj sy [16]). The unyielded regions in Bingham fluid flows can be indicated by the locations of the Apparently Unyielded Regions (AURs). These regions are dependent on the prediction of jsj sy which is a function of the choice of lyield =l. Therefore, the grey regions should not be treated as exact unyielded zones and for this very reason Mitsoulis and Zisis [16] termed them as the Apparently Unyielded Regions.
Fig. 1. Schematic diagram of the simulation domain.
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
2.3. Governing equations The transient flow of an incompressible Bingham fluid is considered. The conservation equations for mass, momentum and energy for incompressible fluids under transient conditions take the following form:
@ui ¼0 @xi
@ui @ui 1 @P @ sij ¼ þ þ g b T T ref di2 þ uj @t @xj q @xi @xj
@T @T @2T ¼a þ uj @t @xj @xj @xj
ð3Þ ð4Þ
ð5Þ
where ui is the ith component of the velocity vector, T and P represent the temperature and pressure respectively, g; b are the acceleration due to gravity, the thermal expansion coefficient and a is the thermal diffusivity. The last term on the right-hand side of Eq. (4) originates due to Boussinesq’s approximation and the temperature difference between horizontal walls is considered to be small enough so that this approximation remains valid. Also, in Eq. (4), the Kronecker delta di2 indicates that buoyancy forces affect the flow only in the vertical direction (i.e. y direction). The reference temperature T ref is taken to be the cold wall temperature (i.e. T ref ¼ T C ). Eqs. (3)–(5) are valid for all fluid flow problems as long as continuum hypothesis remains valid, which is expected to be the case for most natural convection problems involving non-Newtonian fluids in engineering applications. The relation between s and c_ determines whether these equations are valid for Newtonian fluids or non-Newtonian fluids. The mass conservation and energy conservation equations (i.e. Eqs. (3) and (5)) do not change whether the fluid is Newtonian or non-Newtonian as long as the assumption of incompressible flow is valid, which is the case here. The same governing equations were used for the analysis of natural convection of Bingham fluid flows in several previous analyses (e.g. [5,7–14] and references therein) and the same approach has been adopted here. The viscous stress components sij are expressed with the help of Eqs. (1) and (2) (note that sy vanishes for Newtonian fluids). It is worth noting that in Eq. (5), the thermal diffusivity of the Bingham fluid cases (i.e. a) is kept independent of the shear rate. This can be also confirmed by the experimental studies [1,2] on Rayleigh-Bénard convection of yield stress fluids (e.g. Carbopol gel). It should be also noted that the shear rate of the Bingham fluids decreases immensely with increasing Bn for given value of shear stress. Therefore, change in the thermal diffusivity with the shear rate is expected to be negligible for the yield stress fluids. There have been several previous analyses [9–14] where the thermal diffusivity is considered to be independent of shear rate and the same approach has been adopted in this study. 2.4. The non-dimensional numbers According to Buckingham’s pi theorem, it is possible to show that the Nusselt number for natural convection of Bingham fluids in cubic enclosure can be expressed as: Nu ¼ f ðRa; Pr; Bn), where the nominal Rayleigh, Prandtl and Bingham numbers for Bingham fluids can be defined in the following manner:
Ra ¼ qgbðT H T C ÞH3 =ðlaÞ
ð6iÞ
Pr ¼ l=ðqaÞ
ð6iiÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Bn ¼ sy H=ðl gbðT H T C ÞHÞ
ð6iiiÞ
3
Here, the Bingham number represents the ratio of yield stress to viscous stress (i.e. non-dimensional yield stress). The nondimensional local heat transfer coefficient h is defined as:
hhw ¼ kð@T=@yÞy¼0 =ðT y¼0 T y¼H Þ
ð7iÞ
hcw ¼ kð@T=@yÞy¼H =ðT y¼H T y¼0 Þ
ð7iiÞ
where hw and cw represent the values at the hot wall and cold wall respectively, and k is the thermal conductivity. The mean heat transfer coefficient h and mean Nusselt number Nu at the walls are evaluated as:
h¼
1 LW
ZZ
hdxdz; Nu ¼ h H=k
ð8Þ
It is worth noting that hhw and hcw in Eq. (7) can be used in Eq. (8) in order to calculate mean Nusselt number at hot and cold walls (i.e. Nu hw and Nu cw ). 3. Numerical implementations 3.1. Numerical methodology To solve the non-linear set of governing equations in a finitevolume framework, the open-source CFD package OpenFOAM has been utilised. The pressure-velocity coupling has been addressed by the use of the Merged PISO-SIMPLE (PIMPLE) algorithm, which has been optimised for unsteady problems. Convective and diffusive fluxes are evaluated by second-order centered difference schemes. Temporal advancement has been achieved by the second-order Crank-Nicolson scheme with constant timestepping. It has been ensured that the Courant number is always sufficiently below unity so that the underlying physics is captured with sufficient temporal resolution. 3.2. Resolving the fine-scale structure It is vital to resolve all the relevant scales for conducting Direct Numerical Simulations (DNS) since the gradients of the turbulent fields cannot be accurately represented in the absence of sufficient resolution. It is well-known that the resolution of the finest scales in turbulent thermal transport is strongly dependent on the dissipation rates of turbulent kinetic energy and of the temperature variance. The volume averaged turbulent thermal dissipation rate can be expressed as hT i ¼ hað@T 0 =@xi Þ2 i where T 0 are the temperature fluctuations over the mean and the angled bracket suggests a volume averaging procedure. The turbulent kinetic energy dissipa 2 0 0 tion rate can be expressed as hi ¼ hm c_ij i where c_ij is the strain rate component based on the fluctuating velocity field and m is the kinematic viscosity. The mean kinetic energy dissipation rate hi is related to the mean Kolmogorov scale hgK i which is the smallest length scale of turbulence when m a (i.e. Pr 1). In case of m > a (i.e. Pr > 1), the smallest length scale is determined by the dissipation rate of scalar variance, and is represented by the Batchelor scale hgB i. These length scales are defined as:
hgK i ¼
m3=4 hgK i ; hgB i ¼ pffiffiffiffiffi Pr hi1=4
ð9Þ
Thus, DNS becomes excessively expensive for Pr 1, which is often encountered for liquid flows. Therefore, several alternative grid spacing criteria for the DNS of turbulent Rayleigh–Bénard convection have been proposed in the existing literature [20,21]. Shishkina et al. [21] proposed a grid spacing criterion for the DNS of turbulent Rayleigh–Bénard convection of Newtonian
4
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
(i.e. Bn ¼ 0) fluids. Accordingly, the permissible non-dimensional mesh size (Dx=H) can be expressed in terms of both hgK i and hgB i in the following manner:
Dx hgK i ¼ H H
Pr2 ðNu 1ÞRa
Dx hgB i ¼ H H
!1=4
1 Nu 1 Ra
for
Pr 1
ð10iÞ
for
Pr > 1
ð10iiÞ
!1=4
Shishkina et al. [21] expressed that the global measures of thermal transport, such as Nusselt number, time-averaged temperature profiles and volume-averaged dissipation rates, are mostly insensitive to moderate under-resolution of hgB i for turbulent Rayleigh–Bénard convection of Newtonian fluids. Dropkin and Somerscales [22] also experimentally analysed turbulent Rayleigh–Bénard convection of Newtonian fluids (i.e. Bn ¼ 0) in cubic enclosure for 5 104 Ra 7:17 108 and 0:02 Pr 11560. They proposed the mean Nusselt number correlation given by:
Nu ¼ 0:069 Ra 1=3 Pr 0:074
ð11iÞ
Similarly, Xia et al. [23] experimentally analysed turbulent Rayleigh–Bénard convection of Newtonian fluids in cubical enclosures for 2 107 Ra 3 1010 and 4 Pr 1350. They proposed the alternative mean Nusselt number correlation given by:
Nu ¼ 0:14 Ra 0:297 Pr 0:03
ð11iiÞ
In the current analysis, the criteria given in Eq. (10) by Shishkina et al. [21] have been used along with the Nu correlation of Refs. [22,23] for Ra ¼ 107 and Ra ¼ 108 , respectively for the purpose of estimating the grid resolution for both Newtonian fluid (i.e. Bn ¼ 0) and Bingham fluid cases. It is worth noting that the spatial resolution will be sufficient for Bingham fluid cases as long as it is sufficient for the Newtonian (i.e. Bn ¼ 0) fluid. The effective viscosity leff ¼ l þ sy =c_ in Bingham fluids increases with increasing Bingham number for a given value of plastic viscosity l and thus the effective Rayleigh number Raeff ¼ qgbðT H T C ÞH3 =leff a for Bingham fluid cases assumes much smaller values than the nominal value of Rayleigh numberRa. Thus, the effects of instability are much weaker in Bingham fluid cases with high values of Bingham number than in Newtonian (i.e. Bn ¼ 0) fluid cases with same nominal values of Ra and Pr. Based on these criteria, uniform Cartesian grids of 1503 and 2303 have been used for the analysis of Ra ¼ 107 and Ra ¼ 108 cases, respectively. It has been checked a-posteriori that the pffiffiffiffiffi chosen grid resolution was sufficient since yþ Pr (where ffiffiffiffiffiffiffiffiffiffiffi p yþ ¼ qus y=l with us ¼ sw =q, sw and y being the friction velocity, wall shear stress and wall normal distance of the wall adjacent grid point, respectively) remains smaller than unity for both mesh configurations in the case of Ra ¼ 10
7
and
Ra ¼ 108 : 3.3. Benchmarking and validation The present numerical scheme has been validated with DNS data of turbulent Rayleigh–Bénard convection of Newtonian (i.e. Bn ¼ 0) fluids in a cubic enclosure with respect to the results reported by Kaczorowski et al. [24] at Ra ¼ 107 , Ra ¼ 108 and Pr ¼ 4:38. Additionally, Vasco et al. [25] analysed laminar natural
convection of Bingham fluids in cubic enclosure with differentially heated sidewalls for 104 Ra 106 and Pr ¼ 1. The current numerical scheme has been validated with respect to the findings by Kaczorowski et al. [24] and Vasco et al. [25] for Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases, respectively, and an excellent agreement (less than 2% uncertainty on mean Nusselt number) has been obtained, as shown in Table 1. Additionally, Xu et al. [26] presented the mean Nusselt number results which are defined based on the volume and time average of D E @u @u the kinetic energy i:e: heiV ;t ¼ v @xij @xij and thermal dissiV ;t D E @T @T in addition to Reynolds pation rate i:e: heT iV;t ¼ a @x i @xi V ;t
number based on the root-mean-square (rms) velocity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i:e: Rerms ¼ ð hux 2 þ uy 2 þ uz 2 iV ;t HÞ=v for the same configuration at Ra ¼ 107 and Pr ¼ 7. Accordingly, for additional validation, the comparison of the current study and Xu et al. [26] can be seen in Table 2, which shows that the current results agree very well with the data reported in Xu et al. [26]. 4. Results and discussion The effects of Rayleigh and Bingham number on the global flow behaviour and heat transfer characteristics are presented first. This is followed by a detailed discussion of the velocity and temperature spectra for all cases. The results presented in this section are obtained when the quasi-steady state has been reached. In order to demonstrate the quasi-steady state behavior, the temporal evolution of the mean Nusselt number Nu with normalised time is shown in Fig. 2 for Ra ¼ 107 and Pr ¼ 320. Under the quasi-steady state, the mean values of Nu for top and bottom walls become identical to each other, provided the oscillations due to the unsteady flow are neglected in this context. Zhang et al. [5] showed that Bingham fluid flows are unconditionally linearly stable for the quiescent initial condition. Thus, the quasi-steady Newtonian (i.e. Bn ¼ 0) fluid solutions for a given set of values of nominal Rayleigh and Prandtl numbers are used as the initial conditions for Bingham fluid simulations in the current analysis. Since it is possible to modify the yield stress by applying electrical or magnetic fields in magneto- and electro-rheological fluids, the quasi steady-state Newtonian solutions can be considered as realistic initial conditions with practical relevance for the Bingham fluid simulations. 4.1. Effects of variations of nominal Rayleigh number The three-dimensional views of non-dimensional isotherms (i.e. h ¼ ðT T C Þ=ðT H T C Þ) are shown in Fig. 3a for Ra ¼ 107 and 108 at Bn ¼ 0:43 and Pr ¼ 320. It should be noted that 0.05% by weight Table 1 Comparison of the mean Nusselt number for both Newtonian [24] and Bingham [25] fluid cases with the benchmark data. Newtonian (Bn ¼ 0) fluid case Kaczorowski et al. [24]
Present Study
7
Nu (Ra ¼ 10 , Pr ¼ 4:38)
16.1
16.3
Nu (Ra ¼ 108 , Pr ¼ 4:38)
32.2
32.8
Vasco et al. [25]
Present Study
8.946
8.816
Bingham (Bn ¼ 0:1) fluid case
6
Nu (Ra ¼ 10 , Pr ¼ 1:0)
5
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947 Table 2 Comparison of DNS results of the Reynolds and Nusselt numbers between Xu et al. [26] and present study at Ra ¼ 107 and Pr ¼ 7.
Xu et al. [26] Present
Nu Eq. (8)
Nue ¼
16.19 16.16
16.11 16.10
heiV;t H4 Pr2 Rav 3
þ1
Nue ¼
Rerms
heT iV ;t H2
aðDT Þ2
16.03 15.95
86.41 86.22
4.2. Effects of the variation of nominal Bingham number The distributions of h and W are shown in Fig. 5 for different Bn at Ra ¼ 108 and Pr ¼ 320. It is also worth noting that two addi-
Fig. 2.pVariation of Nu with non-dimensional simulation time (i.e. ffiffiffiffiffiffiffiffiffiffi t ¼ ðt a RaPrÞ=H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall at Ra ¼ 107 , Bn ¼ 0 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Carbopol-water solution at Ra ¼ 107 and 108 yields Bn ¼ 0:43: The Bn value are obtained based on sy ¼ 0:007 Pa, q ¼ 961 kg=m3 ,
l ¼ 0:046 kg=ðmsÞ, a ¼ 1:5 107 m2 =s and b ¼ 2 104 1=K (which are taken from Table 1 of Ref. [2]) and different values of H and (T H T c ) according to Eq. (6). It can be seen from Fig. 3a that the iso-surfaces of h are more distorted for Ra ¼ 108 than in 7
Ra ¼ 10 , which is indicative of the strengthening of convection with increasing Ra due to the relative enhancement of buoyancy forces in comparison to the viscous resistance. This can be also confirmed by the path lines coloured by non-dimensional velocity pffiffiffiffiffiffiffiffi magnitude (i.e. W ¼ ui ui H=a) for Ra ¼ 107 and 108 , which are shown in Fig. 3b at Bn ¼ 0:43 and Pr ¼ 320. Fig. 3b indicates that the magnitude of the non-dimensional velocity increases with an increase inRa, which is indicative for the strengthening of convective transport. It also shows the Ra ¼ 108 flow field exhibits smaller flow structures compared to Ra ¼ 107 : The strengthening of convective thermal transport with increasing Ra can also be substantiated from the temporal evolutions of Nu for top and bottom walls in Fig. 4 for different Ra at Bn ¼ 0:43 and Pr ¼ 320. Fig. 4 shows that Nu increases dramatically with increasing Ra as a result of augmented convective thermal transport. The AURs for Ra ¼ 107 and 108 at Bn ¼ 0:43 and Pr ¼ 320 are shown in Fig. 3c (see the grey regions). It is evident from Fig. 3c that the size and the probability of finding AURs decrease with increasing Ra due to diminished viscous resistance in the flow. In the context of bi-viscosity regularisation, the flow does not stop in a true sense in AURs but the flow within the islands of AURs is too weak to influence the heat transfer rate and thus the exact shape and size of AURs do not affect the heat transfer rate and the mean Nusselt number.
tional Bingham numbers (i.e. Bn ¼ 0:215 and 2.15) at Ra ¼ 108 are chosen for the purpose of a numerical experiment. For the sake of brevity the results for Bn ¼ 0:215 are not always shown explicitly due to the monotonic effect of Bn. It can be seen from Fig. 5a that the iso-surfaces become less deformed and thermal plumes from the bottom wall become less apparent as Bn increases. These are the indications of the relative weakening of buoyancy force due to the augmented flow resistance as a result of increased yield stress. This also leads to the weakening of convection with increasingBn, which can be substantiated from Fig. 5b where the path lines colored by non-dimensional velocity magnitude for the corresponding cases are shown. Fig. 5b shows that the non-dimensional velocity magnitude decreases with increasing Bn due to the weakening of the convective thermal transport. The temporal evolutions of Nu for top and bottom walls are shown in Fig. 6 for three different Bn values at Ra ¼ 108 and Pr ¼ 320. All mean Nusselt numbers obtained for different values of Bn for Ra ¼ 107 and Ra ¼ 108 are summarised in Table 3. Eq. (11) suggests a scaling of mean Nusselt number as Nu Ran Pr m but this formula does not take into account the variable fluid properties of Bingham fluids. Introducing the volume and time averaged effective viscosity hleff iV;t , effective Rayleigh and Prandtl numbers can be defined in the following manner:
Raeff ¼ RaBn¼0
l hleff iV;t
;
Preff ¼ PrBn¼0
hleff iV;t
l
ð12Þ
Consequently, it can be expected that the Nusselt number for Bingham fluids NuBn possibly can be scaled with the nominal (Newtonian) Nusselt number NuBn¼0 and the effective viscosity as:
NuBn NuBn¼0
!p
l hleff iV;t
:
ð13Þ
Fig. 7 demonstrates that Eq. (13) approximates the mean Nusselt number for Bingham fluids reasonably well using the concept of effective viscosity and an exponent of p ¼ 0:13: Here, the effective viscosities have been determined after conducting the simulations. The ratio of l=hleff iV;t is expected to be unity for Newtonian fluids (Bn ¼ 0) and to approach small values for larger Bingham numbers when the fluid is unyielded. An analysis of the existing data shows that the ratio of l=hleff iV;t can be reasonably approximated by the following expression:
l=hleff iV;t ¼
1 min ð1 þ b Bnm ; 104 Þ
ð14Þ
with m 2:5 and b e6 ðb e5 Þ for Ra ¼ 107 (Ra ¼ 108 ) where e is Euler’s number. However, more data will be needed to develop a generally valid expression. The far right circle (square) in Fig. 7 represents the Newtonian case for Ra ¼ 107 (Ra ¼ 108 ) and the remaining Nusselt numbers decrease with increasing Bn. The weakening of convective transport, as reflected in the reduction of the mean Nusselt number Nu, with increasing Bn is
6
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Fig. 3. (a) Distributions of iso-surfaces of non-dimensional temperature h, (b) The path lines are coloured by non-dimensional velocity magnitude (i.e. W), (c) Distribution of AURs for different values of Ra at Bn ¼ 0:43 and Pr ¼ 320. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(see Fig. 6). For very large values of Bingham number, the boundary layer thickness becomes of the order of the enclosure size H. At that stage, the fluid flow does not influence the thermal transport and heat transfer takes place principally due to conduction, which is reflected in the unity mean Nusselt number. The distributions of magnitude of normalised shear stress (i.e. pffiffiffiffiffiffiffiffiffiffi s ¼ sij sij H2 =ðqa2 Þ) values on the x z plane at y ¼ H are shown
Fig. 4. Temporal variation of Nu with non-dimensional simulation time (i.e. pffiffiffiffiffiffiffiffiffiffi t ¼ t a RaPr =H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall for different Ra at Bn ¼ 0:43 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
consistent with the distributions of AURs shown in Fig. 8a for different values of Bnat Ra ¼ 108 and Pr ¼ 320. It can be seen from Fig. 8a that the size and the probability of finding AURs increase with increasing Bn due to the weakening of convective transport, and this gives rise to a reduction in the mean Nusselt number
in Fig. 8b for different Bn values at Ra ¼ 108 and Pr ¼ 320. The Rayleigh-Bénard cells can be seen for the Newtonian (i.e. Bn ¼ 0) case and also for small values of Bn (e.g. Bn ¼ 0:43). However, the Rayleigh-Bénard cells do not clearly appear for large values of Bn (i.e. Bn ¼ 2:15) due to the diminished strength of the thermal convection. 4.3. Velocity and temperature spectra It can clearly be seen from Figs. 5 and 8 that the flow structures are strongly affected by the Bingham number. Consequently, it will be of interest to characterise the range of flow structures in terms of the energy spectra of temperature and velocity fluctuations. In this context, it is important to mention that determination of universal scaling laws for inertial range two-point statistics, as done
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
7
Fig. 5. (a) Distributions of iso-surfaces of non-dimensional temperature h, (b) The path lines are coloured by non-dimensional velocity magnitude (i.e. W) for different values of Bn at Ra ¼ 108 and Pr ¼ 320. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
pffiffiffiffiffiffiffiffiffiffi Fig. 6. Temporal variation of Nu with non-dimensional simulation time (i.e. t ¼ ta RaPr =H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall for different Bn at Ra ¼ 108 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
by Kolmogorov in ideal hydrodynamics turbulence, is a long-standing challenge in Rayleigh-Bénard turbulence, even in the context of Newtonian flows [27] and the same holds true for scalar spectra of high Prandtl (Schmidt) number fluids.
In the inertial range, Bolgiano [28] and Obukhov [29] proposed a dual cascade: for j small wave numbers in the range LB < j1 (where LB is the Bolgiano scale) they predicted a dominance of the buoyancy force over the inertial force leading to velocity and
8
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Table 3 Summary of the mean Nusselt number obtained in the current analysis for both Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases at Pr ¼ 320.
Nu
Ra ¼ 107
Ra ¼ 108
Bn ¼ 0 Bn ¼ 0:215 Bn ¼ 0:43 Bn ¼ 2:15
17.43 14.25 11.095 1.0
32.8 27.79 23.35 10.20
Fig. 7. Mean Nusselt number based on the DNS results versus the Nusselt number according to Eq. (13) using an exponent of p ¼ 0:13 for Ra ¼ 107 and 108 and different Bingham numbers.
velocity [30,38] and hence the applicability of this hypothesis is questionable in this configuration. In the light of the above discussion and the remaining uncertainties in regards of universal scaling laws, it will be interesting to analyse velocity and temperature spectra for the flow configurations considered in this work. The present configuration does not have homogeneous directions and hence, despite the aforementioned difficulties, the temporal spectra will be presented. It has been decided not to transform the spectra into the equivalent spatial domain. Nevertheless, a suitable mean velocity scale is needed for determining the Nyquist frequency in the Power Spectral Densities (PSD), which are subsequently discussed. The WienerKhinchin theorem states that the PSD (u0 ) is identical to the fluctuating kinetic energy spectrum in case of homogeneous flows. According to the Nyquist theorem, the maximum detectable frequency in a signal is given by f max ¼ 1=ð2DtÞ where Dt is the sampling interval, which corresponds in the present case to the constant time step in the simulation. In Computational Fluid Dynamics (CFD), Dt typically has an upper bound determined by the Courant number. However, the time step can be arbitrarily small which would result in very high resolvable frequencies. It is intuitively clear that for a given spatial resolution the resolved frequency range cannot be arbitrary large. All simulations have been run with a Courant number smaller than unity. Accordingly, we determine the Nyquist frequency using Taylors hypothesis as f max ¼ 1=ð2Dx=U con Þ where U con is a suitable convection velocity. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Table 4 shows the mean velocity (i.e. U mean ¼ ux 2 þ uy 2 þ uz 2 ), p ffiffiffiffiffiffiffiffiffiffiffi 0 the fluctuation velocity scale (i.e. u ¼ 2=3k) where k is the h i 2 turbulent kinetic energy (i.e. k ¼ 0:5 ðux 0 Þ2 þ uy 0 þ ðuz 0 Þ2 ) and 0
temperature spectra as j11=5 and j7=5 , respectively. For intermediate wave numbers, the temperature field evolves as a passive scalar, and both the velocity and temperature spectra depict a j5=3 range. As reviewed in Boffetta and Mazzino [27], the Bolgiano scale LB coincides with the integral scales L in three dimensions. Further, according to Lohse and Xia [30], there is a fundamental lack of clear scale separation between LB and L which makes the identification of such a regime difficult in practice. In the high wavenumber range, for high Prandtl numbers, there are two dissipation scales, which are the Kolmogorov scale hgK i and the Bachelor scale hgB i [31], respectively. For high Prandtl numbers one obtains hgB i < hgK i and Batchelor [32] predicted that the scalar spectrum should display a j1 dependence in the so called viscous-convective subrange hgB i < j1 < hgK i. In order to better distinguish the spectra for high and low Prandtl numbers, reference simulations with Pr ¼ 1 will be presented as well in this section. According to the review of Warhaft [33], the experimental evidence has been elusive and interested readers are referred to references in Ref. [33] for further information. A detailed discussion has been provided as well by Antonia and Orlandi [34]. In this context, it is worth mentioning that a modified scaling of the form j1p with p Re1=2 has been suggested by Kerstein [35] and that the location of the j1 range has sometimes been observed for j1 > hgKi (Ref. [36]), which is inconsistent with the classical expectation. The lack of resolution especially in early measurements could be an issue, which might explain why the j1 scaling is not well-established [36]. Another problem is that typically time domain results are compared to theoretical predictions for spatial domains by invoking Taylor’s frozen-flow hypothesis [37], which requires the velocity fluctuations to be much smaller than the mean flow velocity. However, it is well-known that the rootmean-square (i.e. rms) velocity fluctuation in Rayleigh-Bénard convection is often of the same order of magnitude as that of the mean
fluctuation intensity (u =U mean ) in the centre of the domain (i.e. x=H ¼ y=H ¼ z=H ¼ 0:5). In order to compare the spectra, two additional simulations have been performed for Newtonian (i.e. Bn ¼ 0) fluids at Pr ¼ 1:0 for Ra ¼ 107 and 108 . It can be seen from Table 4 0 that the ratio u =U mean remains of the order of 75 3% for the 0 Pr ¼ 320 case, whereas this ratio is u =U mean 63% for Pr ¼ 1. This demonstrates the well-known fact that even in the centre of the domain the applicability of Taylor’s hypothesis remains questionable. Hence, in order to determine the Nyquist frequency, U con ¼ U mean u0 has been used. From the foregoing, it becomes clear that the cut-off frequency is not a precise number but rather can be attributed to a certain range. The same holds true for the frequencies corresponding to Kolmogorov and Bachelor scales, which will henceforth be denoted as fhgKi = Ucon/(2hgKi) and fhgBi = Ucon/(2hgBi) and have been estimated using Eqs. (11i) and (11ii), respectively. Fig. 9 shows the PSD of temperature and velocity fluctuation 0
0
(i.e. T and u respectively) for different Pr values at Ra ¼ 107 in the Newtonian (i.e. Bn ¼ 0) case in the centre of the domain based on single probing point. According to Bachelor [28], one would expect a f
1
(i.e.
j1 ) scaling in the range f hgK i < f < f hgB i
1
(i.e. hgB i < j < hgK i) for the temperature spectra for Pr ¼ 320, whereas the velocity spectra should decay much faster in the dissipation range for f > f hgK i . The estimated lower and upper limits of f hgK i are shown with dashed and solid blue lines in Fig. 9, whereas the corresponding Bachelor scales are shown in black. It is important to note that the cases with high Pr and Bn > 0 do not feature a fully developed turbulent flow field. Hence the scaling laws are shown for the purpose of reference. The PSD (T 0 ) for Pr ¼ 320 depicts a decay close to f
1
for f f hgK i but this range does not
reach up to f hgB i . An inertial range scaling for f < f hgK i is not clearly visible. In contrast for Pr ¼ 1, there is a small inertial range for f < f hgK i , whereas the PSD (T 0 ) decays much faster for f > f hgK i : As
9
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Fig. 8. (a) Distribution of three-dimensional AURs, (b) distribution of the magnitude of normalised shear stress (s ¼ s H2 =ðqa2 Þ) on the cold top wall (i.e. x z plane view at y ¼ H) for different Bn at Ra ¼ 108 and Pr ¼ 320.
Table 4 0 0 Summary of the mean velocity (i.e. U mean ), fluctuation velocity scale (i.e. u ) fluctuation intensity (u =U mean ) in the centre of the domain (i.e. x=H ¼ y=H ¼ z=H ¼ 0:5). Cases
0
U mean
u ¼ 3
2.0397 10
Ra ¼ 107 ; Pr ¼ 1:0; Bn ¼ 0
pffiffiffiffiffiffiffiffiffiffiffi 2=3k
0
u =U mean 3
0.634
3
1.2939 10
3
7
5.5003 10
4.1529 10
0.755
8
3.1556 103
1.9778 103
0.627
8
3
8.9840 10
3
6.3900 10
0.711
8
5.7679 103
4.4898 103
0.778
8
3
3
0.722
Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 1:0; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0:215
4.9757 10
Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0:43
expected, the velocity spectra for both Prandtl numbers depict a dissipation range behaviour for scales smaller than the Kolmogorov scales. However, the decay is considerably slower for Pr ¼ 320 which might be a consequence of the slower decay of 0
the temperature fluctuations (i.e. T ) and the active scalar nature of this quantity (back-coupling in the momentum equation via the buoyancy term). The PSD (T 0 ) and PSD (u0 ) look qualitatively similar for Ra ¼ 108 , which is shown in Fig. 10 for Pr ¼ 1 and Pr ¼ 320 in the centre of the domain based on single probing point. Again the PSD (u0 ) decays slower for Pr ¼ 320; possibly due to the coupling between momentum and temperature equation. Finally, it is worth nothing that the features of a viscous-convective range are slightly more pronounced for the Ra ¼ 108 case with Pr ¼ 320. As thermal convection weakens with decreasing (increasing) Ra (Bn), fluid flow starts to exhibit a nearly laminar behaviour for small (high) Ra (Bn) values. This can be easily seen from the
3.5938 10
iso-surfaces of non-dimensional temperature in Fig. 5. Therefore, the discussion of spectra of the Bingham fluids has been limited to the case with high Ra and small Bn values (i.e. Bn ¼ 0:215 and 0.43 for Ra ¼ 108 ). Accordingly, the PSD of rms velocity fluctuation 0
(i.e. u ) and temperature (i.e. T 0 ) for different Bn values at Ra ¼ 108 and Pr ¼ 320 in the centre of the domain based on single probing point are shown in Fig. 11. It should be noted here that the spatial cut-off frequencies based on mean Kolmogorov f hgK i and Batchelor f hgB i length scales have been calculated based on the criterion in Eq. (10) given by Shishkina et al. [21] but using the concept of effective Rayleigh and Prandtl numbers introduced in Eq. (12). As effective viscosity is expected to be higher for Bn > 0, this results in a smaller (higher) effective Rayleigh (Prandtl) numbers. For instance, hleff iV ;t for the Bingham fluid cases (i.e. Bn ¼ 0:215 and Bn ¼ 0:43) was obtained around 3.5 and 27 times bigger than l in the Newtonian case (i.e. Bn ¼ 0) respectively. This implies that
10
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
0
Fig. 9. Power spectral density of (a) temperature (i.e. T ), (b) rms velocity fluctuations (i.e. u ) for Ra ¼ 107 for Newtonian (i.e. Bn ¼ 0) fluid cases in the centre of the domain based on single probing point for different Pr values. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and Batchelor fhgBi 0 0 length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u and U mean u , respectively.
effective values of Rayleigh (Prandtl) numbers based on hleff iV;t values are 3.5 and 27 times smaller (bigger) for the Bingham fluid cases (i.e. Bn ¼ 0:215 and Bn ¼ 0:43) than their respective nominal values. Accordingly, this yields larger Kolmogorov and Batchelor scales or alternatively, smaller values of f hgK i and f hgB i . Thus, it can be observed from Fig. 11 that for the Bingham fluid case similar observations hold true as for the Newtonian fluid (i.e. Bn ¼ 0) with Pr ¼ 320. In particular, there is a small range where viscousconvective behaviour can be observed in the PSD (T 0 ). In order to facilitate easier observation of the viscousconvective range like scaling, Fig. 12 shows the power spectral density of temperature fluctuations multiplied by frequency 0 (i.e. PSD T f ) for two different Ra at Pr ¼ 1 (black) and Pr ¼ 320 (blue) for the Newtonian (i.e. Bn ¼ 0) fluid case as well as for Ra ¼ 108 , Pr ¼ 320, Bn ¼ 0:215 (green) and Ra ¼ 108 , Pr ¼ 320, Bn ¼ 0:43 (red) for Bingham fluid cases in the centre of the domain based on single probing point. Fig. 12 clearly shows that despite having the same nominal Rayleigh and Prandtl num-
bers, the Bingham fluid cases (red and green line) exhibit a different PSD with more energy on large scales compared to the Newtonian fluid (blue line) due to the smaller effective Ra. Further, it can be observed that temperature fluctuation spectra multiplied by frequency decay considerably faster than f
1
for f f hgK i in the
Newtonian (i.e. Bn ¼ 0) fluids (for both Ra ¼ 107 and Ra ¼ 108 ) at Pr ¼ 1, whereas a small plateau tends to appear for Newtonian (Bn ¼ 0 for Ra ¼ 108 ) and Bingham (Ra ¼ 108 , Bn ¼ 0:215 and Bn ¼ 0:43) fluid cases at Pr ¼ 320 at for f f hgK i . This viscousconvective range is most pronounced for the Newtonian case with the highest Rayleigh number for Pr ¼ 320. 5. Conclusions Direct Numerical Simulations (DNS) have been conducted to analyse three-dimensional Rayleigh-Bénard convection (differentially heated horizontal walls with heated bottom wall, i.e. heated from below and cooled from above, whereas the other walls are
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
11
0
Fig. 10. Power spectral density of (a) temperature (i.e. T ), (b) rms velocity fluctuations (i.e. u ) for Ra ¼ 108 for Newtonian (i.e. Bn ¼ 0) fluid cases in the centre of the domain based on single probing point for different Pr values. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and Batchelor fhgBi 0 0 length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u and U mean u , respectively.
adiabatic) of yield stress fluids obeying a Bingham model in cubic domains for high values of nominal Rayleigh number (107 and 108 Þ for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution). It has been found that the mean Nusselt number Nu increases with increasing values of Rayleigh number for both Newtonian and Bingham fluids. However, weaker convective transport in Bingham fluid cases leads to smaller values of Nu than that obtained in the case of Newtonian fluids with the same nominal value of Rayleigh numberRa and a scaling based on volume averaged effective viscosities has been suggested. The flow domain becomes increasingly unyielded and the effects of Rayleigh-Bénard instability become less pronounced as the Bingham number increases. For Bingham fluids, thermal convection decreases considerably with increasing Bingham number. Nevertheless, a small viscous1
convective type range with a f scaling has been observed for the temperature fluctuation spectra, roughly located close to the
Kolmogorov scale, in all cases for Pr ¼ 320. It has been observed that for identical nominal Rayleigh numbers, the Bingham fluid cases exhibit a smaller range of scales with more energy towards the large scales when compared to the corresponding Newtonian fluid case. By comparing simulations for unity Prandtl number and Pr ¼ 320, it has been found that the decay of turbulent kinetic energy in the dissipation range is considerably slower for Pr ¼ 320, which has been attributed to the coupling between temperature and momentum equation combined with the slower decay of thermal variance. To the best of the authors’ knowledge, this is the first analysis of Rayleigh-Bénard convection of yield stress fluids at high nominal Rayleigh numbers. Future investigations are aimed at the analysis of similar configurations at even higher values of Ra: Establishing correlations for the mean Nusselt number for complex nonNewtonian fluids will be of great importance for a variety of technical applications.
12
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
0
Fig. 11. Power spectral density of rms velocity fluctuations (i.e. u ) on the top and temperature fluctuations (i.e. T ) on the bottom, for different Bn values at Ra ¼ 108 and Pr ¼ 320 in the centre of the domain based on single probing point. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and 0 Batchelor fhgBi length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u 0 and U mean u , respectively.
0
Fig. 12. Power spectral density of temperature multiplied by frequency (i.e. PSD ðT Þ f ) for different Ra at Pr ¼ 1 (black), Pr ¼ 320 (blue) for Newtonian (i.e. Bn ¼ 0) fluid case andRa ¼ 108 , Pr ¼ 320 and Bn ¼ 0:215 (green) and Ra ¼ 108 , Pr ¼ 320 and Bn ¼ 0:43 (red) for Bingham fluid cases in the centre of the domain based on single probing point. The mean Kolmogorov fhgKi length scale for all cases is shown by the vertical lines with the respective colours. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgements The authors gratefully acknowledge that N8HPC (North East UK national High-Performance Computing) and Rocket HighPerformance Computing (in-house cluster of Newcastle University) have been used for the current analysis. Computer resources for this project have also been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant no. pn56di. References [1] M. Darbouli, C. Métivier, J.M. Piau, A. Magnin, A. Abdelali, Rayleigh- Bénard convection for viscoplastic fluids, Phys. Fluids 25 (2013) 023101. [2] Z. Kebiche, C. Castelain, T. Burghelea, Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid, J. Nonnewton. Fluid Mech. 203 (2014) 9–23. [3] J.E. Stangroom, Electrorheological fluids, Phys. Technol. 14 (1983) 290–296. [4] N. Wereley, Magnetorheology: Advances and Applications, Royal Society of Chemistry, 2013. [5] J. Zhang, D. Vola, I.A. Frigaard, Yield stress effects on Rayleigh-Bénard convection, J. Fluid Mech. 566 (2006) 389–419. [6] N.J. Balmforth, A.C. Rust, Weakly nonlinear viscoplastic convection, J. Nonnewton. Fluid Mech. 158 (2009) 36–45. [7] A. Vikhansky, On the onset of natural convection of Bingham liquid in rectangular enclosures, J. Nonnewton. Fluid Mech. 165 (2010) 1713–1716. [8] O. Turan, S. Yigit, N. Chakraborty, Critical condition for Rayleigh-Bénard convection of Bingham fluids in rectangular enclosures, Int. Commun. Heat Mass Transf. 86 (2017) 117–125. [9] O. Turan, N. Chakraborty, R.J. Poole, Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure, J. of Nonnewton. Fluid Mech. 171 (2012) 83–96. [10] O. Turan, R.J. Poole, N. Chakraborty, Boundary condition effects on natural convection of Bingham fluids in a square enclosure with differentially heated horizontal walls, J. Comput. Therm. Sci. 4 (2012) 77–97. [11] S. Yigit, R.J. Poole, N. Chakraborty, Effects of aspect ratio on natural convection of Bingham fluids in rectangular enclosures with differentially heated horizontal walls heated from below, Int. J. Heat Mass Transf. 80 (2015) 727– 736. [12] S. Yigit, R.J. Poole, N. Chakraborty, Laminar natural convection of Bingham fluids in inclined differentially heated square enclosures subjected to uniform wall temperatures, J. Heat Transf. 137 (2015) 052504. [13] S. Yigit, S. Chen, P. Quinn, N. Chakraborty, Numerical investigation of laminar Rayleigh-Bénard convection of Bingham fluids in square cross-sectioned cylindrical enclosures, Int. J. Therm. Sci. 110 (2016) 356–368. [14] S. Yigit, N. Chakraborty, Influences of aspect ratio and wall boundary condition on laminar Rayleigh-Bénard convection of Bingham fluids in rectangular enclosures, Int. J. Numer. Meth. Heat Fluid Flow 27 (2017) 310–333.
13
[15] E. O’Donovan, R. Tanner, Numerical study of the Bingham squeeze film problem, J. Nonnewton. Fluid Mech. 15 (1984) 75–83. [16] E. Mitsoulis, T. Zisis, Flow of Bingham plastics in a lid-driven square cavity, J. Nonnewton. Fluid Mech. 101 (2001) 173–180. [17] I. Karimfazli, I.A. Frigaard, A. Wachs, A novel heat transfer switch using the yield stress, J. Fluid Mech. 783 (2015) 526–566. [18] R.R. Huilgol, G.H.R. Kefayati, Natural convection on Bingham fluid using the operator splitting method, J. Nonnewton. Fluid Mech. 220 (2015) 22–32. [19] O. Turan, N. Chakraborty, R.J. Poole, Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls, J. Nonnewton. Fluid Mech. 165 (2010) 901–913. [20] G. Grötzbach, Spatial resolution requirements for direct numerical simulation of the Rayleigh-Bénard convection, J. Comp. Phys. 49 (1983) 241–264. [21] O. Shishkina, R.J.A.M. Stevens, S. Grossmann, D. Lohse, Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution, New J. Phys. 12 (2010) 075022. [22] D. Dropkin, E. Somerscales, Heat transfer by natural convection in liquids confined by two parallel plates which are inclined at various angles with respect to the horizontal, J. Heat Transf. 87 (1965) 1–6. [23] K.-Q. Xia, S. Lam, S-Q, Zhou, Heat-flux measurement in high-Prandtl-number turbulent Rayleigh-Bénard convection, Phys. Rew. Lett. 88 (2002) 064501. [24] M. Kaczorowski, C. Kai-Leong, X. Ke-Quing, Turbulent flow in the bulk of Rayleigh-Bénard convection: aspect-ratio dependence of the small-scale properties, J. Fluid Mech. 747 (2014) 73–102. [25] D.A. Vasco, N.O. Moraga, G. Haase, Parallel finite volume method simulation of three dimensional fluid flow and convective heat transfer for viscoplastic nonNewtonian fluids, Num. Heat Transf. Part A 66 (2014) 990–1019. [26] Xu. Ao, Le Shi, Heng-Dong Xi, Lattice Boltzmann simulations of threedimensional thermal convective flows at high Rayleigh number, Int. J. Heat Mass Transf. 140 (2019) 359–370. [27] G. Boffetta, A. Mazzino, Incompressible Rayleigh-Taylor turbulence, Annu. Rev. Fluid Mech. 49 (2017) 119–143. [28] R. Bolgiano, Turbulent spectra in a stably stratified atmosphere, J. Geophys. Res. 64 (1959) 2226–2229. [29] A.M. Obukhov, On the influence of Archimedean forces on the structure of the temperature field in a turbulent flow, Dokl. Akad. Nauk. SSR 125 (1959) 1246– 1248. [30] D. Lohse, K.Q. Xia, Small-scale properties of turbulent Rayleigh-Bénard convection, Annu. Rev. Fluid Mech. 42 (2010) 335–364. [31] H. Tennekes, J.L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, 1972. [32] G.K. Batchelor, Small-scale variation of convected quantities like temperature in a turbulent field. Part 1. General discussion and the case of small conductivity, J. Fluid Mech. 5 (1959) 113–133. [33] Z. Warhaft, Passive scalars in turbulent flows, Annu. Rev. Fluid Mech. 32 (2010) 203–240. [34] R.A. Antonia, P. Orlandi, Effect of Schmidt number on small-scale passive scalar turbulence, Appl Mech Rev. 56 (2003) 615–632. [35] A.R. Kerstein, Scaling properties of the viscous-convective scalar spectral subrange in turbulent jets, Phys. Fluids A 3 (1991) 1832–1834. [36] D. Bogucki, J.A. Domaradzki, P.K. Yeung, Direct numerical simulations of passive scalars with Pr>1 advected by turbulent flow, J. Fluid Mech. 343 (1997) 111–130. [37] G.I. Taylor, The spectrum of turbulence, Proc. R. Soc. Lond. A 164 (1938) 476– 490. [38] Q. Zhou, K. Xia, Comparative experimental study of local mixing of active and passive scalars in turbulent thermal convection, Phys. Rev. 77 (2008) 056312.
Contents lists available at ScienceDirect
International Journal of Heat and Mass Transfer journal homepage: www.elsevier.com/locate/ijhmt
Effects of Rayleigh-Bénard convection on spectra of viscoplastic fluids Sahin Yigit a, Josef Hasslberger a, Nilanjan Chakraborty b,⇑, Markus Klein a a b
Instute of Applied Mathematics and Scientific Computing, Bundeswehr University Munich, Germany School of Engineering, Newcastle University, Newcastle upon Tyne, United Kingdom
a r t i c l e
i n f o
Article history: Received 19 June 2019 Received in revised form 21 October 2019 Accepted 22 October 2019
Keywords: Direct Numerical Simulation Rayleigh-Bénard convection Non-Newtonian Yield stress fluid High nominal Rayleigh number
a b s t r a c t Direct Numerical Simulations (DNS) have been conducted to analyse velocity and temperature spectra in unsteady three-dimensional Rayleigh-Bénard convection (differentially heated horizontal walls, i.e. heated from below and cooled from above whereas the other walls are adiabatic) at high nominal Rayleigh number of viscoplastic fluids obeying a Bingham model in cubic domains and the corresponding Newtonian cases with same nominal Rayleigh and Prandtl numbers for reference purposes. The simulations have been carried out for nominal Rayleigh numbers of 107 and 108 for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution). The effects of yield stress on velocity and temperature spectra and heat transfer characteristics of the flow fields have been investigated in terms of isotherms, path lines, apparently unyielded regions and mean Nusselt number, along with the power spectral densities of temperature and velocity fluctuations. It has been found that thermal convection strengthens (weakens) with increasing nominal Rayleigh number Ra (Bingham number Bn). Consistent with this observation, the range of flow structures becomes larger with higher (smaller) values of Ra (Bn). Additionally, it has been observed that the decay of the power spectral densities of temperature fluctuations in the viscous-convective range is close to the 1 inverse of frequency f (i.e. it scales as f ) for Pr > 1 in both Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases with moderate values of Bn. The fluid motion (including its unsteady features) weakens and the mean Nusselt number drops progressively with an increase in Bingham number. Ó 2019 Elsevier Ltd. All rights reserved.
1. Introduction Non-Newtonian fluids are immensely important in many industries such as polymer and chemical processing, biochemical synthesis. Yield stress fluids belong to a special type of non-Newtonian fluids, which flow once a critical stress (i.e. yield stress) is surpassed but act like a rigid solid below this threshold stress. Rayleigh-Bénard convection (i.e. the buoyancy-driven flow of a fluid heated from below and cooled from above) of yield stress fluids in enclosed spaces has a wide range of engineering applications, such as preservation of canned food, chemical processing, solar and nuclear power systems – to name only a few. Aqueous solutions of Carbopol are often used as model yield stress fluids in laboratory experiments [1,2]. Some magneto and electro-rheological fluids exhibit yield stress behaviour, where it is possible to modify the yield stress by applying electrical or magnetic fields [3,4]. This can be useful in reducing heat losses or designing new adaptive thermal systems to control convective thermal transfer in enclosed spaces. ⇑ Corresponding author. E-mail address: [email protected] (N. Chakraborty). https://doi.org/10.1016/j.ijheatmasstransfer.2019.118947 0017-9310/Ó 2019 Elsevier Ltd. All rights reserved.
Several recent studies [1,2,5–14] have analysed laminar Rayleigh-Bénard convection of yield stress fluids obeying a Bingham model in two-dimensional rectangular enclosures for small and moderate values of nominal Rayleigh number. These analyses revealed that the flow progressively weakens with increasing Bingham number (i.e. non-dimensional yield stress) due to additional flow resistance arising from the yield stress. This leads to a reduction in convective transfer rate, and thermal transport takes place principally due to conduction for large values of Bingham number. This behaviour arises because the fluid flow practically stops for large values of Bingham number and it does not influence the thermal transport under such a condition. A few studies [1,2,5–8] focused on critical conditions for the onset of Rayleigh-Bénard convection of Bingham fluids, whereas others [9–14] concentrated on heat transfer characteristics under conditions which are far beyond the critical value for the onset of fluid motion. It is worth noting that all the aforementioned studies in the existing literature [1,2,5–14] are restricted to laminar flow (i.e. nominal Rayleigh number from 103 to 105). However, Rayleigh-Bénard convection of yield stress fluids at high nominal values of Rayleigh number for which a steady laminar solution may not exist, has not yet been analysed in the existing literature
2
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
and consequently, little is known about the energy cascade in this kind of flows. This gap has been addressed here by conducting Direct Numerical Simulations (DNS) of unsteady threedimensional Rayleigh-Bénard convection of yield stress fluids obeying a Bingham model in cubic domains at high values of nominal Rayleigh number for which the Newtonian fluids (and also Bingham fluids of small values of Bingham number) are expected to exhibit turbulent natural convection. The simulations have been 7
8
carried out for nominal Rayleigh numbers of 10 and 10 for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution [2]). Accordingly, the main objectives of this study are: i) To demonstrate the effects of Rayleigh and Bingham numbers on Rayleigh-Bénard convection of yield stress fluids in cubic domains for high nominal Rayleigh numbers. ii) To provide physical explanations for the observed Rayleigh and Bingham number dependence of mean Nusselt number in the aforementioned configuration. iii) To analyse the effects of Rayleigh, Prandtl and Bingham numbers on velocity and temperature spectra at the centre of the flow field. The structure of the paper is as follows: Section 2 will introduce the governing equations and the physical background. The Numerical scheme is detailed in Section 3, followed by the discussion of the results in Section 4. The concluding remarks are presented at the end of the paper.
2.1. Modelling of the yield stress fluid The bi-viscosity regularisation proposed by O’Donovan and Tanner [15] has been used to mimic the stress-shear rate characteristics for Bingham fluids:
s ¼ sy c_ =c_ þ lc_ for c_ > sy =lyield
beyond
which
the
linear
stress-strain
rate
relation
(i.e.
s ¼ sy ðc_ =c_ Þ þ lc_ ) according to the Bingham model is invoked. Therefore, the precise shapes and sizes of Apparently Unyielded Regions (AURs) do not affect the mean Nusselt number in the present study. O’Donovan and Tanner [15] indicated that a value of lyield equal to 1000 l mimics the true Bingham model in a satisfactory manner but here lyield =l ¼ 104 is kept fixed to ensure higher fidelity of the simulations. This can also be confirmed by the mean Nusselt number comparison (i.e. Nu) for Ra ¼ 107 , Pr ¼ 320 and Bn ¼ 0:43, as Nu is obtained to be 11.093 and 11.095 for the values of
lyield =l ¼ 103 and 104 respectively.
Additionally, several studies analysed natural convection of Bingham fluids in a square cavity with differentially heated sidewalls obeying the Bingham model without using any regularisation [17,18] with numerical simulations. They compared the results obtained using the bi-viscosity regularisation [5] with the ones without regularisation [17,18]. It has been reported by Refs. [17,18] that the mean Nusselt number values obtained for low Ra (i.e. Ra ¼ 104 ) remain in good agreement, but small differences have been reported for high values of Ra (i.e. Ra ¼ 105 ) between the Bingham model without regularisation and bi-viscosity regularisation (i.e. a maximum of 5% difference in the mean Nusselt number has been reported for Bn ¼ 3, Ra ¼ 105 at Pr ¼ 1 between Refs. [18] and [19] in the Table 2 of Ref. [18]). However, this difference remains within the typical computational uncertainty and may also arise because of the differences in the mesh and numerical schemes.
2. Mathematical background
s ¼ lyield c_ for c_ sy =lyield
Above the critical shear rate (where jsj > sy ), qualitative and quantitative distributions of flow patterns remain independent of the value of lyield =l. This is because, as indicated by Eq. (2), the value of lyield determines the critical shear rate (i.e. c_ c ¼ sy =lyield )
ð1Þ ð2Þ
where c_ ij ¼ ð@ui =@xj þ @uj =@xi Þ are the components of the rate of strain tensor c_ , s is the stress tensor,sy is the yield stress, l is the plastic viscosity and lyield is the yield viscosity. Here, s and c_ are evaluated based on the second invariants of the stress and the rate of strain tensors in a pure shear flow respectively 0 " #1=2 " #1=2 ! @i:e: s ¼ 1=2 s : s _ _ _ , c ¼ 1=2c : c .
2.2. Computational domain The simulation configuration is schematically shown in Fig. 1, which demonstrates that the differentially heated horizontal walls are subjected to constant wall temperature boundary conditions. The bottom wall is taken to be at higher temperature than the top wall (i.e. T bottom > T top ) and all the other walls are considered to be adiabatic (i.e. the temperature gradient in the wall normal direction vanishes at the wall). No-slip and impermeability conditions are specified for all walls. The enclosure is taken to be cubic (i.e. H ¼ L ¼ W).
O’Donovan and Tanner [15] indicated that a value of lyield equal to 1000l mimics the true Bingham model in a satisfactory manner. Eqs. (1) and (2) transform ‘‘unyielded” region to a zone of high viscosity here in order to avoid singularities arising from a pure conduction solution (i.e. where c_ ¼ 0). As noticeable from Eqs. (1) and (2) the value of lyield represents the sensitivity of critical shear rate (i.e. c_ sy =lyield or c_ > sy =lyield ) here in order to mimic satisfactorily the Bingham model. This value affects only the size of the unyielded regions (regions where jsj sy [16]). The unyielded regions in Bingham fluid flows can be indicated by the locations of the Apparently Unyielded Regions (AURs). These regions are dependent on the prediction of jsj sy which is a function of the choice of lyield =l. Therefore, the grey regions should not be treated as exact unyielded zones and for this very reason Mitsoulis and Zisis [16] termed them as the Apparently Unyielded Regions.
Fig. 1. Schematic diagram of the simulation domain.
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
2.3. Governing equations The transient flow of an incompressible Bingham fluid is considered. The conservation equations for mass, momentum and energy for incompressible fluids under transient conditions take the following form:
@ui ¼0 @xi
@ui @ui 1 @P @ sij ¼ þ þ g b T T ref di2 þ uj @t @xj q @xi @xj
@T @T @2T ¼a þ uj @t @xj @xj @xj
ð3Þ ð4Þ
ð5Þ
where ui is the ith component of the velocity vector, T and P represent the temperature and pressure respectively, g; b are the acceleration due to gravity, the thermal expansion coefficient and a is the thermal diffusivity. The last term on the right-hand side of Eq. (4) originates due to Boussinesq’s approximation and the temperature difference between horizontal walls is considered to be small enough so that this approximation remains valid. Also, in Eq. (4), the Kronecker delta di2 indicates that buoyancy forces affect the flow only in the vertical direction (i.e. y direction). The reference temperature T ref is taken to be the cold wall temperature (i.e. T ref ¼ T C ). Eqs. (3)–(5) are valid for all fluid flow problems as long as continuum hypothesis remains valid, which is expected to be the case for most natural convection problems involving non-Newtonian fluids in engineering applications. The relation between s and c_ determines whether these equations are valid for Newtonian fluids or non-Newtonian fluids. The mass conservation and energy conservation equations (i.e. Eqs. (3) and (5)) do not change whether the fluid is Newtonian or non-Newtonian as long as the assumption of incompressible flow is valid, which is the case here. The same governing equations were used for the analysis of natural convection of Bingham fluid flows in several previous analyses (e.g. [5,7–14] and references therein) and the same approach has been adopted here. The viscous stress components sij are expressed with the help of Eqs. (1) and (2) (note that sy vanishes for Newtonian fluids). It is worth noting that in Eq. (5), the thermal diffusivity of the Bingham fluid cases (i.e. a) is kept independent of the shear rate. This can be also confirmed by the experimental studies [1,2] on Rayleigh-Bénard convection of yield stress fluids (e.g. Carbopol gel). It should be also noted that the shear rate of the Bingham fluids decreases immensely with increasing Bn for given value of shear stress. Therefore, change in the thermal diffusivity with the shear rate is expected to be negligible for the yield stress fluids. There have been several previous analyses [9–14] where the thermal diffusivity is considered to be independent of shear rate and the same approach has been adopted in this study. 2.4. The non-dimensional numbers According to Buckingham’s pi theorem, it is possible to show that the Nusselt number for natural convection of Bingham fluids in cubic enclosure can be expressed as: Nu ¼ f ðRa; Pr; Bn), where the nominal Rayleigh, Prandtl and Bingham numbers for Bingham fluids can be defined in the following manner:
Ra ¼ qgbðT H T C ÞH3 =ðlaÞ
ð6iÞ
Pr ¼ l=ðqaÞ
ð6iiÞ
pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Bn ¼ sy H=ðl gbðT H T C ÞHÞ
ð6iiiÞ
3
Here, the Bingham number represents the ratio of yield stress to viscous stress (i.e. non-dimensional yield stress). The nondimensional local heat transfer coefficient h is defined as:
hhw ¼ kð@T=@yÞy¼0 =ðT y¼0 T y¼H Þ
ð7iÞ
hcw ¼ kð@T=@yÞy¼H =ðT y¼H T y¼0 Þ
ð7iiÞ
where hw and cw represent the values at the hot wall and cold wall respectively, and k is the thermal conductivity. The mean heat transfer coefficient h and mean Nusselt number Nu at the walls are evaluated as:
h¼
1 LW
ZZ
hdxdz; Nu ¼ h H=k
ð8Þ
It is worth noting that hhw and hcw in Eq. (7) can be used in Eq. (8) in order to calculate mean Nusselt number at hot and cold walls (i.e. Nu hw and Nu cw ). 3. Numerical implementations 3.1. Numerical methodology To solve the non-linear set of governing equations in a finitevolume framework, the open-source CFD package OpenFOAM has been utilised. The pressure-velocity coupling has been addressed by the use of the Merged PISO-SIMPLE (PIMPLE) algorithm, which has been optimised for unsteady problems. Convective and diffusive fluxes are evaluated by second-order centered difference schemes. Temporal advancement has been achieved by the second-order Crank-Nicolson scheme with constant timestepping. It has been ensured that the Courant number is always sufficiently below unity so that the underlying physics is captured with sufficient temporal resolution. 3.2. Resolving the fine-scale structure It is vital to resolve all the relevant scales for conducting Direct Numerical Simulations (DNS) since the gradients of the turbulent fields cannot be accurately represented in the absence of sufficient resolution. It is well-known that the resolution of the finest scales in turbulent thermal transport is strongly dependent on the dissipation rates of turbulent kinetic energy and of the temperature variance. The volume averaged turbulent thermal dissipation rate can be expressed as hT i ¼ hað@T 0 =@xi Þ2 i where T 0 are the temperature fluctuations over the mean and the angled bracket suggests a volume averaging procedure. The turbulent kinetic energy dissipa 2 0 0 tion rate can be expressed as hi ¼ hm c_ij i where c_ij is the strain rate component based on the fluctuating velocity field and m is the kinematic viscosity. The mean kinetic energy dissipation rate hi is related to the mean Kolmogorov scale hgK i which is the smallest length scale of turbulence when m a (i.e. Pr 1). In case of m > a (i.e. Pr > 1), the smallest length scale is determined by the dissipation rate of scalar variance, and is represented by the Batchelor scale hgB i. These length scales are defined as:
hgK i ¼
m3=4 hgK i ; hgB i ¼ pffiffiffiffiffi Pr hi1=4
ð9Þ
Thus, DNS becomes excessively expensive for Pr 1, which is often encountered for liquid flows. Therefore, several alternative grid spacing criteria for the DNS of turbulent Rayleigh–Bénard convection have been proposed in the existing literature [20,21]. Shishkina et al. [21] proposed a grid spacing criterion for the DNS of turbulent Rayleigh–Bénard convection of Newtonian
4
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
(i.e. Bn ¼ 0) fluids. Accordingly, the permissible non-dimensional mesh size (Dx=H) can be expressed in terms of both hgK i and hgB i in the following manner:
Dx hgK i ¼ H H
Pr2 ðNu 1ÞRa
Dx hgB i ¼ H H
!1=4
1 Nu 1 Ra
for
Pr 1
ð10iÞ
for
Pr > 1
ð10iiÞ
!1=4
Shishkina et al. [21] expressed that the global measures of thermal transport, such as Nusselt number, time-averaged temperature profiles and volume-averaged dissipation rates, are mostly insensitive to moderate under-resolution of hgB i for turbulent Rayleigh–Bénard convection of Newtonian fluids. Dropkin and Somerscales [22] also experimentally analysed turbulent Rayleigh–Bénard convection of Newtonian fluids (i.e. Bn ¼ 0) in cubic enclosure for 5 104 Ra 7:17 108 and 0:02 Pr 11560. They proposed the mean Nusselt number correlation given by:
Nu ¼ 0:069 Ra 1=3 Pr 0:074
ð11iÞ
Similarly, Xia et al. [23] experimentally analysed turbulent Rayleigh–Bénard convection of Newtonian fluids in cubical enclosures for 2 107 Ra 3 1010 and 4 Pr 1350. They proposed the alternative mean Nusselt number correlation given by:
Nu ¼ 0:14 Ra 0:297 Pr 0:03
ð11iiÞ
In the current analysis, the criteria given in Eq. (10) by Shishkina et al. [21] have been used along with the Nu correlation of Refs. [22,23] for Ra ¼ 107 and Ra ¼ 108 , respectively for the purpose of estimating the grid resolution for both Newtonian fluid (i.e. Bn ¼ 0) and Bingham fluid cases. It is worth noting that the spatial resolution will be sufficient for Bingham fluid cases as long as it is sufficient for the Newtonian (i.e. Bn ¼ 0) fluid. The effective viscosity leff ¼ l þ sy =c_ in Bingham fluids increases with increasing Bingham number for a given value of plastic viscosity l and thus the effective Rayleigh number Raeff ¼ qgbðT H T C ÞH3 =leff a for Bingham fluid cases assumes much smaller values than the nominal value of Rayleigh numberRa. Thus, the effects of instability are much weaker in Bingham fluid cases with high values of Bingham number than in Newtonian (i.e. Bn ¼ 0) fluid cases with same nominal values of Ra and Pr. Based on these criteria, uniform Cartesian grids of 1503 and 2303 have been used for the analysis of Ra ¼ 107 and Ra ¼ 108 cases, respectively. It has been checked a-posteriori that the pffiffiffiffiffi chosen grid resolution was sufficient since yþ Pr (where ffiffiffiffiffiffiffiffiffiffiffi p yþ ¼ qus y=l with us ¼ sw =q, sw and y being the friction velocity, wall shear stress and wall normal distance of the wall adjacent grid point, respectively) remains smaller than unity for both mesh configurations in the case of Ra ¼ 10
7
and
Ra ¼ 108 : 3.3. Benchmarking and validation The present numerical scheme has been validated with DNS data of turbulent Rayleigh–Bénard convection of Newtonian (i.e. Bn ¼ 0) fluids in a cubic enclosure with respect to the results reported by Kaczorowski et al. [24] at Ra ¼ 107 , Ra ¼ 108 and Pr ¼ 4:38. Additionally, Vasco et al. [25] analysed laminar natural
convection of Bingham fluids in cubic enclosure with differentially heated sidewalls for 104 Ra 106 and Pr ¼ 1. The current numerical scheme has been validated with respect to the findings by Kaczorowski et al. [24] and Vasco et al. [25] for Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases, respectively, and an excellent agreement (less than 2% uncertainty on mean Nusselt number) has been obtained, as shown in Table 1. Additionally, Xu et al. [26] presented the mean Nusselt number results which are defined based on the volume and time average of D E @u @u the kinetic energy i:e: heiV ;t ¼ v @xij @xij and thermal dissiV ;t D E @T @T in addition to Reynolds pation rate i:e: heT iV;t ¼ a @x i @xi V ;t
number based on the root-mean-square (rms) velocity qffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi i:e: Rerms ¼ ð hux 2 þ uy 2 þ uz 2 iV ;t HÞ=v for the same configuration at Ra ¼ 107 and Pr ¼ 7. Accordingly, for additional validation, the comparison of the current study and Xu et al. [26] can be seen in Table 2, which shows that the current results agree very well with the data reported in Xu et al. [26]. 4. Results and discussion The effects of Rayleigh and Bingham number on the global flow behaviour and heat transfer characteristics are presented first. This is followed by a detailed discussion of the velocity and temperature spectra for all cases. The results presented in this section are obtained when the quasi-steady state has been reached. In order to demonstrate the quasi-steady state behavior, the temporal evolution of the mean Nusselt number Nu with normalised time is shown in Fig. 2 for Ra ¼ 107 and Pr ¼ 320. Under the quasi-steady state, the mean values of Nu for top and bottom walls become identical to each other, provided the oscillations due to the unsteady flow are neglected in this context. Zhang et al. [5] showed that Bingham fluid flows are unconditionally linearly stable for the quiescent initial condition. Thus, the quasi-steady Newtonian (i.e. Bn ¼ 0) fluid solutions for a given set of values of nominal Rayleigh and Prandtl numbers are used as the initial conditions for Bingham fluid simulations in the current analysis. Since it is possible to modify the yield stress by applying electrical or magnetic fields in magneto- and electro-rheological fluids, the quasi steady-state Newtonian solutions can be considered as realistic initial conditions with practical relevance for the Bingham fluid simulations. 4.1. Effects of variations of nominal Rayleigh number The three-dimensional views of non-dimensional isotherms (i.e. h ¼ ðT T C Þ=ðT H T C Þ) are shown in Fig. 3a for Ra ¼ 107 and 108 at Bn ¼ 0:43 and Pr ¼ 320. It should be noted that 0.05% by weight Table 1 Comparison of the mean Nusselt number for both Newtonian [24] and Bingham [25] fluid cases with the benchmark data. Newtonian (Bn ¼ 0) fluid case Kaczorowski et al. [24]
Present Study
7
Nu (Ra ¼ 10 , Pr ¼ 4:38)
16.1
16.3
Nu (Ra ¼ 108 , Pr ¼ 4:38)
32.2
32.8
Vasco et al. [25]
Present Study
8.946
8.816
Bingham (Bn ¼ 0:1) fluid case
6
Nu (Ra ¼ 10 , Pr ¼ 1:0)
5
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947 Table 2 Comparison of DNS results of the Reynolds and Nusselt numbers between Xu et al. [26] and present study at Ra ¼ 107 and Pr ¼ 7.
Xu et al. [26] Present
Nu Eq. (8)
Nue ¼
16.19 16.16
16.11 16.10
heiV;t H4 Pr2 Rav 3
þ1
Nue ¼
Rerms
heT iV ;t H2
aðDT Þ2
16.03 15.95
86.41 86.22
4.2. Effects of the variation of nominal Bingham number The distributions of h and W are shown in Fig. 5 for different Bn at Ra ¼ 108 and Pr ¼ 320. It is also worth noting that two addi-
Fig. 2.pVariation of Nu with non-dimensional simulation time (i.e. ffiffiffiffiffiffiffiffiffiffi t ¼ ðt a RaPrÞ=H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall at Ra ¼ 107 , Bn ¼ 0 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
Carbopol-water solution at Ra ¼ 107 and 108 yields Bn ¼ 0:43: The Bn value are obtained based on sy ¼ 0:007 Pa, q ¼ 961 kg=m3 ,
l ¼ 0:046 kg=ðmsÞ, a ¼ 1:5 107 m2 =s and b ¼ 2 104 1=K (which are taken from Table 1 of Ref. [2]) and different values of H and (T H T c ) according to Eq. (6). It can be seen from Fig. 3a that the iso-surfaces of h are more distorted for Ra ¼ 108 than in 7
Ra ¼ 10 , which is indicative of the strengthening of convection with increasing Ra due to the relative enhancement of buoyancy forces in comparison to the viscous resistance. This can be also confirmed by the path lines coloured by non-dimensional velocity pffiffiffiffiffiffiffiffi magnitude (i.e. W ¼ ui ui H=a) for Ra ¼ 107 and 108 , which are shown in Fig. 3b at Bn ¼ 0:43 and Pr ¼ 320. Fig. 3b indicates that the magnitude of the non-dimensional velocity increases with an increase inRa, which is indicative for the strengthening of convective transport. It also shows the Ra ¼ 108 flow field exhibits smaller flow structures compared to Ra ¼ 107 : The strengthening of convective thermal transport with increasing Ra can also be substantiated from the temporal evolutions of Nu for top and bottom walls in Fig. 4 for different Ra at Bn ¼ 0:43 and Pr ¼ 320. Fig. 4 shows that Nu increases dramatically with increasing Ra as a result of augmented convective thermal transport. The AURs for Ra ¼ 107 and 108 at Bn ¼ 0:43 and Pr ¼ 320 are shown in Fig. 3c (see the grey regions). It is evident from Fig. 3c that the size and the probability of finding AURs decrease with increasing Ra due to diminished viscous resistance in the flow. In the context of bi-viscosity regularisation, the flow does not stop in a true sense in AURs but the flow within the islands of AURs is too weak to influence the heat transfer rate and thus the exact shape and size of AURs do not affect the heat transfer rate and the mean Nusselt number.
tional Bingham numbers (i.e. Bn ¼ 0:215 and 2.15) at Ra ¼ 108 are chosen for the purpose of a numerical experiment. For the sake of brevity the results for Bn ¼ 0:215 are not always shown explicitly due to the monotonic effect of Bn. It can be seen from Fig. 5a that the iso-surfaces become less deformed and thermal plumes from the bottom wall become less apparent as Bn increases. These are the indications of the relative weakening of buoyancy force due to the augmented flow resistance as a result of increased yield stress. This also leads to the weakening of convection with increasingBn, which can be substantiated from Fig. 5b where the path lines colored by non-dimensional velocity magnitude for the corresponding cases are shown. Fig. 5b shows that the non-dimensional velocity magnitude decreases with increasing Bn due to the weakening of the convective thermal transport. The temporal evolutions of Nu for top and bottom walls are shown in Fig. 6 for three different Bn values at Ra ¼ 108 and Pr ¼ 320. All mean Nusselt numbers obtained for different values of Bn for Ra ¼ 107 and Ra ¼ 108 are summarised in Table 3. Eq. (11) suggests a scaling of mean Nusselt number as Nu Ran Pr m but this formula does not take into account the variable fluid properties of Bingham fluids. Introducing the volume and time averaged effective viscosity hleff iV;t , effective Rayleigh and Prandtl numbers can be defined in the following manner:
Raeff ¼ RaBn¼0
l hleff iV;t
;
Preff ¼ PrBn¼0
hleff iV;t
l
ð12Þ
Consequently, it can be expected that the Nusselt number for Bingham fluids NuBn possibly can be scaled with the nominal (Newtonian) Nusselt number NuBn¼0 and the effective viscosity as:
NuBn NuBn¼0
!p
l hleff iV;t
:
ð13Þ
Fig. 7 demonstrates that Eq. (13) approximates the mean Nusselt number for Bingham fluids reasonably well using the concept of effective viscosity and an exponent of p ¼ 0:13: Here, the effective viscosities have been determined after conducting the simulations. The ratio of l=hleff iV;t is expected to be unity for Newtonian fluids (Bn ¼ 0) and to approach small values for larger Bingham numbers when the fluid is unyielded. An analysis of the existing data shows that the ratio of l=hleff iV;t can be reasonably approximated by the following expression:
l=hleff iV;t ¼
1 min ð1 þ b Bnm ; 104 Þ
ð14Þ
with m 2:5 and b e6 ðb e5 Þ for Ra ¼ 107 (Ra ¼ 108 ) where e is Euler’s number. However, more data will be needed to develop a generally valid expression. The far right circle (square) in Fig. 7 represents the Newtonian case for Ra ¼ 107 (Ra ¼ 108 ) and the remaining Nusselt numbers decrease with increasing Bn. The weakening of convective transport, as reflected in the reduction of the mean Nusselt number Nu, with increasing Bn is
6
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Fig. 3. (a) Distributions of iso-surfaces of non-dimensional temperature h, (b) The path lines are coloured by non-dimensional velocity magnitude (i.e. W), (c) Distribution of AURs for different values of Ra at Bn ¼ 0:43 and Pr ¼ 320. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
(see Fig. 6). For very large values of Bingham number, the boundary layer thickness becomes of the order of the enclosure size H. At that stage, the fluid flow does not influence the thermal transport and heat transfer takes place principally due to conduction, which is reflected in the unity mean Nusselt number. The distributions of magnitude of normalised shear stress (i.e. pffiffiffiffiffiffiffiffiffiffi s ¼ sij sij H2 =ðqa2 Þ) values on the x z plane at y ¼ H are shown
Fig. 4. Temporal variation of Nu with non-dimensional simulation time (i.e. pffiffiffiffiffiffiffiffiffiffi t ¼ t a RaPr =H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall for different Ra at Bn ¼ 0:43 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
consistent with the distributions of AURs shown in Fig. 8a for different values of Bnat Ra ¼ 108 and Pr ¼ 320. It can be seen from Fig. 8a that the size and the probability of finding AURs increase with increasing Bn due to the weakening of convective transport, and this gives rise to a reduction in the mean Nusselt number
in Fig. 8b for different Bn values at Ra ¼ 108 and Pr ¼ 320. The Rayleigh-Bénard cells can be seen for the Newtonian (i.e. Bn ¼ 0) case and also for small values of Bn (e.g. Bn ¼ 0:43). However, the Rayleigh-Bénard cells do not clearly appear for large values of Bn (i.e. Bn ¼ 2:15) due to the diminished strength of the thermal convection. 4.3. Velocity and temperature spectra It can clearly be seen from Figs. 5 and 8 that the flow structures are strongly affected by the Bingham number. Consequently, it will be of interest to characterise the range of flow structures in terms of the energy spectra of temperature and velocity fluctuations. In this context, it is important to mention that determination of universal scaling laws for inertial range two-point statistics, as done
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
7
Fig. 5. (a) Distributions of iso-surfaces of non-dimensional temperature h, (b) The path lines are coloured by non-dimensional velocity magnitude (i.e. W) for different values of Bn at Ra ¼ 108 and Pr ¼ 320. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
pffiffiffiffiffiffiffiffiffiffi Fig. 6. Temporal variation of Nu with non-dimensional simulation time (i.e. t ¼ ta RaPr =H2 ) for hot (i.e. Nuhw ) and cold (i.e. Nucw ) wall for different Bn at Ra ¼ 108 and Pr ¼ 320. The time averaged mean Nusselt number is shown with a blue dashed line. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
by Kolmogorov in ideal hydrodynamics turbulence, is a long-standing challenge in Rayleigh-Bénard turbulence, even in the context of Newtonian flows [27] and the same holds true for scalar spectra of high Prandtl (Schmidt) number fluids.
In the inertial range, Bolgiano [28] and Obukhov [29] proposed a dual cascade: for j small wave numbers in the range LB < j1 (where LB is the Bolgiano scale) they predicted a dominance of the buoyancy force over the inertial force leading to velocity and
8
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Table 3 Summary of the mean Nusselt number obtained in the current analysis for both Newtonian (i.e. Bn ¼ 0) and Bingham fluid cases at Pr ¼ 320.
Nu
Ra ¼ 107
Ra ¼ 108
Bn ¼ 0 Bn ¼ 0:215 Bn ¼ 0:43 Bn ¼ 2:15
17.43 14.25 11.095 1.0
32.8 27.79 23.35 10.20
Fig. 7. Mean Nusselt number based on the DNS results versus the Nusselt number according to Eq. (13) using an exponent of p ¼ 0:13 for Ra ¼ 107 and 108 and different Bingham numbers.
velocity [30,38] and hence the applicability of this hypothesis is questionable in this configuration. In the light of the above discussion and the remaining uncertainties in regards of universal scaling laws, it will be interesting to analyse velocity and temperature spectra for the flow configurations considered in this work. The present configuration does not have homogeneous directions and hence, despite the aforementioned difficulties, the temporal spectra will be presented. It has been decided not to transform the spectra into the equivalent spatial domain. Nevertheless, a suitable mean velocity scale is needed for determining the Nyquist frequency in the Power Spectral Densities (PSD), which are subsequently discussed. The WienerKhinchin theorem states that the PSD (u0 ) is identical to the fluctuating kinetic energy spectrum in case of homogeneous flows. According to the Nyquist theorem, the maximum detectable frequency in a signal is given by f max ¼ 1=ð2DtÞ where Dt is the sampling interval, which corresponds in the present case to the constant time step in the simulation. In Computational Fluid Dynamics (CFD), Dt typically has an upper bound determined by the Courant number. However, the time step can be arbitrarily small which would result in very high resolvable frequencies. It is intuitively clear that for a given spatial resolution the resolved frequency range cannot be arbitrary large. All simulations have been run with a Courant number smaller than unity. Accordingly, we determine the Nyquist frequency using Taylors hypothesis as f max ¼ 1=ð2Dx=U con Þ where U con is a suitable convection velocity. pffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffiffi Table 4 shows the mean velocity (i.e. U mean ¼ ux 2 þ uy 2 þ uz 2 ), p ffiffiffiffiffiffiffiffiffiffiffi 0 the fluctuation velocity scale (i.e. u ¼ 2=3k) where k is the h i 2 turbulent kinetic energy (i.e. k ¼ 0:5 ðux 0 Þ2 þ uy 0 þ ðuz 0 Þ2 ) and 0
temperature spectra as j11=5 and j7=5 , respectively. For intermediate wave numbers, the temperature field evolves as a passive scalar, and both the velocity and temperature spectra depict a j5=3 range. As reviewed in Boffetta and Mazzino [27], the Bolgiano scale LB coincides with the integral scales L in three dimensions. Further, according to Lohse and Xia [30], there is a fundamental lack of clear scale separation between LB and L which makes the identification of such a regime difficult in practice. In the high wavenumber range, for high Prandtl numbers, there are two dissipation scales, which are the Kolmogorov scale hgK i and the Bachelor scale hgB i [31], respectively. For high Prandtl numbers one obtains hgB i < hgK i and Batchelor [32] predicted that the scalar spectrum should display a j1 dependence in the so called viscous-convective subrange hgB i < j1 < hgK i. In order to better distinguish the spectra for high and low Prandtl numbers, reference simulations with Pr ¼ 1 will be presented as well in this section. According to the review of Warhaft [33], the experimental evidence has been elusive and interested readers are referred to references in Ref. [33] for further information. A detailed discussion has been provided as well by Antonia and Orlandi [34]. In this context, it is worth mentioning that a modified scaling of the form j1p with p Re1=2 has been suggested by Kerstein [35] and that the location of the j1 range has sometimes been observed for j1 > hgKi (Ref. [36]), which is inconsistent with the classical expectation. The lack of resolution especially in early measurements could be an issue, which might explain why the j1 scaling is not well-established [36]. Another problem is that typically time domain results are compared to theoretical predictions for spatial domains by invoking Taylor’s frozen-flow hypothesis [37], which requires the velocity fluctuations to be much smaller than the mean flow velocity. However, it is well-known that the rootmean-square (i.e. rms) velocity fluctuation in Rayleigh-Bénard convection is often of the same order of magnitude as that of the mean
fluctuation intensity (u =U mean ) in the centre of the domain (i.e. x=H ¼ y=H ¼ z=H ¼ 0:5). In order to compare the spectra, two additional simulations have been performed for Newtonian (i.e. Bn ¼ 0) fluids at Pr ¼ 1:0 for Ra ¼ 107 and 108 . It can be seen from Table 4 0 that the ratio u =U mean remains of the order of 75 3% for the 0 Pr ¼ 320 case, whereas this ratio is u =U mean 63% for Pr ¼ 1. This demonstrates the well-known fact that even in the centre of the domain the applicability of Taylor’s hypothesis remains questionable. Hence, in order to determine the Nyquist frequency, U con ¼ U mean u0 has been used. From the foregoing, it becomes clear that the cut-off frequency is not a precise number but rather can be attributed to a certain range. The same holds true for the frequencies corresponding to Kolmogorov and Bachelor scales, which will henceforth be denoted as fhgKi = Ucon/(2hgKi) and fhgBi = Ucon/(2hgBi) and have been estimated using Eqs. (11i) and (11ii), respectively. Fig. 9 shows the PSD of temperature and velocity fluctuation 0
0
(i.e. T and u respectively) for different Pr values at Ra ¼ 107 in the Newtonian (i.e. Bn ¼ 0) case in the centre of the domain based on single probing point. According to Bachelor [28], one would expect a f
1
(i.e.
j1 ) scaling in the range f hgK i < f < f hgB i
1
(i.e. hgB i < j < hgK i) for the temperature spectra for Pr ¼ 320, whereas the velocity spectra should decay much faster in the dissipation range for f > f hgK i . The estimated lower and upper limits of f hgK i are shown with dashed and solid blue lines in Fig. 9, whereas the corresponding Bachelor scales are shown in black. It is important to note that the cases with high Pr and Bn > 0 do not feature a fully developed turbulent flow field. Hence the scaling laws are shown for the purpose of reference. The PSD (T 0 ) for Pr ¼ 320 depicts a decay close to f
1
for f f hgK i but this range does not
reach up to f hgB i . An inertial range scaling for f < f hgK i is not clearly visible. In contrast for Pr ¼ 1, there is a small inertial range for f < f hgK i , whereas the PSD (T 0 ) decays much faster for f > f hgK i : As
9
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Fig. 8. (a) Distribution of three-dimensional AURs, (b) distribution of the magnitude of normalised shear stress (s ¼ s H2 =ðqa2 Þ) on the cold top wall (i.e. x z plane view at y ¼ H) for different Bn at Ra ¼ 108 and Pr ¼ 320.
Table 4 0 0 Summary of the mean velocity (i.e. U mean ), fluctuation velocity scale (i.e. u ) fluctuation intensity (u =U mean ) in the centre of the domain (i.e. x=H ¼ y=H ¼ z=H ¼ 0:5). Cases
0
U mean
u ¼ 3
2.0397 10
Ra ¼ 107 ; Pr ¼ 1:0; Bn ¼ 0
pffiffiffiffiffiffiffiffiffiffiffi 2=3k
0
u =U mean 3
0.634
3
1.2939 10
3
7
5.5003 10
4.1529 10
0.755
8
3.1556 103
1.9778 103
0.627
8
3
8.9840 10
3
6.3900 10
0.711
8
5.7679 103
4.4898 103
0.778
8
3
3
0.722
Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 1:0; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0 Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0:215
4.9757 10
Ra ¼ 10 ; Pr ¼ 320; Bn ¼ 0:43
expected, the velocity spectra for both Prandtl numbers depict a dissipation range behaviour for scales smaller than the Kolmogorov scales. However, the decay is considerably slower for Pr ¼ 320 which might be a consequence of the slower decay of 0
the temperature fluctuations (i.e. T ) and the active scalar nature of this quantity (back-coupling in the momentum equation via the buoyancy term). The PSD (T 0 ) and PSD (u0 ) look qualitatively similar for Ra ¼ 108 , which is shown in Fig. 10 for Pr ¼ 1 and Pr ¼ 320 in the centre of the domain based on single probing point. Again the PSD (u0 ) decays slower for Pr ¼ 320; possibly due to the coupling between momentum and temperature equation. Finally, it is worth nothing that the features of a viscous-convective range are slightly more pronounced for the Ra ¼ 108 case with Pr ¼ 320. As thermal convection weakens with decreasing (increasing) Ra (Bn), fluid flow starts to exhibit a nearly laminar behaviour for small (high) Ra (Bn) values. This can be easily seen from the
3.5938 10
iso-surfaces of non-dimensional temperature in Fig. 5. Therefore, the discussion of spectra of the Bingham fluids has been limited to the case with high Ra and small Bn values (i.e. Bn ¼ 0:215 and 0.43 for Ra ¼ 108 ). Accordingly, the PSD of rms velocity fluctuation 0
(i.e. u ) and temperature (i.e. T 0 ) for different Bn values at Ra ¼ 108 and Pr ¼ 320 in the centre of the domain based on single probing point are shown in Fig. 11. It should be noted here that the spatial cut-off frequencies based on mean Kolmogorov f hgK i and Batchelor f hgB i length scales have been calculated based on the criterion in Eq. (10) given by Shishkina et al. [21] but using the concept of effective Rayleigh and Prandtl numbers introduced in Eq. (12). As effective viscosity is expected to be higher for Bn > 0, this results in a smaller (higher) effective Rayleigh (Prandtl) numbers. For instance, hleff iV ;t for the Bingham fluid cases (i.e. Bn ¼ 0:215 and Bn ¼ 0:43) was obtained around 3.5 and 27 times bigger than l in the Newtonian case (i.e. Bn ¼ 0) respectively. This implies that
10
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
0
Fig. 9. Power spectral density of (a) temperature (i.e. T ), (b) rms velocity fluctuations (i.e. u ) for Ra ¼ 107 for Newtonian (i.e. Bn ¼ 0) fluid cases in the centre of the domain based on single probing point for different Pr values. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and Batchelor fhgBi 0 0 length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u and U mean u , respectively.
effective values of Rayleigh (Prandtl) numbers based on hleff iV;t values are 3.5 and 27 times smaller (bigger) for the Bingham fluid cases (i.e. Bn ¼ 0:215 and Bn ¼ 0:43) than their respective nominal values. Accordingly, this yields larger Kolmogorov and Batchelor scales or alternatively, smaller values of f hgK i and f hgB i . Thus, it can be observed from Fig. 11 that for the Bingham fluid case similar observations hold true as for the Newtonian fluid (i.e. Bn ¼ 0) with Pr ¼ 320. In particular, there is a small range where viscousconvective behaviour can be observed in the PSD (T 0 ). In order to facilitate easier observation of the viscousconvective range like scaling, Fig. 12 shows the power spectral density of temperature fluctuations multiplied by frequency 0 (i.e. PSD T f ) for two different Ra at Pr ¼ 1 (black) and Pr ¼ 320 (blue) for the Newtonian (i.e. Bn ¼ 0) fluid case as well as for Ra ¼ 108 , Pr ¼ 320, Bn ¼ 0:215 (green) and Ra ¼ 108 , Pr ¼ 320, Bn ¼ 0:43 (red) for Bingham fluid cases in the centre of the domain based on single probing point. Fig. 12 clearly shows that despite having the same nominal Rayleigh and Prandtl num-
bers, the Bingham fluid cases (red and green line) exhibit a different PSD with more energy on large scales compared to the Newtonian fluid (blue line) due to the smaller effective Ra. Further, it can be observed that temperature fluctuation spectra multiplied by frequency decay considerably faster than f
1
for f f hgK i in the
Newtonian (i.e. Bn ¼ 0) fluids (for both Ra ¼ 107 and Ra ¼ 108 ) at Pr ¼ 1, whereas a small plateau tends to appear for Newtonian (Bn ¼ 0 for Ra ¼ 108 ) and Bingham (Ra ¼ 108 , Bn ¼ 0:215 and Bn ¼ 0:43) fluid cases at Pr ¼ 320 at for f f hgK i . This viscousconvective range is most pronounced for the Newtonian case with the highest Rayleigh number for Pr ¼ 320. 5. Conclusions Direct Numerical Simulations (DNS) have been conducted to analyse three-dimensional Rayleigh-Bénard convection (differentially heated horizontal walls with heated bottom wall, i.e. heated from below and cooled from above, whereas the other walls are
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
11
0
Fig. 10. Power spectral density of (a) temperature (i.e. T ), (b) rms velocity fluctuations (i.e. u ) for Ra ¼ 108 for Newtonian (i.e. Bn ¼ 0) fluid cases in the centre of the domain based on single probing point for different Pr values. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and Batchelor fhgBi 0 0 length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u and U mean u , respectively.
adiabatic) of yield stress fluids obeying a Bingham model in cubic domains for high values of nominal Rayleigh number (107 and 108 Þ for a representative value of nominal Prandtl number (i.e. Pr = 320 which corresponds to a 0.05% by weight Carbopol-water solution). It has been found that the mean Nusselt number Nu increases with increasing values of Rayleigh number for both Newtonian and Bingham fluids. However, weaker convective transport in Bingham fluid cases leads to smaller values of Nu than that obtained in the case of Newtonian fluids with the same nominal value of Rayleigh numberRa and a scaling based on volume averaged effective viscosities has been suggested. The flow domain becomes increasingly unyielded and the effects of Rayleigh-Bénard instability become less pronounced as the Bingham number increases. For Bingham fluids, thermal convection decreases considerably with increasing Bingham number. Nevertheless, a small viscous1
convective type range with a f scaling has been observed for the temperature fluctuation spectra, roughly located close to the
Kolmogorov scale, in all cases for Pr ¼ 320. It has been observed that for identical nominal Rayleigh numbers, the Bingham fluid cases exhibit a smaller range of scales with more energy towards the large scales when compared to the corresponding Newtonian fluid case. By comparing simulations for unity Prandtl number and Pr ¼ 320, it has been found that the decay of turbulent kinetic energy in the dissipation range is considerably slower for Pr ¼ 320, which has been attributed to the coupling between temperature and momentum equation combined with the slower decay of thermal variance. To the best of the authors’ knowledge, this is the first analysis of Rayleigh-Bénard convection of yield stress fluids at high nominal Rayleigh numbers. Future investigations are aimed at the analysis of similar configurations at even higher values of Ra: Establishing correlations for the mean Nusselt number for complex nonNewtonian fluids will be of great importance for a variety of technical applications.
12
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
0
0
Fig. 11. Power spectral density of rms velocity fluctuations (i.e. u ) on the top and temperature fluctuations (i.e. T ) on the bottom, for different Bn values at Ra ¼ 108 and Pr ¼ 320 in the centre of the domain based on single probing point. The spatial cut-off frequencies are shown with blue and black lines based on mean Kolmogorov fhgKi and 0 Batchelor fhgBi length scales, respectively and indicate a lower and upper limit depending on the velocity scale used for converting them into frequency space, i.e. U mean þ u 0 and U mean u , respectively.
0
Fig. 12. Power spectral density of temperature multiplied by frequency (i.e. PSD ðT Þ f ) for different Ra at Pr ¼ 1 (black), Pr ¼ 320 (blue) for Newtonian (i.e. Bn ¼ 0) fluid case andRa ¼ 108 , Pr ¼ 320 and Bn ¼ 0:215 (green) and Ra ¼ 108 , Pr ¼ 320 and Bn ¼ 0:43 (red) for Bingham fluid cases in the centre of the domain based on single probing point. The mean Kolmogorov fhgKi length scale for all cases is shown by the vertical lines with the respective colours. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)
S. Yigit et al. / International Journal of Heat and Mass Transfer 147 (2020) 118947
Declaration of Competing Interest The authors declared that there is no conflict of interest. Acknowledgements The authors gratefully acknowledge that N8HPC (North East UK national High-Performance Computing) and Rocket HighPerformance Computing (in-house cluster of Newcastle University) have been used for the current analysis. Computer resources for this project have also been provided by the Gauss Centre for Supercomputing/Leibniz Supercomputing Centre under grant no. pn56di. References [1] M. Darbouli, C. Métivier, J.M. Piau, A. Magnin, A. Abdelali, Rayleigh- Bénard convection for viscoplastic fluids, Phys. Fluids 25 (2013) 023101. [2] Z. Kebiche, C. Castelain, T. Burghelea, Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluid, J. Nonnewton. Fluid Mech. 203 (2014) 9–23. [3] J.E. Stangroom, Electrorheological fluids, Phys. Technol. 14 (1983) 290–296. [4] N. Wereley, Magnetorheology: Advances and Applications, Royal Society of Chemistry, 2013. [5] J. Zhang, D. Vola, I.A. Frigaard, Yield stress effects on Rayleigh-Bénard convection, J. Fluid Mech. 566 (2006) 389–419. [6] N.J. Balmforth, A.C. Rust, Weakly nonlinear viscoplastic convection, J. Nonnewton. Fluid Mech. 158 (2009) 36–45. [7] A. Vikhansky, On the onset of natural convection of Bingham liquid in rectangular enclosures, J. Nonnewton. Fluid Mech. 165 (2010) 1713–1716. [8] O. Turan, S. Yigit, N. Chakraborty, Critical condition for Rayleigh-Bénard convection of Bingham fluids in rectangular enclosures, Int. Commun. Heat Mass Transf. 86 (2017) 117–125. [9] O. Turan, N. Chakraborty, R.J. Poole, Laminar Rayleigh-Bénard convection of yield stress fluids in a square enclosure, J. of Nonnewton. Fluid Mech. 171 (2012) 83–96. [10] O. Turan, R.J. Poole, N. Chakraborty, Boundary condition effects on natural convection of Bingham fluids in a square enclosure with differentially heated horizontal walls, J. Comput. Therm. Sci. 4 (2012) 77–97. [11] S. Yigit, R.J. Poole, N. Chakraborty, Effects of aspect ratio on natural convection of Bingham fluids in rectangular enclosures with differentially heated horizontal walls heated from below, Int. J. Heat Mass Transf. 80 (2015) 727– 736. [12] S. Yigit, R.J. Poole, N. Chakraborty, Laminar natural convection of Bingham fluids in inclined differentially heated square enclosures subjected to uniform wall temperatures, J. Heat Transf. 137 (2015) 052504. [13] S. Yigit, S. Chen, P. Quinn, N. Chakraborty, Numerical investigation of laminar Rayleigh-Bénard convection of Bingham fluids in square cross-sectioned cylindrical enclosures, Int. J. Therm. Sci. 110 (2016) 356–368. [14] S. Yigit, N. Chakraborty, Influences of aspect ratio and wall boundary condition on laminar Rayleigh-Bénard convection of Bingham fluids in rectangular enclosures, Int. J. Numer. Meth. Heat Fluid Flow 27 (2017) 310–333.
13
[15] E. O’Donovan, R. Tanner, Numerical study of the Bingham squeeze film problem, J. Nonnewton. Fluid Mech. 15 (1984) 75–83. [16] E. Mitsoulis, T. Zisis, Flow of Bingham plastics in a lid-driven square cavity, J. Nonnewton. Fluid Mech. 101 (2001) 173–180. [17] I. Karimfazli, I.A. Frigaard, A. Wachs, A novel heat transfer switch using the yield stress, J. Fluid Mech. 783 (2015) 526–566. [18] R.R. Huilgol, G.H.R. Kefayati, Natural convection on Bingham fluid using the operator splitting method, J. Nonnewton. Fluid Mech. 220 (2015) 22–32. [19] O. Turan, N. Chakraborty, R.J. Poole, Laminar natural convection of Bingham fluids in a square enclosure with differentially heated side walls, J. Nonnewton. Fluid Mech. 165 (2010) 901–913. [20] G. Grötzbach, Spatial resolution requirements for direct numerical simulation of the Rayleigh-Bénard convection, J. Comp. Phys. 49 (1983) 241–264. [21] O. Shishkina, R.J.A.M. Stevens, S. Grossmann, D. Lohse, Boundary layer structure in turbulent thermal convection and its consequences for the required numerical resolution, New J. Phys. 12 (2010) 075022. [22] D. Dropkin, E. Somerscales, Heat transfer by natural convection in liquids confined by two parallel plates which are inclined at various angles with respect to the horizontal, J. Heat Transf. 87 (1965) 1–6. [23] K.-Q. Xia, S. Lam, S-Q, Zhou, Heat-flux measurement in high-Prandtl-number turbulent Rayleigh-Bénard convection, Phys. Rew. Lett. 88 (2002) 064501. [24] M. Kaczorowski, C. Kai-Leong, X. Ke-Quing, Turbulent flow in the bulk of Rayleigh-Bénard convection: aspect-ratio dependence of the small-scale properties, J. Fluid Mech. 747 (2014) 73–102. [25] D.A. Vasco, N.O. Moraga, G. Haase, Parallel finite volume method simulation of three dimensional fluid flow and convective heat transfer for viscoplastic nonNewtonian fluids, Num. Heat Transf. Part A 66 (2014) 990–1019. [26] Xu. Ao, Le Shi, Heng-Dong Xi, Lattice Boltzmann simulations of threedimensional thermal convective flows at high Rayleigh number, Int. J. Heat Mass Transf. 140 (2019) 359–370. [27] G. Boffetta, A. Mazzino, Incompressible Rayleigh-Taylor turbulence, Annu. Rev. Fluid Mech. 49 (2017) 119–143. [28] R. Bolgiano, Turbulent spectra in a stably stratified atmosphere, J. Geophys. Res. 64 (1959) 2226–2229. [29] A.M. Obukhov, On the influence of Archimedean forces on the structure of the temperature field in a turbulent flow, Dokl. Akad. Nauk. SSR 125 (1959) 1246– 1248. [30] D. Lohse, K.Q. Xia, Small-scale properties of turbulent Rayleigh-Bénard convection, Annu. Rev. Fluid Mech. 42 (2010) 335–364. [31] H. Tennekes, J.L. Lumley, A First Course in Turbulence, MIT Press, Cambridge, 1972. [32] G.K. Batchelor, Small-scale variation of convected quantities like temperature in a turbulent field. Part 1. General discussion and the case of small conductivity, J. Fluid Mech. 5 (1959) 113–133. [33] Z. Warhaft, Passive scalars in turbulent flows, Annu. Rev. Fluid Mech. 32 (2010) 203–240. [34] R.A. Antonia, P. Orlandi, Effect of Schmidt number on small-scale passive scalar turbulence, Appl Mech Rev. 56 (2003) 615–632. [35] A.R. Kerstein, Scaling properties of the viscous-convective scalar spectral subrange in turbulent jets, Phys. Fluids A 3 (1991) 1832–1834. [36] D. Bogucki, J.A. Domaradzki, P.K. Yeung, Direct numerical simulations of passive scalars with Pr>1 advected by turbulent flow, J. Fluid Mech. 343 (1997) 111–130. [37] G.I. Taylor, The spectrum of turbulence, Proc. R. Soc. Lond. A 164 (1938) 476– 490. [38] Q. Zhou, K. Xia, Comparative experimental study of local mixing of active and passive scalars in turbulent thermal convection, Phys. Rev. 77 (2008) 056312.