Large-eddy simulation of natural convection in an asymmetrically-heated vertical parallel-plate channel: Assessment of subgrid-scale models

Computers & Fluids 59 (2012) 101–116

Contents lists available at SciVerse ScienceDirect

Computers & Fluids j o u r n a l h o m e p a g e : w w w . e l s e v i e r . c o m / l o c a t e / c o m p fl u i d

Large-eddy simulation of natural convection in an asymmetrically-heated vertical parallel-plate channel: Assessment of subgrid-scale models G.E. Lau a, G.H. Yeoh a,b, V. Timchenko a,⇑, J.A. Reizes a a b

School of Mechanical and Manufacturing Engineering, University of New South Wales, Sydney, NSW 2052, Australia Australian Nuclear Science and Technology Organisation (ANSTO), PMB 1, Menai, NSW 2234, Australia

a r t i c l e

i n f o

Article history: Received 6 December 2010 Received in revised form 29 November 2011 Accepted 16 January 2012 Available online 2 February 2012 Keywords: Natural convection Large-eddy simulation Subgrid-scale models

a b s t r a c t The performance of four different large-eddy simulation subgrid-scale models has been examined a posteriori for natural convection in an asymmetrically-heated vertical parallel-plate channel with a high aspect ratio. The compressible three-dimensional Favre-filtered mass, momentum and energy conservation equations have been closed using the Smagorinsky, dynamic, approximate localised dynamic and Vreman models. A two-stage predictor–corrector numerical methodology for low-Mach-number compressible flows was adopted to strongly couple the density with the Navier–Stokes equations. Based on the comparison with experimental data, it has been shown that the Smagorinsky model predicts inaccurate near-wall flow dynamics and delayed transitional behaviour while both dynamic procedures to compute the Smagorinsky model coefficient result in over prediction of wall temperatures, suggesting an under estimation of subgrid-scale dissipation. The time extrapolation procedure utilised in the approximate localised dynamic model has been shown to produce better adaptation towards the local flow behaviour when compared with the standard dynamic model. At the same time, time-averaged wall temperature and velocity field profiles have been well captured by the Vreman model, demonstrating its superiority when compared to the rest of the models. Ó 2012 Elsevier Ltd. All rights reserved.

1. Introduction Success of large-eddy simulation (LES) in predicting turbulent flows has been demonstrated in its vast employment in diverse engineering applications. For example, the application of LES in studying turbulent buoyant jets [1–6], rotating flows [7,8], vortex breakdown and recirculation flows [9,10], and flow over an obstruction [11,12] have shown reasonably accurate prediction of turbulent flows. Furthermore, success of LES has recently been extended to accurately simulate complex turbulent flows such as shock waves and supersonic flows [13–17], bubbly flows [18], binary particle systems [19], and flows in human lungs and intra-thoracic airways [20,21]. Despite its application to various types of turbulent flows, the performance of various subgrid-scale (SGS) models has not been thoroughly investigated for natural convective flows. For convection-dominated flows with high Reynolds numbers, the inertia force is much larger than the viscous force, leading to a reasonable assumption of turbulent homogeneity; the statistical properties of turbulence do not vary spatially. On the other hand, natural convective flows represent flows where the viscous forces may be comparable in magnitude to the buoyancy force especially when ⇑ Corresponding author. Tel.: +61 2 9385 4148; fax: +61 2 9663 1222. E-mail address: [email protected] (V. Timchenko). 0045-7930/$ - see front matter Ó 2012 Elsevier Ltd. All rights reserved. doi:10.1016/j.compfluid.2012.01.006

the Rayleigh number is moderately high and the flow is undergoing transition and turbulence. This suggests that most natural convective flows are heterogeneous and anisotropic. In addition, the anisotropic behaviour of most buoyancy-driven flows is also partly attributed to the small time- and length scales typically found in such flows. This results in a rather tough implementation of LES because most SGS models are developed based on the homogeneity framework. Hence, the question beckons in the capability of LES in aptly predicting natural convective flows where both buoyancy production and viscous dissipation are significant. In addition, the non-linear interaction between the viscous, pressure and buoyancy forces in such flows depicts even more vigorous evaluation for the SGS models. For inhomogeneous and complex turbulent flows, it is essential that the SGS model predicts accurate SGS dissipation to correctly represent the physical flow behaviour by adapting to the local flow dynamics. Regrettably, this prerequisite has not been satisfied by the well-known Smagorinsky model [22] even though its implementation in practical LES is robust. The Smagorinsky model has some prominent shortcomings: (i) it requires an a priori model coefficient CS which is flow dependent; (ii) it predicts incorrect asymptotic flow behaviour near a wall; and (iii) it predicts inaccurate eddy viscosity in a laminar or transitional flow. Extensive effort has been made to improve on the Smagorinsky model by developing various damping functions to account for the near-wall

102

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

effects [23,24]. These functions are of limited universality as they presume a particular form of the flow dynamics in the region considered. In cases such as buoyant thermal plumes and freestanding fires where the flow is not bounded by the presence of a wall, the Smagorinsky model has been successfully applied [25–28]. However, Yan and Nilsson [29] who studied natural convection adjacent to a vertical flat plate by using the Smagorinsky model have found local discrepancies in the profile of the maximum normalised temperature fluctuation particularly at locations downstream of the vertical plate where the flow exhibited complex transitional dynamics. This finding highlights the inaccuracy of the Smagorinsky model in predicting the weak transitional boundary layer flow particularly for natural convection. Natural convection in either horizontal or vertical cylinders has also attracted considerable attention by various researchers. Miki et al. [30] who studied turbulent natural convection in a concentric horizontal annuli by using the Smagorinsky model concluded that at high wave numbers, the amount of SGS dissipation depended strongly on CS. The same model which has also been used by Barhaghi et al. [31] in an a priori study on natural convection boundary layer in a vertical cylinder has resulted in discrepancies particularly in the turbulent shear and normal stresses as well as temperature fluctuations profiles when compared against experimental data. It is also worth noting that, in a priori studies such as that of Barhaghi et al. [31], the effects of the modelling errors are neglected because a comparison is made between the exact subgrid stresses with those predicted by a subgrid model evaluated on the basis of the filtered exact solution, thus resulting in validations which are only relative in value [32]. Other similar physical space static models which have also been developed include the structure function model [33] and the mixed scale model [34] but they have not been applied extensively in natural convective flows. Dynamic procedures have been proposed to compute the Smagorinsky model coefficient locally within the flow domain. These dynamic procedures, in general, do not change the form of the model, but aim to adapt the model better to the local structure of the flow. The dynamic procedure formulated by Germano et al. [35] has been tested by Padilla and Silveira-Neto [36] and have resulted in reasonably well prediction of the transition to turbulence in natural convection in a horizontal annular cavity. Nevertheless, the dynamic model is not without drawbacks. The homogeneous averaging and ad hoc clipping procedures to prevent unstable numerical calculations have been found to be impractical in complex flows such as natural convection. At the same time, the Lagrangian dynamic model [37] and the dynamic localisation model [38] which have then been developed to overcome these drawbacks have resulted in huge computational overhead. This disadvantage inspired Piomelli and Liu [39] to propose an approximate localised dynamic procedure based on a time extrapolation technique. However, numerical experiments carried out by the authors have demonstrated the necessity of constant clipping to yield a well-behaved numerical calculation. Several other different SGS models have been developed to account for the shortcomings of the dynamic models, for instance the MILES approach which uses the dissipation terms introduced in the framework of flux-limiting finite volume discretizations [40]. Lately, an eddy-viscosity model developed by Vreman [41] has been shown to be able to predict small SGS dissipation for various laminar and transitional shear flows. It has been shown that the model with a fixed coefficient is able to predict accurate turbulent behaviours for wall-bounded channel shear flow without introducing any wall-damping functions. Even though subsequent tests have shown that the coefficient is not universal for turbulent channel flow [42–44], but the successful application of the Vreman model in premixed turbulent combustion [45] seems to suggest its useful application in natural convection.

