Elliptic relaxation model for stably stratified turbulence

International Journal of Heat and Fluid Flow 74 (2018) 173–186

Contents lists available at ScienceDirect

International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff

Elliptic relaxation model for stably stratified turbulence Sandipan Kumar Das

T

ExxonMobil Research and Engineering, 22777 Springwoods Village Parkway, Spring, TX 77389, USA

ARTICLE INFO

ABSTRACT

Keywords: Turbulent Flows Stable Stratification Second Moment Closure (SMC) Models Reynolds Averaged Navier Stokes (RANS) Models Elliptic Relaxation

The article presents a new elliptic relaxation procedure for Second Moment Closure (SMC) turbulence modeling of wall-bounded, stably stratified flows. Specifically, it directly introduces buoyancy effects in the elliptic relaxation equation. The new formulation results in an additional length scale that produces a coupled set of Helmholtz equations. The study then extends the new technique to Reynolds scalar fluxes. Furthermore, it examines the budgets of the different terms of their transport equations near the wall and systematically derives the relevant boundary conditions. Comparison of model predictions with Direct Numerical Simulation (DNS) data for stably stratified plane channel flows available in the literature rounds off this paper.

1. Introduction Turbulent flows are ubiquitous in nature and are almost always the norm in engineering applications. Unsurprisingly, turbulence has historically attracted a lot of attention from fluid dynamicists and engineers alike. The non-linear nature of the governing equations makes it very difficult to develop a practical predictive tool that engineers can readily use for design purposes. Nevertheless, high performance computing has made it possible to numerically simulate the full Navier–Stokes equations for fundamental flows and understand the intricacies of turbulence. This has invigorated the turbulence modeling research community especially in the realms of model development and validation. Turbulence anisotropy, inhomogeneity, wall-effects and stratification are some of the formidable challenges encountered in turbulence modeling. The last two topics are the focus of the current study. The research on modeling the effects of the wall on turbulence has seen quite a bit of progress in the last two decades or so beginning with the seminal work of Durbin (1991) on v 2 f turbulence modeling. That paper addressed the wall blocking effect through an elliptic relaxation method. In the process, Durbin introduced a new rationale for modeling the inhomogeneity of the pressure-strain term. In a later paper (Durbin, 1993), he successfully applied the procedure for a full second moment closure on a wide variety of canonical flows. However, this resulted in a prohibitively large number of equations as there was a separate elliptic equation for each of the Reynolds stresses. On the other hand, the v 2 f turbulence model while adequately and economically considering the wall-blocking effects could not properly capture the turbulence anisotropy as compared to the full Second Moment Closure (SMC). Manceau and Hanjalić (2002) recognized the need for an

economically practical model to effectively address both the wallblocking and anisotropic effects and introduced the so called elliptic blending model that needed just one elliptic relaxation equation for a full Reynolds stress model. The literature has already witnessed the successful application of the model in different applications (Thielen et al., 2005; Fadai-Ghotbi et al., 2008). However, it is important to note that the geophysical fluid dynamics community has never really embraced the above elliptical concept. This is primarily due to two reasons. Firstly, the huge length scales in the order of hundreds or thousands of kilometers in their applications render these models prohibitively expensive. Secondly, the geometries in their computational domains contain walls at most on one side and hence the walleffects generally are not as pronounced as in internal flows encountered in a variety of engineering applications. The above cited research (Thielen et al., 2005; Fadai-Ghotbi et al., 2008) catered mainly to flows without heat-transfer. These models have also been applied to forced convective heat transfer and have generally yielded good agreement with experimental results. Manceau et al. (2000), Behnia et al. (1999) and Ooi et al. (2002) have applied the v 2 f model along with a turbulent diffusivity assumption for turbulent heat flux to generate data that compared reasonably well against experiments. Earlier, Durbin (1993) had successfully applied a scalar flux model for heat transfer to curved wall boundary layer flow. However, the problem becomes more complicated when buoyancy starts playing an active role in the flow. In a review paper, Hanjalić (2002) reiterated the modeling challenges associated with buoyancy like strong anisotropy and the possible existence of multiple regions of stagnant fluid, laminar circulation, transition and fully turbulent regimes in stably stratified flows. Clearly, these physics require more sophisticated treatment. Kenjereš et al. (2005) have turned to quasi-

E-mail address: [email protected]. https://doi.org/10.1016/j.ijheatfluidflow.2018.09.007 Received 1 May 2018; Received in revised form 8 July 2018; Accepted 14 September 2018 Available online 16 October 2018 0142-727X/ © 2018 Elsevier Inc. All rights reserved.

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

linear, algebraic heat flux and buoyancy-extended stress models in this regard. Their model was along the lines of v 2 f model that included a transport equation for the scalar variance. Dehoux et al. (2012) applied the elliptic blending technique to obtain an algebraic model for the turbulent heat fluxes. They showed that their model predictions depended significantly on the value of the length scale, Lθ that was used for the elliptic blending for the heat fluxes. The current study endeavors to modify the elliptic relaxation technique for stratified flows. The starting point is Durbin’s formulation of the pressure-velocity correlation in terms of a turbulence length scale. The standard derivation of the elliptic relaxation equations expresses the pressure-velocity correlation ui

ui

p 1 (x) = xj 4

ui (x)

S (x xj

x x

)

p xj

2. Model development This section systematically develops the elliptical relaxation models used in the study. It starts with presenting the general linear form of the Reynolds Stress Model along with the associated equations for the scalar fluxes and variance. The paper then proceeds to develop the general elliptic relaxation model applicable to the above set. The final subsection examines the budget of the different terms of the second moment transport equations near the wall and systematically derives the relevant boundary conditions. 2.1. General linear Reynolds Stress Model

as

d3x

In the Reynolds Stress Model, the equations for the Reynolds Stress, ui uj are

(1)

ui uj

where S is a term containing all the mean and fluctuating velocity gradients. The exact nature of the term can be found in standard texts on turbulence modeling like that by Durbin and Pettersson Reif (2001). The interested reader is referred to that text for further details. Durbin argued that the non-homogeneity of the wall would require S a turbulence length scale for modeling the correlation, ui (x) x (x ) . He then expressed the correlation as Qij (x )exp

(

t

x LI

), where L

I

= Pij

xk

u i uj (gi uj

uk ui uj

+

ui uk

k xk + gj ui ) +

2u

Uj

Ui xk

uj uk

xk

i uj

(2) where t is the time coordinate, xi is the spatial coordinate, Ui is the mean velocity, ui is the turbulent fluctuating velocity, θ is the turbulent

j

x

ui uj

+ Uk

fluctuating temperature, Pij

is the