Natural convection boundary layer in an enclosed cavity has attracted particular attention by various researchers due to the complex physics involving the simultaneous presence of laminar, transitional and turbulent flow regimes in a thermally-stratified flow. Interestingly, the dynamic model has been found by Peng and Davidson [46] to incorrectly capture the flow structure and hence producing huge discrepancies for the mean velocity field as compared with experimental data. In a later study, Peng and Davidson [47] who included a SGS buoyant production term based on the formulation of Eidson [48] in the dynamic model have found reasonable agreement with experimental measurement except in the outer part of the natural convective boundary layer. Meanwhile, in the natural convection studies of indoor airflow by Zhang and Chen [49,50], the dynamic model has been found to behave inadequately in predicting the turbulent kinetic energy near a wall. More recently, natural convection boundary layer in a confined cavity has been studied by Barhaghi and Davidson [51] through three SGS models; the Smagorinsky model, the WALE model, and the dynamic model. Through a clipping procedure for the Smagorinsky model coefficient, the dynamic model has been found to be better as compared to the Smagorinsky and WALE models in predicting the transition location. In this article, we focus on the study of natural convection in an open-ended vertical parallel-plate channel which represents even more challenging configuration for the performance of the SGS models. This configuration presents more complication as compared to the aforementioned natural convection in an enclosed cavity due to the presence of open boundaries where the natural entrainment of ambient fluid has to be appropriately treated. In the limited investigations reported for this problem [52–57], periodic boundary conditions either in the span-wise or stream-wise direction were imposed, resulting in the assumption of homogenous turbulence. Nonetheless, the results of Smith [58] who studied turbulent free convection adjacent to a vertical flat plate has revealed that a periodic flow structure, which is present in the early stages of transition, disappears as the intensities of temperature and velocity increase to a maximum in the mid-stage of transition. In addition, the strong turbulent energy production triggered by the interaction of the buoyancy and viscous forces due to the presence of solid walls obviates the homogeneity and isotropic turbulence assumption. This heterogeneous turbulence inherent in such flows results in a rather tough examination for the dynamic model in particular because, unlike forced turbulent channel flow where plane averaging can be performed in homogeneous directions, there is no identified homogenous direction in which averaging of the Smagorinsky model coefficient can be found. Thus, to our knowledge, no LES study has been reported for natural convection in an open-ended vertical parallel-plate channel without the assumption of periodic boundary condition. Therefore, in the present study, in order to investigate the effects associated with the presence of open boundaries and inhomogeneous flow behaviour, four SGS models have been examined a posteriori for natural convection in an asymmetrically-heated vertical parallelplate channel with a high aspect ratio. The SGS models studied are: (i) Smagorinsky, (ii) dynamic, (iii) approximate localised dynamic and (iv) Vreman. 2. Governing equations In LES, the Navier–Stokes equations governing the mass, momentum and energy conservation are spatially filtered. This filtering is represented mathematically in physical space as a convolution product:

 tÞ ¼ /ðx;

Z

1

1

/ðn; tÞGðx  nÞdn

ð1Þ

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

 tÞ is the resolved part of an arbitrary space–time variable where /ðx; while the convolution kernel G is the characteristic of the filter used. For the finite volume formulation considered in the present study, the box or top-hat filter is commonly used. Also, in threedimensional computations with grid cells of different sizes in the different Cartesian coordinate directions, the filter width of the subgrid length D has often been taken to be:

D ¼ ðDxDyDzÞ1=3 :

ð2Þ

In the present investigation, in order to account for the variation of density within the flow, the Favre-averaging is utilised, i.e.

~ tÞ ¼ q/ðx; tÞ: q /ðx;

ð3Þ

103

where PrT is given a value of 0.4. This value has been found to adequately model the turbulent thermal diffusion for buoyancy-driven flows where in the study of Barhaghi and Davidson [51], a similar value has been adopted for confined natural convection flows with sufficiently high Rayleigh numbers. Furthermore, in the numerical study of buoyancy-driven plumes by Pham et al. [25], similar range of values have also been revealed through direct numerical simulations. In the present work, the SGS Prandtl number has not been varied dynamically because the main objective of this work is to purely assess the SGS viscosity field in order to present a fair comparison between the SGS models. In addition, the focus is also on the prediction of SGS viscosity and its effect on various turbulent parameters.

The instantaneous space–time variable may now be written as 2.1. Smagorinsky model

~ tÞ þ /0 ðx; tÞ /ðx; tÞ ¼ /ðx;

ð4Þ

~ tÞ represents the Favre-averaged filtered or resolvable where /ðx; component and /0 (x, t) is the Favre-averaged subgrid-scale component which represents the unresolved spatial variables at a length smaller than the filter width D. The low-Mach-number Favre-filtered compressible mass, momentum and energy conservation equations in a Cartesian coordinate frame obtained when the filtering operation is applied are:

 @ðq u ~j Þ @q ¼0 þ @xj @t

ð5Þ

 @r ~ ij @M ij ~j Þ @ðq ~i u ~j Þ u u ! @ðq @p   qref Þ g þ ¼ þ þ þ ðq @t @xi @xj @xi @xi

ð6Þ

@ Te @ Te @ q C p þ q C p u~ i ¼ @t @xi @xi

lC p @ Te Pr @xi

! þ

u;T @s @xi

  ~i @ u ~j ~k @u 2 @u  lðTÞ þ : 3 @xj @xi @xk

ð7Þ

ð8Þ

The dynamic viscosity is calculated according to the Sutherland formula:

l ¼ 19:65  106

Te 298:15 þ C e 298:15 T þC

!3=2 ð9Þ

where C is the Sutherland constant. In the present investigation, C and Pr are prescribed with values of 120 and 0.7 respectively for air. u;T in Eqs. (6) and (7) represent the unknown The terms Mij and s SGS correlations obtained when the Favre filtering operation is applied to the governing equations. In the present investigation, they are modelled using the commonly used physical space models. The unresolved turbulent SGS momentum stress tensor Mij is modelled according to:



1 1  ui uj  s kk dij ¼ 2lT e M ij ¼ s S kk dij S ij  e 3 3

 ð10Þ

where lT is the modelled turbulent eddy-viscosity and e S ij ¼   ~j @u ~i 1 @u þ @xi . 2 @xj Meanwhile, the thermal flux vector due to SGS motion is correlated accordingly to the turbulent Prandtl number PrT:

~ i Te Þ  su;T ¼ qðg ui T  u

lT @ Te PrT @xi

lT ¼ q ðC S DÞ2 ð2eS ij eS ij Þ1=2

ð11Þ

ð12Þ

where Cs is the model coefficient. In practice, CS is typically in the range of 0.1–0.2 for most flows and is adjustable to improve results. For instance, Clark et al. [59] used CS = 0.2 for a case of isotropic homogeneous turbulence, while Deardorff [60] used CS = 0.1 for a plane channel flow. In the present study, due to the much debated dissipative nature of the Smagorinsky model [32], CS was prescribed a lower value of 0.125. Furthermore, in order to account for the flow behaviour near the wall, a wall function was used, resulting in the following formulation for the SGS eddy viscosity:

lT ¼ q L2S ð2eS ij eS ij Þ1=2

where q and T are the density and temperature respectively, while qref is the reference density at the reference temperature Tref, u~i is  is the pressure, Pr is the molecular the velocity in the ith direction, p ~ ij is the Prandtl number and l is the dynamic viscosity. In Eq. (6), r stress tensor based on the Stokes’ hypothesis:

r~ ij ¼ lðTÞ

In the Smagorinsky model, lT is expressed as:

ð13Þ

where LS = min (jd, flCSD) with fl representing a damping function and d denoting the distance from the wall. In the present investigation, the Van Driest damping where damping function can be written as:

fl ¼ 1  expðxþ =25Þ has been used while j was prescribed a value of 0.42. 2.2. Dynamic model Germano et al. [35] proposed a dynamic formulation for the Smagorinsky model coefficient by automatically adapting the model for laminar and transitional flows. This procedure results in a time- and space-dependent dynamic Smagorinsky constant Cd. This procedure eliminates the use of an ad hoc and a priori prescription of the Smagorinsky constant in the use of a SGS model for transitional and turbulent flows. The Germano model is based on an algebraic identity between the subgrid Reynolds stress tensor and the resolved turbulent tensors for two levels of the filtering process which is written in the form:

^~ij Lij ¼ T ij  s

ð14Þ

in which

^ ud ^ c g uei c uej T ij ¼ q i uj  q

ð15Þ

and

1 dd d ^ ue ej  q  uei q  uej Lij ¼ q iu q^

ð16Þ

~ij and The symbol (^) designates the second level filter. The tensors s Tij are the subgrid tensors corresponding, respectively, to the first and second filtering levels. The latter filtering level, termed the test b and can be filter, is associated with the characteristic length D obtained for the box filter as:

104

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

Z

1 ^~ ~ / /ðxÞdx i ¼ b bD D

ð17Þ

b is normally taken to be D b ¼ 2D. Assuming that the two where D subgrid tensors can be modelled by the same constant Cd for both filtering levels, the Eq. (14) can now be rewritten based on a least-squares method by Lilly [61]:

Cd ¼

M ij Ldij

ð18Þ

M kl M kl

in which

ð19Þ

and

M ij ¼ aij  c bij

ð20Þ 

1^ ^ e ^ b 2q ^ je S kk dij aij ¼ 2 D Sj S ij  e 3



  1  je and bij ¼ 2D2 q S kk dij Sj e S ij  e 3

Negative values can be obtained for the constant Cd, so the model is allowed to have an anti-dissipative effect locally. This characteristic is often interpreted as the modelling of the backward energy cascade mechanism, also known as backscatter phenemenon. Furthermore, it is not bounded, since it appears in the form of a fraction whose denominator can reach small values, resulting in large Cd. These two properties tend to give rise to unstable numerical computations. Hence, in order to ensure good numerical properties, many investigators have suggested an averaging procedure in homogeneous directions. Nevertheless, in the present investigation where no homogeneous directions could be identified, a local smoothing procedure is applied to Cd by averaging over neighbouring test-filtered grid cells. A clipping procedure is also applied, setting Cd to zero when the SGS viscosity is negative, thus preventing the anti-dissipative effect, i.e.,

! hMij Ldij i C d ¼ max 0; hMkl Mkl i

 Ldij ¼ C d aij  Cd bij

ð22Þ

where on the right-hand side, C⁄ is an estimate of the dynamic constant Cd, which is assumed to be known. The dynamic constant Cd is now evaluated as:  M ij ðLdij þ Cd bij Þ M kl M kl

ð23Þ

A backward extrapolation scheme has been suggested for the computation of C⁄, viz.,

C ¼ C n1 þ Dt

 @C  þ  @t n1

The huge computational overhead required by most dynamic models has motivated Vreman to formulate a different eddy-viscosity model in which its dissipation is relatively small in transitional and near-wall regions. This model does not involve explicit filtering, averaging, or clipping procedures. Vreman [41] proposed the following eddy viscosity:

lT ¼ q C v

sffiffiffiffiffiffiffiffiffiffi Bb

aij aij

ð27Þ

~ @u

with aij ¼ @xij ; bij ¼ D2m ami amj ; and Bb ¼ b11 b22  b212 þ b11 b33  b213 þ b22 b33  b223 . The model coefficient Cv is related to the Smagorinsky constant CS by C v  2:5C 2S and Dm is the filter width of subgrid length in the m direction. Similar to the Smagorinsky model, this model is easy to compute in actual LES, since it does not need more than the local filter width and the first-order derivatives of the velocity field. In the present investigation, Cv is prescribed a value of 0.1. As reported by Vreman [41], for tiny values of the eddy viscosity, machine precision errors accumulated from floating point operations may contaminate the calculation of Bb. For this reason, the simulations in the present study used lT = 0 if Bb < 108. 3. Numerical model

The numerical code considered herein is a Navier–Stokes solver for low-Mach number Newtonian, compressible and three-dimensional flows. The finite volume formulation was employed to discretize the filtered governing equations on a collocated grid. A non-diffusive fourth-order central differencing scheme was applied to approximate the convective fluxes so as to allow the natural propagation of disturbances and hence, a better representation of the SGS dissipation whilst a second-order central differencing scheme was adopted for the other spatial derivatives. The solution was advanced in time using an explicit two-step predictor–corrector approach, which handles the disparities of density effectively due to the strong coupling between density effects and fluid motion in the flow considered in the present study. This approach has been adopted from the original formulation of Knio et al. [62]. It involves a second-order Adams–Bashforth time integration scheme for the predictor stage and a second-order quasi Crank–Nicolson integration for the corrector stage. Pressure correction steps that were incorporated in both the predictor and corrector stages involve the inversion of pressure correction Poisson equations, which are solved using Krylov methods. Implementation of the numerical algorithm is summarised in Appendix A.

ð24Þ 3.2. Computational domain and boundary conditions

In the present investigation, a zeroth-order approximation was used for the extrapolation:

C ¼ C

In practice, the resulting dynamic procedure is fully local, and does not induce large extra computational effort as the original localised procedure described by Ghosal et al. [38].

3.1. Predictor–corrector methodology

In order to adapt the dynamic model better to the local structure of the flow, the approximate localised dynamic model has been formulated by Piomelli and Liu [39]. In this model, the expression in Eq. (19) can be recast in the following form:

n1

ð26Þ

ð21Þ

2.3. Approximate localised dynamic model

Cd ¼

* +!  Mij ðLdij þ Cd bij Þ C d ¼ max 0; Mkl Mkl

2.4. Vreman model

1 Ldij ¼ Lij  Lkk dij ¼ C d aij  Cd d bij 3

with

To ensure numerical stability, the coefficient Cd was locally averaged over the test cell:

ð25Þ

The numerical methodology described above has been applied to natural convective flow in an asymmetrically-heated vertical parallel-plate channel. This flow was investigated experimentally

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

105

(a)

(b)

Top wall

Top right opening

Top left opening

Span-wise wall

D Span-wise wall

Adiabatic wall

Heated wall H W Walls Openings

(0,0,0)

Span-wise wall y Bottom right opening x

Fig. 2. Comparison of axial velocity profiles between two different mesh resolution using the Vreman model at (a) y/H = 0.531 and (b) y/H = 0.774.

z Free-slip wall Bottom wall Span-wise wall

Fig. 1. Studied configuration: (a) the experimental apparatus used by Miyamoto et al. [63], (b) computational domain.

by Miyamoto et al. [63] and their experimental setup is shown in Fig. 1a. Data are available for the time-averaged flow field as well as various turbulent fluctuations. As shown in Fig. 1a, the experiment utilised two vertical channels consisting of three vertical parallel plates, with the plate in the centre heated. The apparatus was placed in an enclosure in such a way that flow was not induced in the span-wise direction. Different heat fluxes and channel widths were investigated; the case of qw = 104 W/m2 and W = 0.1 m were adopted in the present investigation.

The computational domain and the flow configuration considered here are represented in Fig. 1b. With reference to Fig. 1a, the domain dimensions are as follow: D = 0.96 m, W = 0.1 m and H = 4.98 m. They are exactly the same as those investigated by Miyamoto et al. [63] except for the two additional imaginary computational domains introduced at the top and bottom of the vertical channel to allow the natural entrainment of ambient air. In the present investigation, the x-, y-, and z-directions are the wallnormal, stream-wise and span-wise directions respectively. Owing to the great computational resources that may be required to simulate both channels, computational calculations have been restricted to only one channel. Thus, a free-slip wall was imposed at the bottom left end (x/W = 0, 0.092 < y/H < 0, 0 < z/D < 1) to mimic the effect of the other channel. The heated wall at x/W = 0 was imposed a temperature gradient of 104 W/m2. At the top and bottom boundaries where flow was induced or expelled, an opening boundary condition as described by Gresho [64,65] was employed.

106

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

This non-Dirichlet and traction-free boundary condition can be written as following:

Table 1

SGS model Smagorinsky Dynamic Approximate localised dynamic Vreman

ðrij þ pdij Þ  nij ¼ 0 where rij is the stress tensor as depicted by Eq. (8) while nij is the unit vector normal to the boundary. The advantage of this boundary condition is that the natural entrainment of ambient fluid is allowed. All the other boundaries were specified to be no-slip and adiabatic.

x+ 1.4 0.88 0.90 1.2

0:35 : maxðju=Dxj þ jv =Dyj þ jw=DzjÞ

3.3. Grid and time step resolution

Dt ¼

In order to accurately investigate and assess the effect of different SGS models in representing the losses of kinetic energy due to the filtering operation, it is important to study the effect of grid size in LES. In this study, two mesh distributions of 31  555  50 and 41  695  60 cells in the x-, y-, and z-directions were tested within the computational domain. A final mesh distribution of 41  695  60 cells was chosen because it was found that the difference between these two mesh distributions was not significant. Fig. 2a shows the axial velocity profile at y/H = 0.531 obtained using the Vreman model for both mesh distributions. It can be seen that when the number of grid points were almost doubled, the numerical results between the coarse and fine grid resolutions for the time-averaged axial velocity profile was found to be less than 3%. The discrepancies can be explained by the transition phenomenon which presents great complexity in modelling for turbulent natural convection. This mesh dependency of the numerical results within the transition region could be a result of, among many factors, the SGS model used, and has been highlighted by several other investigators [51,53]. In addition, the transition which is measured experimentally is also subjected to great difficulties. Meanwhile, a further investigation of the axial velocity profile at a different height (y/H = 0.774) is presented in Fig. 2b for the coarse and fine grid resolutions. A difference of no more than 0.1% has been achieved. Thus, it could be verified that the numerical solution obtained with these mesh distributions is sufficiently independent of grid spacing and hence filter width. In addition, the sufficiency of the grid spacing was validated against existing scaling relationships for natural convection adjacent to a vertical plate. In natural convection, instabilities which develop in the laminar boundary layers quickly propagate and trigger the onset of turbulence if certain critical conditions are met [66]. Hence, these instabilities in the boundary layer should at least be captured by the smallest grid size. For a non-stratified flow, the boundary layer thickness scale in the quasi-steady phase of the flow is given by the following relationship [67]:

In the present article, the minimum time step captured is in the or-

d

y Ray1=5

:

ð28Þ

For the Rayleigh number considered in the present investigation, the smallest boundary layer thickness scale is approximately of the order of 3  10–3 m. Thus, the mesh distribution of 41  695  60 which resulted in the minimum grid size of 2.5  10– 3 m was considered sufficiently fine for practical applications. Presented in Table 1 is the comparison between different SGS models for the shear velocity in terms of wall units. Hence, the mesh distribution used in the current numerical study is considered practically sufficient with adequate nodes placed within the boundary layer. A transient analysis was performed with approximately a total of 60,000 time steps for all SGS models in order to fully capture the turbulent statistics of the quasi-steady flow. Stable time stepping was accomplished through the acoustic CFL condition. The time step was determined by employing a CFL number of 0.35 according to:

der of 8  10–4 s. When compared with the scaling relation t 

ð29Þ

y2

aRa2=5 y

[67], it has been demonstrated that the present computer code fully captured the time scale growth of the thermal boundary layer. 4. Results and discussion In the present investigation, four sampling lines were placed at different axial locations within the channel at the centreline (z/ D = 0.5) i.e. at y/H = 0.267, 0.531, 0.774, and 0.982, with three of them corresponded to the locations sampled by Miyamoto et al. [63]. Simulation was started with a quiescent fluid at an ambient temperature of 25 °C. Great care was then taken to observe the development of the flow. Once the system reached a quasi-steady state, statistical averaging was performed till sufficient amount of samples were collected. Fig. 3a and b shows the axial velocity and temperature profiles respectively at y/H = 0.778 and z/D = 0.5 obtained using the Vreman model. As can be seen, when the profiles are time-averaged separately for 30 s < t < 40 s, 40 s < t < 50 s, and 50 s < t < 60 s, very little differences are found between the results for both the velocity and temperature values. Again, it can be seen that the simulation has reached the fully developed state. Hereafter, the statistics and time-averaged results are obtained from t = 30 s till t = 60 s. 4.1. Time-averaged wall temperature distribution The existence of a solid wall represents a challenging aspect for the performance of SGS model since a good model ought to be able to predict the correct near-wall behaviour in wall-bounded shear flows. Furthermore, the coupled effect of the buoyancy production and viscous dissipation adjacent to the heated wall which drives the onset of turbulence has to be predicted well in order to capture the correct transition location. The time-averaged wall temperature excess distribution (T  Ti) where Ti is the inlet air temperature, obtained in the present simulation at the centre-line of both the heated and adiabatic wall (z/ D = 0.5) is shown in Fig. 4 for all SGS models considered together with experimental results previously cited. As observed experimentally, the heated wall temperature excess distribution shows a gentle, yet observable peak at about y/H = 0.4, after which the wall temperature decreases. Further downstream, the wall temperature shows an increasing trend, albeit at a lower rate. This peak in temperature distribution along the heated plate has been characterised by Miyamoto et al. [63] to be the transition point towards turbulence based on a previous study on natural convection adjacent to a vertical flat plate [68]. On the other hand, for the insulated wall, a slight increase in the wall temperature is observed both due to the strong radiation losses as reported experimentally and the rapid growth of the turbulent thermal boundary layer. The wall temperature excess distribution presented in Fig. 4 allows us to analyse in more detail the behaviour of the heated wall

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

107

Fig. 4. Comparison between numerical and experimental wall temperature excess distribution along heated and adiabatic walls.

Fig. 3. Comparison between contiguous time-averaged data at y/H = 0.778: (a) axial velocity profile and (b) axial temperature profile.

temperature excess distribution predicted by the SGS models considered in the present work. As can be seen in Fig. 4, the Smagorinsky model under predicts the temperature values notably. Furthermore, the observed negative temperature gradient in the experiment has not been well captured by the Smagorinsky model. This result implies that the Smagorinsky model which has been found to be too dissipative for wall-bounded shear flows does not predict the correct near-wall asymptotic behaviour, resulting in the lower temperatures predicted. At the same time, when a comparison is made between the temperature values obtained numerically using the dynamic model against experimental data, it is found that the temperature values at locations upstream of the channel for y/H < 0.2 are well predicted. However, the negative temperature gradient at about y/H = 0.4 in the temperature excess distribution as observed experimentally has not been captured well by the dynamic model. On the contrary, over prediction of wall temperature is observed at all locations above y/H = 0.2. This finding suggests that the SGS dis-

sipation computed using the dynamic model based on the local volume averaging of neighbouring cells has been predicted to be lower, resulting in under estimation of the total physical dissipation, thus producing unexpected higher temperature values at the heated wall. When compared to the standard dynamic model, the approximate localised dynamic model not only can capture the correct range of temperature values at the upstream regions (y/H < 0.2), but also the negative temperature gradient as measured by Miyamoto et al. [63]. However, the negative temperature gradient in the temperature excess distribution is found to be predicted at a higher location (y/H  0.6) while temperature values downstream of the channel are still over predicted in a similar manner to the results predicted by the dynamic model. As pointed out by Sagaut [32], these dynamic procedures of computing the Smagorinsky constant do not generally change the formulation of the model, but only allow a better adaptation to the local flow dynamics. This reason explains the discrepancy found in the wall temperature prediction which is a result of the inaccuracy inherent from the original formulation of the Smagorinsky model in predicting near-wall behaviour. On the other hand, the Vreman model has resulted in good agreement with experimental data without introducing any wall model or dynamic procedure. This is mainly due to the different formulation used in modelling the SGS eddy viscosity which is able to predict vanishing or small eddy viscosity in near-wall regions. This favourable property of the model has a significant effect particularly when predicting transitional flows which exhibit weak turbulence, resulting in the well-captured transition point as characterised by the negative temperature gradient. However, a discrepancy between experimental and numerical temperature values of more than 1 °C is observed in the vicinity of the inlet of the channel and in the transitional region which. In addition, when comparing the numerical values obtained from the Vreman model with those obtained using both dynamic models, a slight over prediction can be seen at the inlet. This observation suggests that the value of the model coefficient may have been over prescribed in the present study which could have over dampened the disturbances in upstream laminar regions, resulting in slightly higher temperature values predicted.

108

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

As to the adiabatic wall, the temperature excess distribution shows a discrepancy of order 1–3 °C between experimental measurements and numerical values obtained using all SGS models. As can be seen in Fig. 4, the largest discrepancy was found when the Smagorinsky model was used while both dynamic models have been able to predict higher temperature values downstream. At upstream laminar regions, all SGS models have failed to account for the temperature increase as measured experimentally. This under prediction of the wall temperatures could be attributed to the radiation heat exchange which is not considered in the present work.

4.2. Time-averaged axial velocity and temperature profiles To assess the SGS models in terms of mean velocity field, numerically obtained axial dimensionless velocity profiles at y/ H = 0.165, 0.531, and 0.774 were compared with experimental data and plotted in Fig. 5. As can be seen from Fig. 5a, experimental measurements at y/H = 0.165 reveal that the velocity adjacent to the heated wall is slightly higher due to the asymmetric heating configuration. The profile then flattens a little before it starts to increase slightly towards the insulated wall. This increase in velocity at the insulated wall can be accounted for by radiation heat exchange which appears to be significant especially for laminar regions. Further downstream at y/H = 0.531 (Fig. 5b) and y/ H = 0.774 (Fig. 5c), it can be noted that the velocity peak adjacent to the heated wall exhibits higher values. In terms of the evaluation of the performance of different SGS models, as presented in Fig. 5a, the Smagorinsky model results in large discrepancies when compared to both experimental measurements and numerical results obtained by the other SGS models. Firstly, the velocity peak adjacent to the heated wall has been over predicted considerably. Secondly, the velocity values for the bulk flow and the shape of the profile do not match well at all with experimental data. This observation can be explained by the large eddy viscosity values obtained by the Smagorinsky model which produce unphysical dissipation towards the flow at laminar regions where dissipation is expected to be small. On the other hand, the dimensionless axial velocity profile at y/H = 0.165 has been reasonably well predicted by the Vreman, dynamic and approximate localised dynamic models when compared against experimental data. The slightly lower velocity values predicted adjacent to the adiabatic wall could be explained by radiation heat exchange which is not considered in the present investigation. Again, in the laminar region where the fluid exhibits a well-behaved pattern, the strong radiation loss from the heated plate has caused the temperature of the fluid adjacent to the insulated plate to be increased slightly, resulting in the higher velocity as measured experimentally. Further downstream at y/H = 0.531 (Fig. 5b), the Smagorinsky model results in an over prediction of the velocity peak near the heated plate while adjacent to the adiabatic plate, velocity values are slightly under predicted. This can be explained from the fact that in transitional flows where turbulence is rather weak particularly in regions outside of the boundary layer, the dissipative formulation inherent in the Smagorinsky model which results in isotropic or uni-directional energy transfer from the large scales to the subgrid scales tends to laminarise the flow in the boundary layer region. Furthermore, in the near-wall boundary layer where the onset of transitional and turbulent flow is triggered, the dissipative character of the Smagorinsky model might have dampened the turbulence which is expected to occur at that location. Consequently, the model has a tendency to predict a more laminar profile with a relatively narrow boundary layer and a high velocity