integral length scale of turbulence. He reasoned that the exponential term represented the inhomogeneity through the length scale, LI. The current paper modifies the argument in the context of stably stratified flows. The author believes that buoyancy would play a significant role in the correlation. Hence an additional length scale is in order. Stable stratification is known to produce internal gravity waves. An argument can be made that these waves would make the correlation behave like a damped oscillatory function instead of an exponentially decaying one. The entire model is built on these two premises. Damped oscillatory behavior has previously been reported in the literature for two-point statistics in the context of Lagrangian fluid particle dispersion in stably stratified turbulence (Das and Durbin, 2005; van Artrijk and Clercx, 2007). Das and Durbin (2005) developed a Lagrangian Stochastic Model (LSM) from a second moment closure turbulence model that predicted oscillatory behavior for vertical particle velocity autocorrelation and vertical dispersion. van Artrijk and Clercx (2007) studied fluid particle dispersion in stably stratified flows through direct numerical simulations and witnessed the same behavior for vertical autocorrelation function. The objective of the current study is to incorporate this behavior in the elliptic relaxation equations. The rationale is that whenever a fluid particle moves in the vertical direction due to turbulence, buoyancy would tend to move it back to its original position. This should manifest as oscillations in the pressure-velocity correlation functions. Though the physical reasoning is valid for stable stratification, the later part of the exercise extends the model for unstably stratified flows. It is pertinent to reiterate that this elliptic relaxation approach is not really practically suited for applications with large length scales. Even mild turbulence inhomogeneity generated in external flows does not warrant this technique. Hence, the atmospheric and oceanographic sciences communities can easily rule out this method of turbulence modeling. However, this method will be appropriate for simulating engineering equipment like heat exchangers, electronics cooling, turbine blades heat transfer etc. The article develops the model in the next section. At first it lays down the general framework for the Reynolds Stress Model (RSM) before proceeding to the elliptic relaxation procedure. A separate subsection is devoted completely to the derivation of the new model in its general form. The focus then shifts to the formulation of an elliptic relaxation equation for scalar fluxes. Thereafter, the paper describes the test problems for the new model and its predictions. The article ends with a few concluding remarks and provides some possible directions for future research.

tensor,

1

ij

(

p uj x i

+

( (

p ui x j

ij

+

ij

2 3

ij

) ) is the redistribution

), ρ is the density, p is the pressure, ϵ

ij

is

the dissipation rate tensor of the Reynolds stress ui uj , δij is the

(

)

Kronecker Delta, k is the turbulent kinetic energy, is its dis2 ii sipation rate, β is the co-efficient of thermal expansion, gi is the acceleration due to gravity and ν is the kinematic viscosity. The above form is quite general and is often the starting point in research studies on second moment closure modeling. The interested reader is referred to suitable texts like those by Durbin and Pettersson Reif (2001) or by Pope (2000) for further details. The principal challenge is in the modeling of the redistribution tensor, Pij . The traditional procedure is to separate it into the slow and the rapid parts. The slow part does not contain mean velocity gradients and is usually modeled as a linear relaxation of the stress anisotropy

(

ui uj

1

)

tensor, bij towards zero, as originally proposed by k 3 ij Rotta (1951). The rapid part explicitly involves the mean velocity gradients, Ui and initially one resorts to quasi-homogeneous conditions 2

xj

to evaluate Pij before making requisite corrections for inhomogeneity. The homogeneous component of the redistribution tensor, Pijh is modeled by

Pijh =

c1

k

cs k

ui uj c5

ij

2 k 3 ij

c2

ij

2 3

ij

2 3

ij

c3

ij

2 3

ij

ij

(3) where

1 2

Uj Ui + xj xi

ij

=

ij

=

ui uk

ij

=

ui uk

ij

= 1 = 2 1 = 2

Uj xk Uk xj

Ui xk U uj uk k xi

uj uk

(gi uj + gj ui ) kk kk

c1, c2, c3, cs and c5 are different modeling constants. The above equation 174

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

includes the slow model by Rotta (1951). It is the term containing the constant, c1. Following Daly and Harlow (1970), the turbulent transport term in Eq. (2)is modeled as

uk ui uj xk

=

CS Tuk ul

xk

They had arrived at a decision to equate c2θ and c2θ to c2 and c3 respectively as a means to simplify their general stochastic model of turbulent dispersion. The value of was also fixed at 1.5 in their calibration studies. Warhaft (2000) in his review article on passive scalars noted that the value of generally laid between 1.0 and 3.0. He also to be nearly constant at 1.5 in the presence on a uniform observed mean shear gradient. It is worthwhile to mention that while the focus of the research of Das and Durbin (2005) was the homogeneous model, the effects of wall inhomogeneity is the subject in the current endeavor.

ui uj xl

where CS is a constant and T is a turbulent time-scale. The turbulent kinetic energy can be readily evaluated by halving the sum of the normal Reynolds stresses. Sometimes, especially when the values of all the normal stresses are not available, the turbulent kinetic energy is obtained from the following

k k + Uk = t xk

+

+

+

xk

k xk

t k

2.2. New elliptic relaxation model As mentioned before, the first step towards the model development is to introduce a buoyancy length scale, LII (say). The following is proposed

(4)

By analogy with the k equation, the ϵ equation is modeled as

+ Uk

t

c 1(

=

xk

+

)

c2

+

T

+

xk

1

LII =

t

(5)

xk

where σk, c 1, c 2 and σϵ are modeling constants and νt is a turbulent diffusivity that will be defined in a later section. In a similar vein, the equations for the Reynolds scalar flux and the scalar variance are ui t

ui xk

+ Uk

= Pi + +

+

i

1 ( 2

)

i

xj

i

ui

xk xj

(ui uk ) +

1 ( + ) 2

2u

t

2

+ Uk

=

xk

2uk

xk

uk 2 + xk

i

ui xj

2 2

gk

i

=

ui uk

(7)

=

gi

xk

uk

=

1 4

Qij (x )

=

1 4

Qij (x )

x x x x

=

c1

k

ui + c2 uk

Ui U + c3 uk k + c4 ui uk + c5 gi xk xi xk

Qij (x )/2 4 x x

2

k

)

d3x

exp

x x LI

cos

x x LII

exp

x x LI

1 x x exp i 2 LII

d3x

d3x exp

1 LI

i x LII

x

d3x +

exp

1 i + x LI LII

x

d3x

(12)

1 . Before progressing forward, a brief physical explanawhere i = tion for the new mathematical modification is in order. The cosine function provides oscillations to the correlation. The above expression assumes the same wavenumber for all the components of the correlation tensor although an argument can be made that it should be different. This article makes that assumption anyway primarily to keep the formulation simple at this point. However it does not mean that stable stratification modifies the correlation in a directionally independent way. Indeed, buoyancy will result in anisotropic changes to Qij, thereby

2

(9)

=

x x LII

4 x x

whereas the homogeneous form of the pressure scalar correlation, Pih and ϵθ are modeled as

Pih

i

Qij (x )/2

=

(8)

S (x xj

x x

+ exp

Ui xk

2

ui (x)

p 1 ui (x) = xj 4

where Pi is the pressure scalar correlation, i is the rate of production of scalar fluxes, i is the production term due to buoyancy, α is the scalar diffusivity, ϵiθ is the scalar dissipation rate which is generally zero and ϵθ is the variance dissipation rate. The production rate terms are given by i

(11)

xk

where CL2 and C 2 are constants that need to be calibrated. It is to be noted that the denominator in the above expression is actually the buoyancy frequency which is well known as the Brunt-Väisälä frequency in the atmospheric sciences and the oceanographic research communities. The choice seems apt as this frequency plays a central role in internal gravity waves during stable stratification. The velocity-pressure gradient correlation can therefore be modeled as

(6) 2

1

max CL2 k 2 , C 2 ( ) 4

(10)

where c1θ, c2θ, c3θ, c4θ, c5θ and are the modeling constants associated with scalar transport. This paper adopts a simple relation for ϵθ in Eq. (10) as suggested by Durbin and Pettersson Reif (2001). There are instances in the literature where researchers like Kenjereš and Hanjalić (2000) have solved a full transport equation for ϵθ. The last term in Eq. (6) which is a part of the molecular transport is assumed to be included in the model for Pih in Eq. (9). The scalar dissipation vector, ϵiθ is also included in that model. Durbin (1993) had advocated this approach for both the terms in his paper. The literature contains quite a few turbulence models that give different values for the constants. The current study study adopts the model developed by Das and Durbin (2005). That model was specifically designed for stratified flows and calibrated against a bunch of experimental and numerical data on stable stratification in the literature. The values of their model constants were also guided by realizability constraints through a set of stochastic differential equations.

inflicting overall directional adjustments on ui

p (x) . xj

The above expresses the velocity-pressure gradient correlation as a

sum of two integrals. Let us split ui

ui

p xj 1

=

ui

p xj 2

=

Qij (x ) / 2 4

x

x

Qij (x ) / 2 4

x

x

exp exp

( (

p xj

1 LI 1 LI

into ui

i LII

+

i LII

)x )x

p xj 1

and ui

x

d3x

x

d3x

p xj 1

such that

(13)

The kernels in the above integrals are Green’s functions for the modified Helmholtz equation, provided LI and LII are constant. In real life applications they are not but this analysis makes that assumption anyway. Therefore, the following equations can be formulated 175

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

2u 2u

( (

p i x j1

p i x j2

) )

1 LI

i 2 LII

ui

1 LI

i 2 LII

p ui x j2

+

p xj 1

conditions can make f1Rij = f2Rij a valid solution throughout the domain.

Qij

=

Substituting it in Eq. (21) results in f1Iij + f2Iij = 0 . Therefore fij = 2f1Rij . It essentially means that though in general f1ij and f2ij are both complex, their sum, fij is always real. The entire exercise dictates that the physical Dirichlet boundary conditions at the wall for fij should be split equally between f1Rij and f2Rij . Furthermore, the paper selects zero as the wall

2 Qij

=

(14)

2

Following Durbin (1991), an intermediate variable fij is introduced. It is again split into two parts f1ij and f2ij such that fij = f1ij + f2ij . Eq. (14)can be recast in terms of f1ij and f2ij as

( (

2f 1ij 2f 2ij

) )

1 LI

i 2 f1ij LII

1 LI

2

i LII

+

( (

=

f2ij =

h i 2 Pij + bij LII 2k s

1 LI 1 LI

) )

boundary values for both f1Iij and f2Iij . This choice would be consistent with the above identified behavior. The entire analysis makes it clear that only Eqs. (16) and (17) need to be solved for elliptic relaxation. One can also reach the above conclusion by only examining the real

+

h i 2 Pij + bij LII 2k s

Pijh + bij ks

exp

exp

the intermediate variable, fij is in fact equal to s for the general ink homogeneous case. Thus the above Helmholtz equations can be utilized to evaluate the redistribution tensor in the SMC turbulence models for stably stratified flows. The equations suggest that f1ij and f2ij will be complex and they need to be dealt with accordingly. Let us first focus on the equation for fij. f1ij can be expressed as a sum of the real and the imaginary parts as follows Pij

f1ij =

+

2f I 1ij

2 I f = LI LII 1ij

h 1 Pij + bij 2 2k s LII

1 R f LII2 1ij

1 LI2

1 LI2

Pijh + bij 1 I 2 R f + f = LI LII 1ij k sLI LII LII2 1ij

1 LI2

1 R 2 I f + f = LI LII 2ij LII2 2ij

2f I 2ij

1 LI2

1 I f LII2 2ij

2 R f = LI LII 2ij

1 LI2

1 LI2

1 LII2

f1Rij

L

1 LI2

1 LII2

f1Iij + f2Iij f1Iij

+

f2Iij

=

(18)

k sLI LII

=

2 LI LII

2

(20)

f2Rij f1Rij

f2Rij

1 4 + 2 2 LII2 LI LII

2

1 LI2

1 LII2

1

f1Rij

x

p xj 2

) sin (

x

1 LII

1 LI

1 LII

1 LI

x

1 LII

p i x j1

are complex conjugate of each other.

1 L 2II

f1Ri

2 fI L I L II 1i

=

2f I 1i

1 L 2I

1 L 2II

f1Ii +

2 fR L I L II 1i

=

I

II

k2

= max CL1 =

max CL

2

,C

1

() 3

1

1 L 2I

Pih

k sL I L II

2f1Ri

1 L 2II

Pih

2k s

(23)

. The length scales, LθI and LθII

1 4

1

k2 , C 2 ( ) 4 gk

(24)

xk

2

are constants whose values are given in

( (

1 LI 1 LI

+

1 LII 1 LII

) )

2

2

f1ij = f2ij =

( (

1 LI 1 LI

+

1 LII 1 LII

) )

2 Pijh + bij 2k s 2 Pijh + bij 2k s

(25)

The above needs to be solved. The boundary conditions for f1ij and f2ij here will be the same as given for f1Rij . The model was initially built for stable stratification but the mathematics holds for unstably stratified flows too. It must be mentioned that the correlation function did not develop from a physical reasoning in this case. So this is given here just as a mathematical extension. A more elaborate analysis on unstable flows is outside the scope of this paper. Therefore, the current study refrains from making any claim on the performance of the above model in unstable stratification. Unstratified turbulence makes fi = f1i + f2i , the equations being

Substituting Eq. (21) in Eq. (20) the following is obtained

1 LI2

1 LI

) and exp ( x x ) cos ( x x ) x x ) sin ( x x ) and are exp ( x ) . Therefore, it is obvious that u

1 L 2I

2f 2ij 1

(21)

2