Fig. 5. Comparison between numerical and experimental dimensionless axial velocity profile at three different heights: (a) y/H = 0.165, (b) y/H = 0.531 and (c) y/ H = 0.774.

109

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

(z/D = 0.5) for all SGS models considered in the present work. It is noted that the true scale of this figure has not been retained due to the high aspect ratio considered in the present study. As can be seen, the thermal boundary layer which starts to develop adjacent to the heated plate increases in thickness along the height of the channel. Assuming that transition to turbulence is expected to occur at about y/H = 0.4 as suggested by Miyamoto et al. [63], then the thickness of the thermal boundary layer would be expected to grow rapidly thereafter. Consequently, where the aspect ratio of the channel is high, the thermal boundary layer is expected to hit the opposite adiabatic wall. Regrettably, as can be seen from Fig. 6a, the thermal boundary layer predicted by the Smagorinsky model shows slower growth as compared to the other SGS models, thus emphasising the delayed transitional behaviour of the model. The Vreman model (Fig. 6d) has predicted a rapid growth of the thermal boundary layer from approximately y/H = 0.5 onwards while interestingly, both dynamic models have resulted in a much earlier growth of the boundary layer, suggesting the unnecessary chaotic structures produced at upstream regions as will be further detailed in Sections 4.3 and 4.4. 4.3. Turbulent parameters Fig. 7a and b presents the comparison of numerically computed stream-wise and wall-normal RMS velocity for qw = 104 W/m2 with experimentally measured data at y/H 0.774 for qw = 101 W/ m2. As can be seen from Fig. 7a, the trend in the profile of vrms has not been well predicted by the Smagorinsky model particularly

y/H

y/H

y/H

y/H

peak which is exactly as shown in Fig. 5b. Meanwhile, both dynamic models have slightly under predicted the velocity peak as compared to experimental data. The reason for this can be that the small SGS dissipation predicted by both dynamic models has unnecessarily generated high unsteadiness within the flow, causing the flow to fluctuate in an anisotropic manner as will be detailed in Section 4.3. On the contrary, the Vreman model results in good agreement with experimental data due to the accurate prediction of SGS dissipation which closely resembles the expected physical dissipation. It is only at y/H = 0.774 (Fig. 5c) which corresponds to the fully turbulent regime as suggested by Miyamoto et al. [63] that the Smagorinsky model has been found to reasonably predict the profile with a slight under prediction towards the adiabatic wall. This finding suggests that not only the model behaves inadequately in near-wall regions, but it also predicts delayed transitional flow behaviour. At the same time, both dynamic models have produced slightly lower velocity values, particularly at the heated wall when compared against experimental data. This under prediction of velocity near the heated wall indirectly explains the over prediction of wall temperature as observed previously in Fig. 4 where the heat generated has not been well convected by the bulk flow. Once again, it can be observed that the Vreman model results in the best agreement with experimental data among all SGS models tested. The delayed transitional behaviour of the flow predicted by the Smagorinsky model is further emphasised in Fig. 6. Fig. 6a–d presents the time-averaged temperature contour at the centre plane

x/W

x/W

x/W

x/W

(a)

(b)

(c)

(d)

Fig. 6. Comparison between time-averaged temperature contours obtained using four different SGS models: (a) Smagorinsky model, (b) dynamic model, (c) approximate localised dynamic model, and (d) Vreman model.

110

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

vrms

urms

(a)

Fig. 8. Comparison between resolved-scale turbulent kinetic energy k at y/ H = 0.778 between different SGS models.

wrms

v 'u '

(b) Fig. 7. Comparison between numerical and experimental RMS velocities: (a) and urms, and (b) wrms.

vrms Fig. 9. Comparison between the profile of turbulent production by shear H = 0.778 between different SGS models.

in regions adjacent to the heated wall. The profile predicted by this model shows a peak at the heated wall which has not been observed both experimentally and numerically in all the other SGS models considered. As for the profile predicted by both dynamic models, the peak of vrms has been found to be rather large near the adiabatic plate. The approximate localised dynamic model has been found to be able to capture the RMS velocity values better within the bulk flow as compared to the standard dynamic model. This has demonstrated that the model is able to adapt better to the local flow dynamics as compared to the standard dynamic model. At the same time, the general trend in the profile of vrms has been reasonably well predicted by the Vreman model where a flat profile within the bulk flow is observed. Slightly higher values found towards the adiabatic plate could be a result of the higher Rayleigh number considered in the current numerical study. As to the prediction of urms, experimentally with a lower heat flux reported, the profile exhibits a rather flattened peak slightly towards the heated wall. The described ‘loose peak’ by Miyamoto

v 0 u0

at y/

et al. [63] in the centre has been predicted to be slightly higher for all SGS models due to the higher heat flux imposed in the present study. As can be seen, almost all SGS models yield reasonably close values with each other. However, the slightly higher fluctuations given by both dynamic models as compared to the static models appear to demonstrate their less dissipative behaviour. The ad hoc clipping procedure used in these models might introduce unnecessary perturbations in the flow field, resulting in the higher random fluctuations seen in Fig. 7 as compared to the static models. In addition, during the numerical calculation, it was also observed that the fluctuations of the Smagorinsky constant in the dynamic procedures were large. This ad hoc property of the dynamic models suggests limited robustness as compared to the Vreman model which only requires a fixed constant. Even though not reported experimentally, we present here the profile of wrms at y/H = 0.774 across the channel width (Fig. 7b).

111

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

Trms

y/H

‘Eruption’

(a) v 'T '

‘Inrush’

x/W Fig. 11. Instantaneous velocity vector and temperature contour at z/D = 0.5 for the section 0.65 < y/H < 0.85 using the Vreman model.

(b) u 'T '

(c) Fig. 10. Comparison between different SGS models at y/H = 0.778: (a) Trms, (b) and (c) v 0 T 0 .

v0T0

As can be seen, the Vreman, dynamic and approximate localised dynamic models have predicted reasonably close results, while the Smagorinsky model yields lower values near the heated wall. This lower values predicted can be verified by virtue of the conservation of mass if we consider the higher peak in the profile of vrms predicted by the Smagorinsky model comparing with the other SGS models. Evident in the profiles of the RMS velocities shown previously in Fig. 7a and b, all three components of the RMS velocity show different values obtained in the three different directions, suggesting that the turbulence generated might be anisotropic. This explains the discrepancy obtained between experimental measurement and numerical results using the Smagorinsky model for the wall temperature excess distribution. Again, it is shown that the formulation of the Smagorinsky which is based on the isotropic framework produces incorrect flow behaviour especially when the flow exhibits strong anisotropy. Reasonably large magnitude of the span-wise velocity fluctuations wrms as shown in Fig. 7b suggests that three-dimensional effects are not negligible. In LES simulations, it is always of interest to estimate the turbulent kinetic energy carried by the resolved or large scales. Fig. 8 presents the profile of the resolved-scale turbulent kinetic energy profile across the channel width at y/H = 0.774 for all four SGS models considered in the present study. As can be seen from the figure, the resolved-scale turbulent kinetic energy predicted by the Smagorinsky model depicts one which has a relatively higher peak towards the heated wall. This trend is not observed for the

112

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

profile predicted by the other SGS models. On the other hand, both of the dynamic models, together with the Vreman model has resulted in almost uniform profile across the channel width. Also, it is observed in Fig. 8 that the Vreman model and the approximate localised dynamic model result in almost similar profile in the bulk of the flow, except adjacent to the heated wall where the Vreman model predicts slightly smaller values. Meanwhile, the standard dynamic model is seen to yield slightly larger resolved-scale turbulent kinetic energy as compared to the approximate localised dynamic model and the Vreman model. This finding further suggests that the ad hoc clipping procedure used in the dynamic model might have resulted in unnecessary fluctuations within the flow field due to the large variation of the Smagorinsky constant. Fig. 9 presents the profile of v 0 u0 which represents the turbulent production by shear at y/H = 0.774 across the channel width for all four SGS models considered. It can be seen that all four models predicted a slight negative value in regions adjacent to the heated wall, consistent with the findings of various researchers in the case of natural convective boundary layer flow [47,51,53–57]. These slight negative values are due to the buoyancy generated which produces an upward rise for the warmer plumes of air adjacent to the heated wall. By the conservation of mass, air in the outer part of the boundary layer then exhibits a tendency to be induced towards the heated wall, thus resulting in the slight near-wall negative values. As can be seen, when comparing the profile obtained

using different SGS models, the Smagorisnky model, unlike the rest of the SGS models, under predicts the Reynolds stress considerably for the bulk flow. This can be explained by the dissipative formulation of the model which tends to dampen instabilities instead of propagating them towards the outer part of the boundary layer. In addition, as for both dynamic models, it can be seen that the approximate localised dynamic model has predicted lower turbulent shear production as compared to the standard dynamic model. Fig. 10a presents the temperature fluctuation Trms distribution across the channel width at y/H = 0.778 for all four SGS models considered. It can be seen that both the dynamic and the approximate localised dynamic models predict high temperature fluctuations as compared to the static Vreman and Smagorinsky models. In addition, it can also be seen that the local maxima predicted by the dynamic models are relatively closer to the wall as compared to the Vreman model while the Smagorinsky model predicts the local maxima further away from the solid boundary. Further away from the boundary layer, the profiles of Trms predicted by all four models are approximately in the same range with a slight over estimation by the dynamic models. As for the prediction of the turbulent heat fluxes, Fig. 10b and c show the turbulent heat flux distribution u0 T 0 and u0 T 0 across the channel width at y/H = 0.778 in the y- and x-direction respectively. As can be seen, similar trends as compared to the distribution of Trms have been predicted by the SGS models. Both versions of the dynamic model predict higher turbulent heat fluxes as compared

y/H

y/H

y/H

y/H

µSGS µref

x/W

x/W

x/W

x/W

(a)

(b)

(c)

(d)

Fig. 12. Instantaneous contour plot llSGS of using different SGS models: (a) Smagorinsky model, (b) dynamic model, (c) approximate localised dynamic model, and (d) Vreman ref model.

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

This complex interacting mechanism thus forms the ‘hairpin’ or ‘horseshoe’ vortex in most turbulent boundary layer flows, resulting in the generation of large Reynolds stresses. 4.4. Subgrid-scale dissipation

y/H

The successful implementation of LES essentially lies in the accurate computation of the SGS dissipation to represent the total physical turbulent dissipation. To investigate this aspect, the ratio of the instantaneous SGS turbulent viscosity to the reference dynamic viscosity llSGS was computed and presented in Fig. 11 at the ref centre plane (z/D = 0.5) of the channel for the Smagorinsky model (Fig. 12a), dynamic model (Fig. 12b), approximate localised dynamic model (Fig. 12c), and the Vreman model (Fig. 12d). As can be observed in Fig. 12a, the Smagorinsky model has yielded large ratio of llSGS in near-wall regions as well as the inlet of the channel. ref This finding further explains the erroneous prediction of wall temperature values as shown previously in Fig. 4. In addition, since the transition towards turbulence in natural convective flows is often triggered by the instabilities which propagate from the inlet and near-wall boundary layer, the large eddy viscosity values obtained have dampened the necessary perturbations. This clearly explains the delayed transitional behaviour as seen previously in the axial dimensionless velocity profiles (Fig. 5b). This theoretical aspect implies that a refined grid is desired when using a dissipative model such as the Smagorinsky model for the prediction of the onset of

y/H

y/H

to the Vreman and Smagorinsky models. These predictions by both dynamic models confirm that the large fluctuation of the model coefficient due to the ad hoc clipping procedure has resulted in higher random fluctuations in the velocity and temperature fields. These higher turbulent heat fluxes could explain the slightly lower velocity profile obtained when these dynamic models are used (Fig. 5c) due to the less convection of heat generated within the bulk flow at the heated wall. Furthermore, the higher turbulent heat fluxes obtained through these dynamic models could also explain the observed higher temperature values at the adiabatic wall because of the anisotropic mixing yielded as a result of the variation of the model coefficient. The observation of various turbulent fluctuations sheds light into the turbulence dynamics within the channel. Fig. 11 shows a snapshot for the instantaneous temperature contour and velocity vector at the centre plane (z/D = 0.5) within the channel for the section 0.65 < y/H < 0.85 obtained using the Vreman model. This figure shows an ‘eruption’ and an ‘inrush’ coherent structures which are important to the understanding of turbulent boundary layers. As can be seen, an eruption originates close to the wall where there are separated regions of fluid with different velocities due to the density differences. In the present study, this eruption is also found to be characterised by a relatively higher temperature. Subsequently, cooler air in the regions away from the wall has a tendency to be induced in the direction of the heated wall by virtue of the conservation laws, resulting in the observation of an ‘inrush’.

113

x/W

x/W

x/W

(a)

(b)

(c)

Fig. 13. Streaks originating from near-wall regions as illustrated by instantaneous vorticity iso-surface coloured with temperature obtained using the Vreman model: (a) 0 < y/H < 0.32, (b) 0.32 < y/H < 0.64, and (c) 0.64 < y/H < 0.96.

114

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

turbulence, in which the amplification of the large-scale turbulence needs to be well accounted for with appropriate energy scatter from the resolve large scales to the subgrid scales. Such a grid refinement, nonetheless, could substantially increase the computational requirements of LES similar to those of DNS, resulting in limited practical applications. At the same time, both the standard dynamic and the approximate localised dynamic models result in small ratios of llSGS near the ref wall. This favourable feature has allowed both dynamic models to capture the inlet temperatures reasonably well. However, as can be seen in Fig. 12b and c, the predicted SGS viscosity shows large variation within the computational domain which has resulted in unnecessary disturbances being generated particularly at upstream laminar regions. This result could be due to the ad hoc clipping procedure used in these models. On the other hand, the Vreman model has been capable of predicting low SGS eddy viscosity in upstream regions, thus accurately representing the onset towards turbulence in wall-bounded natural convective shear flow and its near-wall boundary layer behaviour. This observation further suggests possible introduction of implicit damping effect when the model is used. This aspect of the model may be due to the fundamental derivation of the Vreman model itself which results in the correct prediction of the near-wall flow dynamics. Unlike the Smagorinsky model, which requires either a wall function or a dynamic procedure to provide the implicit damping effect, the Vreman model can adapt itself towards the wall, suggesting promising application in various complex turbulent problems. Interestingly, as presented in the preceding paragraphs, rather large values of the SGS viscosity have been predicted at the inlet where flow is entrained into the channel by all SGS models, except for the Smagorinsky model. This observation suggests that the flow exhibits rather high unsteadiness even at the inlet due to a recirculation zone of a reasonable size. This is demonstrated by Fig. 13 which shows the streaks in near-wall regions illustrated by the instantaneous iso-surface of vorticity x = 12.5 s–1 coloured with temperature at three sections along the height of the channel obtained using the Vreman model. As can be seen from Fig. 13a, vortices initially are generated near both the heated and adiabatic walls as expected because instabilities and disturbances for boundary layer flows are mainly induced by the presence of solid walls. These instabilities then propagate downstream, causing the formation of chaotic flow structures responsible for the enhanced swirling and mixing of the fluid. This is clearly presented by Fig. 13c where the temperature of the vortices generated adjacent to the adiabatic wall is found to be higher in comparison to the results shown in Fig. 13a and b. In addition, it can also be observed that the vortices generated adjacent to the adiabatic wall tend to travel towards the heated wall, resulting in air with relatively lower temperature being drawn towards the heated wall, thus producing the negative temperature gradient as observed previously in Fig. 5 which has been considered to be the onset of turbulence.

5. Conclusion A thorough assessment of the performance of different LES SGS models for the simulation of buoyancy-driven flow in an open-ended vertical parallel-plate channel has been carried out a posteriori. Based on the comparison with experimental results, the Smagorinsky model has been found to yield large SGS dissipation, resulting in the prediction of inaccurate near-wall flow profiles and delayed transitional behaviour. At the same time, both of its dynamic formulations, however, have resulted in under prediction of SGS dissipation, thus producing excessive fluctuations in the flow field. This finding has further confirmed the

inadequacy inherent in the formulation of the Smagorinsky model. On the other hand, it was shown that the Vreman model has produced excellent agreement with experimentally measured data for both the temperature and velocity field profiles. The slight discrepancies in wall temperatures observed at the channel inlet seem to suggest that the Vreman model coefficient is not universal in all flow types and thus may necessitate the implementation of a dynamic procedure to compute the model coefficient. Appendix A The implementation of the numerical algorithm used in the present investigation is summarised here. Based on the ideal gas e , where  RT law, the equation of state can be expressed as pð0Þ ¼ q (0) the zeroth-order pressure p is assumed to be equivalent to the atmospheric pressure and R is the gas constant. For the purpose of numerical implementation, the time rate of change of density is found by differentiating the equation of state:

 @q q @ Te ¼ @t Te @t