(

2f 1ij

1 LII2

x

It is worthwhile to mention that for neutral flows the current model 1 yields L = 0 . As a result the real and imaginary parts of f1ij and f2ij II become decoupled and the elliptic equations revert back to the original procedure introduced by Durbin (1991) for neutral flows. It is to be noted that Durbin’s model had s = 1.0 . For unstable stratification, LII becomes imaginary. In this scenario, the equations for f1ij and f2ij do not contain any imaginary parts. Therefore, the additional damping that was introduced for stable stratification naturally vanishes. The equations are now simply as follows

(19)

1 LI2

and comparing them. In parti-

2.3. Neutral and unstable stratification

Pijh + bij

2 fR LI LII 1ij

x

where CL1 , C 1 , CL2 and C Table 3.

h 1 Pij + bij 2k s LII2

Adding Eqs. (17) and (19) results in 2

1 LII

2f R 1i

L

(17)

2 f I + f2Iij = 0 LI LII 1ij

f2Rij

) cos (

x

3

(16)

The above equations with appropriate boundary conditions need to be solved. However, before proceeding towards the boundary conditions, further mathematical analysis needs to be carried out to simplify the above system. Subtracting Eq. (18) from Eq. (16) gives 2

x

for stably stratified flows with fi = are given by

A similar exercise for f2ij would yield 2f R 2ij

p xj 2

A similar analysis for the scalar flux would yield elliptic relaxation equations for the vector variables, fiθ instead of tensors. The equations would be exactly the same excepting the source terms on the right that would have Pijh + bij replaced by Pih instead. The equations are

if1ijI

1 LI2

1 LI

and ui

where f1Rij and f1Iij are the real and imaginary parts. Substituting the above in the Helmholtz equation for f1ij and separating the real and imaginary parts, the following is arrived at 2f R 1ij

(

while the imaginary parts

. The numerator,

Pijh + bij is the quasi-homogeneous form of the redistribution term. So

f1Rij

and ui

cular, the exponential functions should be focused on. The real parts are

where Pijh is the homogeneous part of the redistribution tensor, Pij and s is a constant. From the above equation it is clear that in the absence of inhomogeneity far away from the walls, fij =

p xj 1

and imaginary parts of ui

(15)

f2Rij = 0 (22)

The above equation suggests that an appropriate choice of boundary 176

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

( (

2f 1i 2f 2i

1 L I

+

1 L I

1 L II

) )

2

2

1 L II

f1i

=

f2i

=

( (

1 L I 1 L I

+

1 L II 1 L II

) )

2 Ph i

the wall

2 Ph i

0=

2k s

(26)

2k s

0 = Pi +

Dreeben and Pope (1997) investigated the Reynolds stress budget near the wall when they extended probability density methods to include near-wall modeling while Manceau and Hanjalić (2002) have examined it later in greater detail in the context of developing their elliptic blending model. The literature also contains studies that have applied elliptic blending to scalar fluxes. The papers by Choi and Kim (2008), Shin et al. (2008) and Dehoux et al. (2017) are some of the more prominent ones in this regard. In the context of the current model, the Reynolds stress budget analysis will be very similar to the ones cited above. However, the treatment of the terms in the scalar flux model, especially the molecular diffusion and dissipation terms is slightly different. Unlike the current model, Shin et al. (2008) modeled the last term of Eq. (6) and the dissipation rate separately instead of including them in the model for the redistribution terms. Hence, a more detailed examination of the scalar flux budget becomes necessary. The article will also dwell on the Reynolds stress budget, albeit briefly for completeness. This study follows the standard procedure of expanding the turbulent variables ui, p and θ in Taylor series near the wall as follows

1

+

p xi

1 ( + ) 2

2u

i

y2

)

y

ui

y

i

(30)

1)

+

w

k w1 + s

4

ui

(31)

y

0 y

The above analysis suggests different values of s for different components of fiθ, depending on whether fiθ is normal or tangential to the wall. However, it may not be possible to follow such guidelines in complex geometries where a single component of fiθ may become normal and tangential at different wall boundaries. Therefore, it seems pertinent to sacrifice the correct asymptotic behavior for some components. It is also prudent to give preference to the normal component. 1 Hence, the paper adopts a value of s = 2 and n = 3 for all components of fiθ. A brief discussion on the boundary conditions for the Reynolds stress equations will round off this analysis. A similar procedure would reduce to the following at the walls

(27)