ðA1Þ

A1. Predictor stage 1. Recalling the conservation equation for energy and Eq. (A1), the e =@tjn and @ q  =@tjn can be evaluated. local time derivatives of @ T 2. Predicted values for the density are determined accordingly by:

q   q n Dt

¼

 n 1 @q  n1 3 @q j  j 2 @t 2 @t

ðA2Þ

and the predicted temperature distribution is found from the equation of state:

  R Te  ¼ pð0Þ =q

ðA3Þ

3. An intermediate velocity field is then determined by integrating the pressure-split momentum equations:

q  u~0j  qn u~ nj Dt

q  u^~j  q  u~ 0j Dt

¼

n 3 n 1 n1 @ p Ri  Ri  2 2 @xj

ðA4Þ

¼

n @p @xj

ðA5Þ

where

    ~ @u ~ ~ ~i u ~j Þ u @ðq @ @u 2 @u þ l iþ j  l k @xi @xi 3 @xk @xj @xi ui uj @s !   qref Þ g :  þ ðq @xi

Ri ¼ 

ðA6Þ

Note that when Eqs. (A4) and (A5) are combined, the original formulation of Knio et al. [62] for the intermediate velocity without the pressure gradient is recovered which represents the standard fractional-step technique frequently employed. Establishing boundary conditions for these intermediate velocities is a common source of ambiguity. In time-splitting methods, only the boundary conditions for the velocity field are given at each complete time-step and those of the intermediate velocity field are unknown. Kim and Moin [69] demonstrated that proper boundary conditions are required to be consistent with the governing equations or else the solution suffers from appreciable numerical errors. Nonetheless, the introduction of an additional step of Eq. (A5) removes the need to construct boundary conditions for the intermediate velocity field since the pressure gradient is accounted in Eq. (A4); boundary

115

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

conditions that are employed for the velocity field can be immediately applied for the intermediate velocity field. 4. The predicted velocity field is then obtained using the projection step:

q  u~ j  q  u^~ j Dt

 @p ¼ @xj

ðA7Þ

  p n , Eqs. (A5) 5. By introducing the pressure correction p0 ¼ p and (A7) are combined to yield:

q  u~ j  q  u~ 0j Dt

¼

@p0 @xj

ðA8Þ

The predicted velocity field is now obtained through Eq. (A8) instead of Eq. (A7) in the projection step: The pressure correction p0 is the solution of the Poisson equation:

   ~0 uj Þ @ q   @ 2 p0 1 @ðq ¼ þ j Dt @t @xj @x2j

ðA9Þ

 =@tj is given by a first-order discretization: where @ q

~  @q 1 n  n1   q j ¼ ½q @t Dt

ðA10Þ

A2. Corrector stage 6. The temporal derivatives of the scalar fields at the new time level tn+1 are estimated based on the predicted values, and corrected values for the scalar fields are obtained through a mixed non-spilt scheme. The density in this corrector step can be obtained as:

q nþ1  q n Dt

1 n 1 ðC þ Dnq þÞ þ ðC q þ Dq þÞ 2 q 2

¼

ðA11Þ

where

Cq ¼

q ~ @ Te Te

ui

@xi

"  ui T 1 @s @ Dq ¼   @xi C p Te @xi

lC p @ Te Pr @xi

!# ;

The new temperature distribution is found using the equation of state:

 nþ1 Rnþ1 Te nþ1 ¼ pð0Þ =q

ðA12Þ

7. A second intermediate velocity field is then determined using the pressure-split momentum equations,

q  u~ 0j  qn u~ nj Dt

¼

n 3 n 1 n1 @ p Ri  Ri  2 2 @xj

ðA13Þ

8. The pressure distribution at the new time level is obtained through the solution of the Poisson equation:

   ~0 uj Þ @ q  nþ1 @ 2 p0 1 @ðq ¼ þ j Dt @t @xj @x2j

ðA14Þ

  nþ1 @q 1 nþ1 ¼ q  q n j @t Dt

ðA15Þ

and

9. Finally, the predicted velocity field at the new time level is determined using:

~ 0j u q nþ1 u~ nþ1 q j

Dt

¼

@p0 @xj

ðA16Þ

In the present study, the intermediate face velocities that are required to be evaluated prior to solving the Poisson equation for the pressure correction can be determined by the following:

~ 0j Þf ð u

q

n ~n uj Þf

¼ ðq

"  n #  3 n 1 n1 @p þ Dt Ri;f  Ri;f  2 2 @xj f

ðA17Þ

~ nj Þf ; Rni;f and Rn1 The terms represented by ðqn u can simply be obi;f tained by linear interpolation whilst the face pressure gradient n =@xj Þf is approximated by a second-order central difference beð@ p tween the cell centres. References [1] DesJardin PE, Frankel SH. Large eddy simulation of a nonpremixed reacting jet: application and assessment of subgrid-scale combustion models. Phys Fluids 1998;10:2298–314. [2] Satish M, Ma F, Islam MR. Large eddy simulation of a three dimensional buoyant jet. In: 2005 ASME international mechanical engineering congress and exposition. Orlando, Florida USA; 2005. p. 1321–30. [3] Addad Y, Benhamadouche S, Laurence D. The negatively buoyant wall-jet: LES results. Int J Heat Fluid Flow 2004;25:795–808. [4] Muller SB, Kleiser L. Large-eddy simulation of vortex breakdown in compressible swirling jet flow. Comput Fluids 2008;37:844–56. [5] Labbe O, Maglaras E, Garnier F. Large-eddy simulation of a turbulent jet and wake vortex interaction. Comput Fluids 2007;36:772–85. [6] Dinesh KKJR, Kirkpatrick MP. Study of jet precession, recirculation and vortex breakdown in turbulent swirling jets using LES. Comput Fluids 2009;38:1232–42. [7] Raufeisen A, Breuer M, Botsch T, Delgado A. LES validation of turbulent rotating buoyancy- and surface tension-driven flow against DNS. Comput Fluids 2009;38:1549–65. [8] Alkishriwi N, Meinke M, Schroder W. Large-eddy simulation of streamwiserotating turbulent channel flow. Comput Fluids 2008;37:786–92. [9] Zang Y, Street RL, Koseff JR. A dynamic mixed subgrid-scale model and its application to turbulent recirculating flows. Phys Fluids 1993;5:3186–96. [10] Kobayashi H, Ham F, Wu X. Application of a local SGS model based on coherent structures to complex geometries. Int J Heat Fluid Flow 2008;29:640–53. [11] Kim D, Yang K, Senda M. Large eddy simulation of turbulent flow past a square cylinder confined in a channel. Comput Fluids 2004;33:81–96. [12] Ouvrard H, Koobus B, Dervieux A, Salvetti MV. Classical and variational multiscale LES of the flow around a circular cylinder on unstructured grids. Comput Fluids 2010;39:1083–94. [13] Ghosh S, Sesterhenn J, Friedrich R. Large-eddy simulation of supersonic turbulent flow in axisymmetric nozzles and diffusers. Int J Heat Fluid Flow 2008;29:579–90. [14] Genin F, Menon S. Studies of shock/turbulent shear layer interaction using large-eddy simulation. Comput Fluids 2010;39:800–19. [15] Garnier E, Sagaut P, Deville M. Large eddy simulation of shock/homogeneous turbulence interaction. Comput Fluids 2002;31:245–68. [16] Wang T, Bai JS, Li P, Liu K. Large-eddy simulation of the Richtmyer–Meshkov instability of rectangular interfaces accelerated by shock waves. Sci China Phys Mech Astron 2010;53:905–14. [17] Bai J-s, Liu J-h, Wang T, Zou L-y, Li P, Tan D-w. Investigation of the Richtmyer– Meshkov instability with double perturbation interface in nonuniform flows. Phys Rev E 2010;81:056302. [18] Darmana D, Deen NG, Kuipers JAM, Harteveld WK, Mudde RF. Numerical study of homogeneous bubbly flow: influence of the inlet conditions to the hydrodynamic behavior. Int J Multiphase Flow 2009;35:1077–99. [19] He Y, Deen NG, Annaland MvS, Kuipers JAM. Gassolid turbulent flow in a circulating fluidized bed riser: numerical study of binary particle systems. Ind Eng Chem Res 2009;48:8098–108. [20] Choi J, Tawhai MH, Hoffman EA, Lin C. On intra- and intersubject variabilities of airflow in the human lungs. Phys Fluids 2009:21. [21] Lin C-L, Tawhai MH, McLennan G, Hoffman EA. Characteristics of the turbulent laryngeal jet and its effect on airflow in the human intra-thoracic airways. Respir Physiol Neurobiol 2007;157:295–309. [22] Smagorinsky J. General circulation experiments with the primitive equations. Mon Weather Rev 1963;91:99–164. [23] Driest ERv. On turbulent flow near a wall. J Aerosp Sci Technol 1956;23:1007–11. [24] Piomelli U, Ferziger JH, Moin P, Kim J. New approximate boundary conditions for large eddy simulations of wall-bounded flow. Phys Fluids 1989;1:1061–8. [25] Pham MV, Plourde F, Doan KS. Direct and large-eddy simulations of a pure thermal plume. Phys Fluids 2007:19. [26] Wang HY, Coutin M, Most JM. Large-eddy simulation of buoyancy-driven fire propagation behind a pyrolysis zone along a vertical wall. Fire Safety J 2002;37:259–85. [27] Dusing M, Kempf A, Flemming F, Sadiki A, Janicka J. Combustion LES for premixed and diffusion flames. Prog Comput Fluid Dynam 2005;5:363–74. [28] Cheung SCP, Yeoh GH, Cheung ALK, Yuen RKK, Lo SM. Flickering behaviour of turbulent buoyant fires using large-eddy simulation. Numer Heat Transfer Part A – Appl 2007;52:679–712.

116

G.E. Lau et al. / Computers & Fluids 59 (2012) 101–116

[29] Yan ZH, Nilsson EEA. Large eddy simulation of natural convection along a vertical isothermal surface. Heat Mass Transfer 2005;41:1004–13. [30] Miki Y, Fukuda K, Taniguchi N. Large eddy simulation of turbulent natural convection in concentric horizontal annuli. Int J Heat Fluid Flow 1993;14:210–6. [31] Barhaghi DG, Davidson L, Karlsson R. Large-eddy simulation of natural convection boundary layer on a vertical cylinder. Int J Heat Fluid Flow 2006;27:811–20. [32] Sagaut P. Large eddy simulation for incompressible flows. 2nd ed. Springer; 2004. [33] Leisuer M, Metais O. New trends in large-eddy simulations of turbulence. Annu Rev Fluid Mech 1996;28:45–82. [34] Sagaut P. Numerical simulations of separated flows with subgrid-models. Rech Aero 1996;1:51–63. [35] Germano M, Piomelli U, Moin P, Cabot WH. A dynamic subgrid-scale eddyviscosity model. Phys Fluids A 1991:3. [36] Padilla ELM, Silveira-Neto A. Large-eddy simulation of transition to turbulence in natural convection in a horizontal annular cavity. Int J Heat Mass Transfer 2008;51:3656–68. [37] Meneveau C, Lund TS, Cabot WH. A Lagrangian dynamic subgrid-scale model of turbulence. J Fluid Mech 1996;319:353–85. [38] Ghosal S, Lund TS, Moin P, Akselvoll K. A dynamic localization model for largeeddy simulation of turbulent flows. J Fluid Mech 1995;286:229–55. [39] Piomelli U, Liu J. Large-eddy simulation of rotating channel flows using a localized dynamic model. Phys Fluids 1994;7:839–48. [40] Fureby C, Grinstein FF. Monotonically integrated large eddy simulation of free shear flows. AIAA J 1999:37. [41] Vreman AW. An eddy-viscosity subgrid-scale model for turbulent shear flow: algebraic theory and applications. Phys Fluids 2004;16:3670–81. [42] You D, Moin P. A dynamic global-coefficient subgrid-scale model for largeeddy simulation of turbulent scalar transport in complex geometries. Phys Fluids 2009:21. [43] You D, Moin P. A dynamic global-coefficient subgrid-scale eddy-viscosity model for large-eddy simulation in complex geometries. Phys Fluids 2007:19. [44] Park N, Lee S, Lee J, Choi H. A dynamic subgrid-scale eddy viscosity model with a global model coefficient. Phys Fluids 2006:18. [45] Vreman AW, Bastiaans RJM, Geurts BJ. A similarity subgrid model for premixed turbulent combustion. Flow Turbul Combust 2009;82:233–48. [46] Peng S, Davidson L. Large eddy simulation for turbulent buoyant flows induced by differentially heated vertical walls. Doktorsavh Chalmers Tek Hogsk 1998;1391:1–18. [47] Peng S-H, Davidson L. Large eddy simulation for turbulent buoyant flow in a confined cavity. Int J Heat Fluid Flow 2001;22:323–31. [48] Eidson T. Numerical simulation of the turbulent Rayleigh–Benard problem using subgrid modelling. J Fluid Mech 1985;158:245–68. [49] Zhang W, Chen Q. Large eddy simulation of indoor airflow with a filtered dynamic subgrid scale model. Int J Heat Mass Transfer 2000;43:3219–31. [50] Zhang W, Chen Q. Large eddy simulation of natural and mixed convection airflow indoors with two simple filtered dynamic subgrid scale models. Numer Heat Transfer Part A 2000;37:447–63.

[51] Barhaghi DG, Davidson L. Natural convection boundary layer in a 5:1 cavity. Phys Fluids 2007:19. doi:10.1063/1.281574. [52] Barhaghi DG, Davidson L. On the validity of the Boussinesq approximation in a developing mixed convection boundary layer. In: Thermal issues in emerging technologies: theory and application, THETA 2007 international conference on 2007; 2007. p. 183–90. [53] Barhaghi DG, Davidson L. Large-eddy simulation of mixed convection– radiation heat transfer in a vertical channel. Int J Heat Mass Transfer 2009;52:3918–28. [54] Lee JS, Pletcher RH. Large eddy simulation of variable property turbulent flow in a vertical channel with buoyancy effects and heat transfer. In: 35th National heat transfer conference, Anaheim, California; 2001. p. 1839–50. [55] Lee JS, Pletcher RH. Large eddy simulations of variable property turbulent flow in horizontal and vertical channels with buoyancy and heat transfer effects. Proc Inst Mech Eng, Part C: J Mech Eng Sci 2007:429–41. [56] Yin J, Wang BC, Bergstrom DJ. Large-eddy simulation of combined forced and natural convection in a vertical plane channel. Int J Heat Mass Transfer 2007;50:3848–61. [57] Peng S, Davidson L. On a subgrid-scale heat flux model for large eddy simulation of turbulent thermal flow. Int J Heat Mass Transfer 2002;45:1393–405. [58] Smith R. Characteristics of turbulence in free convection flow past a vertical plate. Ph.D. thesis. University of London; 1972. [59] Clark RA, Ferziger JH, Reynolds WC. Evaluation of subgrid-scale models using an accurately simulated turbulent flow. J Fluid Mech 1997;91:1–16. [60] Deardorff JW. A numerical study of three-dimensional turbulent channel flow at large Reynolds numbers. J Fluid Mech 1970;41:453–65. [61] Lilly DK. A proposed modification of the Germano subgrid-scale closure method. Phys Fluids 1992;4:633–5. [62] Knio OM, Najm HB, Wyckoff PS. A semi-implicit numerical scheme for reacting flow: II. Stiff, operator-split formulation. J Comput Phys 1999;154:428–67. [63] Miyamoto M, Katoh Y, Kurima J, Sasaki H. Turbulent free convection heat transfer from vertical parallel plates. In: Tien CI, Carey VP, Ferrel JK, editors. Proceedings of the eighth international heat transfer conference, Hemisphere, Washington, DC; 1986. p. 1593–8. [64] Gresho PM. Some current CFD issues relevant to the incompressible NavierStokes equations. Comput Methods Appl Mech Eng 1991;87:201–52. [65] Gresho PM. Incompressible fluid dynamics: some fundamental formulation issues. Annu Rev Fluid Mech 1991;23:413–53. [66] Gebhart B. Instability, transition, and turbulence in buoyancy-induced flows. Annu Rev Fluid Mech 1973;5:213–46. [67] Armfield SW, Patterson JC, Lin W. Scaling investigation of the natural convection boundary layer on an evenly heated plate. Int J Heat Mass Transfer 2007;50:1592–602. [68] Miyamoto M, Kajino H, Kurima J, Takanami I. Development of turbulence characteristics in a vertical free convection boundary layer. In: Grigull U, Hahne E, Stephan K, Straub J, editors. Proc 7th int heat transfer conf. Munchen: Fed. Rep. of Germany; 1982. p. 323–8. [69] Kim J, Moin P. Application of a fractional-step method to incompressible Navier–Stokes equations. J Comput Phys 1985;59:308–23.