ui 1 + ( y y 2

i

y2

where n = 2 for i = 1, 3 whereas n = 3 for i = 2 and the subscript w refers to the value at the wall. The value of s should be so chosen that the singularity of fiθ at the wall is avoided. This means that s should take the value 1 for i = 2 and the value 0.0 for i = 1, 3. The reason is that k 2 goes as y2 at the wall. For all practical purposes kw can taken as value of k at the grid point nearest to the wall. The above derivation has as2 k sumed w = lim 2 .

+ + + +

( + )

n (n

fi =

Pij

where ai, bi and ci (i = 1, 2, 3, p and θ) are random functions of x, z and t, y is the direction normal to the wall, the subscript i = 2 is again in the y direction whereas i = 1 and 3 are the wall tangential directions. No slip condition at the wall makes a1 = a2 = a3 = 0 . Invoking the continuity equation near the wall renders b2 = 0 . The paper further assumes that the wall suppresses any turbulent temperature fluctuation on it, so a = 0 . Let us begin the analysis with the Reynolds scalar flux equations. A balance of terms near the wall boils down to the following as they contain the lowest powers of y

0=

2u

1 ( + ) 2

The above should ensure the correct asymptotic behavior of the Reyy n , it is easy to arrive at nolds scalar fluxes at the wall. Assuming ui the following wall boundary condition for fiθ

2.4. Boundary conditions

= a + b y + c y2 +

(29)

p / y in the The above yields bp = 2 c2 that is used to evaluate 1/ table. It is evident that the model needs to predict the following near the wall.

One of the advantages of the elliptic procedure is to predict accurate near-wall behavior. This is ensured by correct specification of the boundary conditions at the walls for the variables fij and fiθ. Budgets of the Reynolds stresses and the scalar fluxes close to the wall can elucidate further on this issue.

u1 = a1 + b1 y + c1 y 2 u2 = a2 + b2 y + c2 y 2 u3 = a3 + b3 y + c3 y 2 p = ap + bp y + cp y 2

2u 2 y2

1 p + y

ui uj k

+

2u

i uj

y2

=0

(32)

The above would lead to the following wall boundary conditions

fij =

2

n (n 2

1)

w

k w1 + s

ui uj

(33)

where n = 3, for i = 1, 3 and j = 2, whereas n = 4 for i = 2 and j = 2 . Finally, n = 2 for i = 1, 3 and j = 1, 3 in the above. The exponent s should take the value 1 for i = 1, 3 and j = 2, the value 1.0 for i = 2 and 2 j = 2 and finally the value 0 for i = 1, 3 and j = 1, 3. Previous research in the literature Manceau and Hanjalić (2002) has suggested a value of n = 4 for i = 1, 3 and j = 2 . This paper also adopts this approach and takes s = 1 accordingly. Additionally, s adopts the same value for the tangential components too as anyway n = 2 makes the boundary value of fij equal to 0 for i = 1, 3 and j = 1, 3.

ui y (28)

3. Numerical test cases

The first three terms on the right hand side of the above equation are the pressure-scalar correlation, the scalar dissipation, ϵiθ and a molecular transport term. As mentioned before, here they are collectively modeled as a single term, Pi which is analogous to the redistribution term in the Reynolds Stress equations (Durbin, 1993). It should be noted that the trace of Pi does not go to zero and hence it is not strictly redistributive. The last term in the above equation signifies molecular diffusion. Table 4 shows the values of the leading order terms in the Reynolds scalar flux equations near the wall. The budget clearly shows the balance of the relevant terms. The calculation uses the following form of the Navier–Stokes equation for the turbulent velocity fluctuations near

This section describes the test cases simulated by the above models in this paper. The literature contains detailed numerical studies (Armenio and Sarkar, 2002; García-Villalba and Del Álamo , 2011) on problems relevant to the current work. Armenio and Sarkar (2002) performed large eddy simulations of stably stratified channel flows. More recently, García-Villalba and Del Álamo (2011) re-examined the same problem and determined that the previous computational domain 4 was too small (4πH in the streamwise direction and 3 H in the spanwise direction, where H is the channel half-height) to obtain meaningful data especially at high stratification levels. Then they proceeded to correct it by performing DNS in sufficiently large domains (8πH and 177

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Fig. 1. Geometry of the flow domain under consideration.

3πH in the streamwise and spanwise directions respectively). They further concluded that small domains could induce spurious laminarization effects and that the work by Armenio and Sarkar (2002) likely suffered from this transgression. The flow domain is a plane horizontal channel with translational periodic conditions in the streamwise x direction while the gravity is in the negative y direction. Symmetric conditions are assumed in the spanwise z direction. Fig. 1 illustrates the geometrical details of the problem pursued here. The following subsections will list the equations describing the flow and the thermal physics using the Reynolds stress and the scalar flux model. It is to be noted that the problem is completely specified by the following non-dimensional parameters

Reynolds Number Re Richardson Number Ri

= =

kfuv

uv

kfvv

v2

k

k

+

d duv ( + CS Tv 2) dy dy

+

d dv2 ( + CS Tv 2) dy dy

dU +g u =0 dy

v2

(39)

+ 2g v = 0

(40)

The elliptic relaxation equations for fuv and fvv are needed here. For stable stratification, fuv = 2f1Ruv and fvv = 2f1Rvv . The equations for solving f1Ruv are d2f1Ruv

1 LI2

1 2 LII

f1Ruv

2 fI LI LII 1uv

d2f1Iuv

1 LI2

1 2 LII

f1Iuv

2 fR LI LII 1uv

dy 2 dy 2

u H

+

1 LI2

= =

h + uv / k P 12

1 2 LII

2k

h + uv / k P 12

(41)

kLI LII

where

g H(

bot )

top 2

u

h P12 =

for stable stratification (34)

dU c1 uv + c2 v 2 k dy

c5 gu

(42)

Table 1 Values of the constants used in the Reynolds Stress model.

where Θtop is the top wall temperature, Θbot is the bottom wall temw perature, u is the friction velocity and τw is the shear stress at the wall. It is to be noted that for all simulations in this article, Θtop ≥ Θbot. The equality sign is for neutral stratification.

c1

c2

c3

c5

cs

c1θ

c2θ

c3θ

c4θ

c5θ

1.8

0.6

0.0

1/3

0.0

2.5

0.6

0.0

0.0

1/3

1.5

3.1. Governing equations Table 2 Values of the constants used in the k–ϵ equations and the turbulent diffusivity expression.

The governing equations for the mean streamwise velocity, U and the mean temperature, Θ are

d dy

dU dy

uv =

d dy

d dy

v

1 dP o dx

(35)

=0

(36)

c1

(

dU +g v dy dU

uv dy + g v

+

)

d dy c2

T

+

t k

+

d dy

dk dy

+

=0

t

d dy

CT

c

0.22

6.0

1.6

c

1

2

1.9

σk

σϵ

σθ

1.2

1.5

2.0

Table 3 Values of the constants used in the elliptic relaxation equations.

where u is the fluctuating velocity in the streamwise direction, v is the fluctuating velocity in the vertical direction and ρ is the mean density. The turbulent kinetic energy, k and its dissipation rate, ϵ equations are

uv

CS

CL1

CL2

C

0.2

0.2

80.0

C

1

2

CL1

CL2

C

50.0

0.2

0.2

250.0

C

1

2

8.0

(37)

=0

Table 4 Relevant terms of the Reynolds scalar flux equation near the wall as functions of y to the leading order.

(38)

where the turbulent viscosity, νt is calculated as CS v 2T .

i

3.2. Reynolds stress model

1

O(y)

( + ) b1 b

O(y)

3

2 b c2 y O(y)

2( + ) c2 b y

( O(y)

The equations for the relevant Reynolds stresses, uv and follows

v2

2

are as 178

1

p xi

( + )

ui y y

( + ) b3 b

1 ( 2

)

y

ui

y

ui y

1 ( 2

+ )

2u i y2

( + ) b1 b

) b c2 y

3( + ) c2 b y

( + ) b3 b

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Fig. 2. Dimensionless mean velocity profile across half channel at Re = 550 and Ri = 0 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 4. Dimensionless mean velocity profile across half channel at Re = 550 and Ri = 120 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

The expressions for the length scales, LI and LII are

3.3. Scalar flux model

1

3

LI = maxCL1

k2

,C

1

1

LII =

The required scalar fluxes, u and v and the variance, obtained by solving

3 4

1

maxCL2 k 2 , C 2 ( ) 4

1

k 2 fu

d g dy

dy 2

d2f1Ivv dy 2

1 LI2 1 LI2

1 2 LII 1 2 LII

f1Rvv f1Ivv +

1

k 2 fv + g

2 fI LI LII 1vv 2 fR LI LII 1vv

1 LI2

= =

h + P 22

k

(v

h + P 22

1 2 LII 2

2 k 3

kLI LII

k

(v

2

2 k 3

)

2v

2k

)

h P22 =

(44)

d2f1Ru dy 2

c1

k

v2

v

2 dU c2 uv 3 dy

4 c5 gv 3

2

v2

dU d + dy dy

d d + dy dy

d d + dy dy

+

t

Prt

+ 2 + 2

+

d 2 dy

dv dy

t

+

Prt 2

k

du dy

t

Prt

=0

=0

can be

(46)

(47)

= 0

(48)

As before, fuθ and fvθ need equations for their solution. For stable stratification, fu = 2f1Ru and fv = 2f1Rv . The equations are

where

2 k 3

d dy

(43)

where CL1 and C 1 are constants from the original elliptic relaxation model of Durbin (1991). The equations for f1Rvv are d2f1Rvv

uv

2

d2f1Iu

(45)

dy 2

Fig. 3. Dimensionless mean velocity profile across half channel at Re = 550 and Ri = 60 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

1 L 2I 1 L 2I

1 L 2II 1 L 2II

f1Ru f1Iu +

2 fI L I L II 1u 2 fR L I L II 1u

1 L 2I

= =

1

1 L 2II

P 1h

1

2k 2

P 1h

k 2 L I L II

(49)

Fig. 5. Dimensionless mean velocity profile across half channel at Re = 550 and Ri = 480 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model. 179

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Fig. 9. Dimensionless mean temperature profile across half channel at Re = 550 and Ri = 120 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dashdot, current model without temperature relaxation; dotted, standard model.

Fig. 6. Dimensionless mean velocity profile across half channel at Re = 550 and Ri = 960 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

where

P1h =

c1

k

u + c2 v

dU dy

(50)

and d2f1Rv

1 L 2I

dy 2

d2f1Iv dy 2

1 L 2II

1 L 2I

1 L 2II

2 fI L I L II 1v

f1Rv f1Iv +

2 fR L I L II 1v

1 L 2I

= =

1 L 2II

P2h

1

2k 2

P 2h 1 k 2 LI LII

(51)

where

P2h =

c1

k

v

c5 g

2

(52)

Finally, the following scalar variance equation is solved Fig. 7. Dimensionless mean temperature profile across half channel at Re = 550 and Ri = 0 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dashdot, current model without temperature relaxation; dotted, standard model.

2v

Fig. 8. Dimensionless mean temperature profile across half channel at Re = 550 and Ri = 60 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dashdot, current model without temperature relaxation; dotted, standard model.

d d + dy dy

+

t

d 2 dy

2

k

=0

(53)

Fig. 10. Dimensionless mean temperature profile across half channel at Re = 550 and Ri = 480 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dashdot, current model without temperature relaxation; dotted, standard model. 180

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

3.5. Numerics It is well documented (Durbin and Pettersson Reif, 2001) that turbulence model equations are in general strongly coupled together and their boundary conditions can make them particularly stiff to solve numerically. A plethora of solution techniques have evolved over the years to circumvent the numerical issue. Solving the equations in a coupled manner along with backward implicit transient techniques are some of the more popular approaches that have been practiced in the literature (Deshpande and Giddens, 1977; Deshpande, 1984; Das, 2017). Deshpande and Giddens (1977) developed an algorithm for directly inverting a tridiagonal matrix with diagonal submatrices on either side in the context of second order spatial discretization of the k turbulence model as a coupled system. Later Deshpande (1984) extended the algorithm for an arbitrary number of diagonal submatrices. Recently, Das (2017) developed a numerical algorithm for directly inverting a pentadiagonal matrix containing an arbitrary number of tridiagonal submatrices. He showed that a fourth order coupled spatial discretization of the turbulence models like the k and the v 2 f would result in such types of matrices and later proceeded to solve them using his technique in that paper. The current study adheres to second order spatial discretization of the above equations with a control volume approach. The k and the ϵ equations, the f1Ruv and the f1Iuv equations and the f1Rvv and the f1Ivv equations are grouped together as three equation sets. Each set is then solved individually as a coupled system. As a further aid to convergence, a temporal derivative term is added to all the governing differential equations described in the section except those for f1Ruv , f1Iuv , f1Rvv and f1Ivv while a backward implicit temporal discretization is deployed. Convergence is declared when the solution reaches a steady state. All the simulations employ a uniform mesh of 551 points. This places the nearest point away from the wall at a distance of y+ = 1. As part of the mesh convergence study, additional simulations for Ri = 960 were conducted with mesh sizes of 651, 451 and 251 points. These translate to y+ of 0.85, 1.22 and 2.20 respectively. There was not much difference in the results. In fact, the predictions from the finest three meshes (corresponding to y+ of 0.85, 1 and 1.22) were almost indistinguishable.

Fig. 11. Dimensionless mean temperature profile across half channel at Re = 550 and Ri = 960 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dashdot, current model without temperature relaxation; dotted, standard model.

3.4. Boundary conditions The computational domain consists of half of the horizontal channel starting from the bottom wall to the channel centerline. No-slip walls ensure 0 value for U, k, uv, v 2, u and v . 2 also takes 0 value at the wall. Dirichlet boundary condition is given for the temperature at the 2k wall. The value of w is set at 2w , where kw and y2 are the values of k y2

and the y coordinate at the mesh point nearest to the wall. The values of 5 w 2 1 + w w u w v , 4 f1Ruu , f1Ruv , f1Ru and f1Rv are specified as 3 uvw , 2 2 w k

and

3 4

+

w 3 kw2

v

w

where uvw, vw2 , u

kw

kw2

w

and v

w

w

are the values at the grid

point closest to the wall. At the centerline, the temperature value is set at 0.5( bot + top) while uv takes the value 0. The values of f1Iuu , f1Iuv , f1Iu and f1Iv are each set at 0 at both the walls and the centerline. The first order derivatives of the rest of the quantities are set to 0 at the centerline.

Fig. 12. Dimensionless uv profile across half channel at Re = 550 and Ri = 0 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

181

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Fig. 13. Dimensionless uv profile across half channel at Re = 550 and Ri = 60 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 15. Dimensionless uv profile across half channel at Re = 550 and Ri = 480 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

4. Results and discussion The paper tests three different models. The first is the model with the newly developed elliptic relaxation procedure for both the sets of the Reynolds stress and the thermal flux equations. The second is the model with the new elliptic relaxation equations for the Reynolds stress only. The literature contains arguments both in favor and against the use of elliptic relaxation for the scalar fluxes. Durbin (1993) reported that the additional equations for scalar fluxes did not make any difference to the results and published heat transfer predictions that used elliptic relaxation for only the Reynolds Stresses. On the other hand Ji (2006) observed that some models needed elliptically relaxed scalar fluxes for better accuracy. For example, he advocated its use for the model that Ji and Durbin (2004) had developed earlier in the context of stable stratification. The third model that is used here is the standard elliptic relaxation model that has only one length scale. This model does not enforce elliptic relaxation for the thermal flux equations. It should be noted that the basic Reynolds stress model used in this paper is the same and the constants are given in Tables 1–3. The constants for the Reynolds stresses and the scalar fluxes are those that were used by Das and Durbin (2005). That model compared very well against

homogeneous stably stratified turbulence experiments. Therefore, the current exercise is essentially a test for the new elliptic relaxation technique. The study presents model predictions for stably stratified flows at Re = 550 and Ri = 0, 60, 120, 480 and 960. The predictions are compared against the DNS data of García-Villalba and Del Álamo (2011). This dataset has already found usage in validation studies of turbulence models by Karimpour and Venayagamoorthy (2014). They derived a simple two equation k model for stably stratified flows and compared their mean velocity and mean temperature predictions against the DNS data for Ri = 60 and 120. They did not test their model at higher stratification levels. The main focus of this paper is also on the mean velocity and the temperature profiles across the channel. The important second moments that contribute to these mean profiles are the turbulent momentum and the buoyant thermal fluxes, uv and v . These will also be compared in detail against the DNS data. Finally, the section will exhibit a brief comparative study of the normal stresses and the thermal variance. The Richardson Number, Riτ is gradually increased from 0 to 960 in all the plots. As mentioned before, the Reynolds Number, Reτ is kept constant at 550 for all the simulations reported in the paper.

Fig. 14. Dimensionless uv profile across half channel at Re = 550 and Ri = 120 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 16. Dimensionless uv profile across half channel at Re = 550 and Ri = 960 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model. 182

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Fig. 19. Dimensionless v profile across half channel at Re = 550 and Ri = 120 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 17. Dimensionless v profile across half channel at Re = 550 and Ri = 0 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Figs. 2–6 show the mean velocity profiles for the different Richardson Numbers. The models do not show any difference in predictions for neutral flows as evidenced in Fig. 2. It is expected as the new model reduces to the standard model at Ri = 0 . As the Richardson Number is increased, the higher stratification level alleviates turbulence and the resulting laminarization effect manifests in the velocity profiles tending to move towards a more parabolic shape. All the models capture this trend, however their predictions show variation amongst themselves, especially at higher Richardson Numbers Ri = 480 and 960. Figs. 5 and 6 capture this discrepancy. The standard model significantly underpredicts the centerline velocity whereas the new model without the elliptic relaxation of the thermal flux does quite better. The predictions of the new model with elliptically relaxed thermal flux match extremely well with the DNS data. Figs. 7–11 present the mean temperature profiles across the channel at different Richardson Numbers. The profiles at neutral stratification in Fig. 7 show some variation amongst the results of the models with and without elliptical relaxation of the thermal flux. The predictions of the former are better throughout the domain. The laminarization effects at increasing Richardson number change the sign of the curvature of the temperature profile as highlighted in Figs. 8–11. This shift is captured

Fig. 20. Dimensionless v profile across half channel at Re = 550 and Ri = 480 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 18. Dimensionless v profile across half channel at Re = 550 and Ri = 60 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model.

Fig. 21. Dimensionless v profile across half channel at Re = 550 and Ri = 960 . Circles: DNS data of García-Villalba and Del Álamo (2011). Lines: solid, current model with temperature relaxation; dash-dot, current model without temperature relaxation; dotted, standard model. 183

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

Figs. 17–21 represent the predictions of the turbulent buoyant flux, v . A couple of points need to be indicated. Except for the neutral stratification case in Fig. 17, none of the models predict the v as well as the prediction of the turbulent momentum flux that was shown before. High stratification level reduces the turbulent flux and makes it go towards zero near the centerline in Figs. 20 and 21. Though the models seem to qualitatively agree with this phenomenon, not a single model actually predicts zero value at the centerline. Even then the new models show more obviation of turbulence than the standard model. Furthermore, the new model with the elliptic relaxation of the thermal flux seems to do a better job near the wall. The overall nature of the v profiles are reasonably well represented at all Richardson Numbers. The final segment of this section examines the behavior of the v2 ) normal stresses, u2 and v 2 as profiles of urms ( u2 ) and vrms ( across the channel in Figs. 22 and 23. Furthermore, they are non-dimensionalized by the channel bulk mean velocity, Ub. This shows the trends more clearly as the stratification is varied. Fig. 24 explores the 2 non-dimensioprofiles of the thermal variance, 2 through rms ( top ) nalized by . Each of these figures contains four plots. bot Part (a) plots the DNS data. The rest three are model predictions. Parts (b) and (c) refer to the predictions from the current model with and without elliptic relaxation for the thermal flux respectively while part (d) displays the standard model results. All the models seem to capture the overall nature of the DNS results. Fig. 22 shows a peak for urms near the wall before it decays towards the centerline. The peaks reduce in height at increasing stratification levels. The models predict less decay at higher Riτ compared to the DNS data. Figure 23 suggests good agreement of the model predictions for the vrms profiles. The peak values of the stress shifts away from the wall as compared to urms, all the models predict this shift quite well. Stratification suppresses the vrms values throughout the channel. Though all the models display this trend, these curves show more clustering towards each other when compared against the DNS data. However, parts (b) and (c) evince that the current model both with and without elliptically relaxed thermal flux does a better job in this aspect. Fig. 24 illustrates the θrms profiles. The DNS data shows a small localized peak near the wall, especially for the neutral flow. The profile then displays higher values towards the

Fig. 22. Dimensionless urms profile across half channel at Re = 550 for Ri = 0, 60, 120, 480 and 960. (a) DNS data of García-Villalba and Del Álamo (2011). (b) Results from current model with elliptic relaxation for thermal flux. (c) Results from current model without elliptic relaxation for thermal flux. (d) Results from standard model. Lines: thick solid, Ri = 0 ; long dash, Ri = 60 ; thin solid, Ri = 120 ; short dash, Ri = 480 ; dash-dot, Ri = 960 .

by all the models. The new models predict better than the standard model at Ri = 60 and 120 with the elliptic relaxation of the thermal flux showing closer match with the DNS results. Figs. 8 and 9 demonstrate this aspect. Figs. 10 and 11 inform us that there is not much to choose between the performances of the models at very high stratification levels. Figs. 12–16 exhibit the comparison of the turbulent momentum flux, uv . All the models agree very well with the DNS data. The laminarization evinces itself very subtly towards the centerline of the channel at high Riτ. Figs. 15 and 16 show the DNS uv values plateauing to zero just before the centerline. The new models capture this feature better than the standard model.

Fig. 23. Dimensionless vrms profile across half channel at Re = 550 for Ri = 0, 60, 120, 480 and 960. (a) DNS data of García-Villalba and Del Álamo (2011). (b) Results from current model with elliptic relaxation for thermal flux. (c) Results from current model without elliptic relaxation for thermal flux. (d) Results from standard model. Lines: thick solid, Ri = 0 ; long dash, Ri = 60 ; thin solid, Ri = 120 ; short dash, Ri = 480 ; dash-dot, Ri = 960 .

184

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das

kind of analysis can indeed be insightful for the current model but it is outside the scope of this article. Furthermore, the author is of the opinion that further research is needed before the model can be made ready for practical engineering applications. The current model in its most general form would contain a total of eighteen Helmholtz equations for the wall treatment of the Reynolds stresses and the scalar fluxes. Hence, the obvious next step is to simplify the model to reduce the number of equations. Previous research on elliptic blending by Manceau and Hanjalić (2002) or the v 2 f modeling by Durbin (1991) provides possible future approaches towards simplification. Extension of the simplified model to unstably stratified flows is also an important pre-requisite before it can be confidently applied to more complex flows. Acknowledgment The author gratefully acknowledges the financial support provided by ExxonMobil. Fig. 24. Dimensionless θrms profile across half channel at Re = 550 for Ri = 0, 60, 120, 480 and 960. (a) DNS data of García-Villalba and Del Álamo (2011). (b) Results from current model with elliptic relaxation for thermal flux. (c) Results from current model without elliptic relaxation for thermal flux. (d) Results from standard model. Lines: thick solid, Ri = 0 ; long dash, Ri = 60 ; thin solid, Ri = 120 ; short dash, Ri = 480 ; dash-dot, Ri = 960 .

Supplementary material Supplementary material associated with this article can be found, in the online version, at 10.1016/j.ijheatfluidflow.2018.09.007 . References

channel centerline. It also shows strange behavior at increasing Riτ. Highest centerline value is noted at Ri = 60, then the value at the centerline starts decreasing as Riτ increases further. It should be noted that the localized peaks near the wall show a monotonic decrease as Riτ increases from zero. All the models capture these phenomena. The models without the elliptic relaxation of the thermal flux slightly overpredict the peak near the wall. The DNS data shows faster increase in value near the center whereas the models predict a gradual rise though the centerline values are more or less similar. The above exercise unequivocally demonstrates the better predicting capability of the newly formulated elliptic relaxation technique at high Richardson Numbers of 480 and 960 for the key flow and thermal statistics. At low levels of stratification there is not too much of a difference in the performance as compared to the standard model. As mentioned before, it is not surprising as the new model reduces exactly to the standard model at zero stratification.

Armenio, V., Sarkar, S., 2002. An investigation of stably stratified turbulent channel using large-eddy simulation. J. Fluid Mech. 459, 1–42. van Artrijk, M., Clercx, H., 2007. The effect of stable stratification on fluid particle dispersion. Particle-Laden Flow. Springer, pp. 151–163. Behnia, M., Parneix, S., Shabany, U., Durbin, P., 1999. Numerical study of turbulent heat transfer in confined and unconfined impinging jets. Int. J. Heat Fluid Flow 20, 1–9. Choi, S.-K., Kim, S.-O., 2008. Treatment of turbulent heat fluxes with the elliptic-blending second-moment closure for turbulent natural convection flows. Int. J. Heat Mass Transf. 51, 2377–2388. Daly, B., Harlow, F., 1970. Transport equations of turbulence. Phys. Fluids 13, 2634–2649. Das, S., 2017. Direct solver for pentadiagonal matrix containing tridiagonal submatrices. Numer. Heat Transf. 72, 1–20. Das, S., Durbin, P., 2005. A Lagrangian stochastic model for dispersion in stratified turbulence. Phys. Fluids 17. 025109 Dehoux, F., Benhamadouche, S., Manceau, R., 2017. An elliptic blending differenetial flux model for natural, mixed and forced convection. Int. J. Heat Fluid Flow 63, 190–204. Deshpande, M., 1984. Non-iterative solution of tridiagonal-type matrix coupled by diagonal submatrices. Int. J. Numer. Methods Eng. 20, 297–304. Deshpande, M., Giddens, D., 1977. Direct solution of two linear systems of equations forming coupled tridiagonal-type matrices. Int. J. Numer. Methods Eng. 11, 1049–1052. Dreeben, T., Pope, S., 1997. Probability density function and Reynolds-stress modeling of near-wall turbulent flows. Phys. Fluids 9, 154–163. Durbin, P., 1991. Near-wall turbulence closure modeling without damping functions. Theor. Comput. Fluid Dyn. 3, 1–13. Durbin, P., 1993. A Reynolds stress model for near wall turbulence. J. Fluid Mech. 249, 465–498. Durbin, P., Pettersson Reif, B., 2001. Statistical Theory and Modeling for Turbulent Flows. Wiley, West Sussex, England. F., Dehoux, Lecocq, Y., Benhamadouche, S., Manceau, R., Brizz, L.-E., 2012. Algebraic modeling of the turbulent heat fluxes using the elliptic blending approach – application to forced and mixed convection regimes. Flow Turbul. Combust. 88, 77–100. Fadai-Ghotbi, A., Manceau, R., Borée, J., 2008. Revisiting URANS computations of the backward-facing step flow using second moment closures. influence of the numerics. Flow Turbul. Combust. 88, 395–414. García-Villalba, M., Del Álamo, J.C., 2011. Turbulence modification by stable stratification in channel flow. Phys. Fluids 23. 045104 Hanjalić, K., 2002. One-point closure models for buoyancy-driven turbulent flows. Annu. Rev. Fluid Mech. 34, 321–347. Ji, M., 2006. Second Moment Closure Modeling for Rotating Stably Stratified Turbulent Shear Flow. Stanford Univeristy Ph.D. thesis. Ji, M., Durbin, P., 2004. On the equilibrium states predicted by second moment models in rotating, stably stratified homogenous shear flow. Phys. Fluids 16, 3540–3556. Karimpour, F., Venayagamoorthy, S., 2014. A simple turbulence model for stably stratified wall-bounded flows. J. Geophys. Res. Oceans 119, 870–880. Kenjereš, S., Gunarjo, S., Hanjalić, K., 2005. Contribution to elliptic relaxation modelling of turbulent natural and mixed convection. Int. J. Heat Fluid Flow 26, 569–586. Kenjereš, S., Hanjalić, K., 2000. Convective rolls and heat transfer in finite length Rayleigh-Benard convection: a two-dimensional study. Phys. Rev. E 62, 7987–7998.

5. Conclusion This article tackles the challenges associated with RANS modeling of stably stratified wall-bounded flows. The new approach introduces a buoyancy related length scale in the elliptical relaxation procedure. It results in a coupled set of twin Helmholtz equations. The physical basis for incorporating the buoyancy length scale is valid for stable stratification, so the model is tested for only that class of flows. Neutral flows reduce the elliptical relaxation model to the standard elliptical relaxation model. The budgets of the different terms in the transport equations of the Reynolds stresses and the scalar fluxes yield the boundary conditions that ensure the correct near-wall behavior of the corresponding second moments. The model compares quite well against the DNS results of García-Villalba and Del Álamo (2011), especially the mean velocity and temperature at high stratification levels. In the past, researchers like Dehoux et al. (2017) have performed extremely detailed a priori tests of their model with the help of available DNS data before validating it against more complicated flows. Specifically, they modified the elliptic blending scalar flux model of Shin et al. (2008) by introducing new length and time scales and proceeded to successfully check their hypothesis by a thorough comparison of different redistribution terms against previous DNS studies. Such 185

International Journal of Heat and Fluid Flow 74 (2018) 173–186

S.K. Das Manceau, R., Hanjalić, K., 2002. Elliptic blending model: a new near-wall Reynolds-stress turbulence closure. Phys. Fluids 14, 744–754. Manceau, R., Parneix, S., Laurence, D., 2000. Turbulent heat transfer predictions using v 2 f model on unstructured meshes. Int. J. Heat Fluid Flow 21, 320–328. Ooi, A., Iaccarino, G., Durbin, P., Behnia, M., 2002. Reynolds averaged simulation of flow and heat transfer in ribbed ducts. Int. J. Heat Fluid Flow 23, 750–757. Pope, S., 2000. Turbulent Flows. Cambridge, Cambridge, United Kingdom. Rotta, J., 1951. Statistische theorie nichthomogener turbulenz. Z. Phys. 129, 547–572.

Shin, J., An, J., Choi, Y., Kim, Y., M.S., K., 2008. Elliptic relaxation second moment closure for the turbulent heat fluxes. J. Turbul. 9, 1–29. Thielen, L., Hanjalić, K., Jonker, H., Manceau, R., 2005. Predictions of flow and heat transfer in multiple impinging jets with an elliptic-blending second moment closure. Int. J. Heat Mass Transfer 48, 1583–1598. Warhaft, Z., 2000. Passive scalars in turbulent flows. Annu. Rev. Fluid Mech. 32, 203–240.

186