RANS Modeling of Stably Stratified Turbulent Boundary Layer Flows in OpenFOAM

Quantifying mixing processes relating to the transport of heat, momentum, and scalar quantities of stably stratified turbulent geophysical flows remains a substantial task. In a stably stratified flow, such as the stable atmospheric boundary layer (SABL), buoyancy forces have a significant impact on the flow characteristics. This study investigates constant and stability-dependent turbulent Prandtl number (Prt) formulations linking the turbulent viscosity (νt) and diffusivity (κt) for modeling applications of boundary layer flows. Numerical simulations of plane Couette flow and pressure-driven channel flow are performed using the Reynolds-averaged Navier-Stokes (RANS) framework with the standard k-ε turbulence model. Results are compared with DNS data to evaluate model efficacy for predicting mean velocity and density fields. In channel flow simulations, a Prandtl number formulation for wall-bounded flows is introduced to alleviate overmixing of the mean density field. This research reveals that appropriate specification of Prt can improve predictions of stably stratified turbulent boundary layer flows.


Introduction
Most environmental and geophysical flows are characterized by turbulence and stratification such as the atmosphere, oceans, lakes, estuaries, and rivers.A substantial feature of turbulence is that it enhances mixing.Stable stratification generally suppresses mixing as turbulent scales are reduced and kinetic energy is transferred to potential energy.It is essential to understand the mixing processes which influence the transport of momentum, heat, and scalars (i.e., pollutants, sediments, and nutrients) in boundary layer flows such as the atmospheric boundary layer (ABL).These transport processes are very active in the ABL and have direct implications on wind energy, weather, climate change, and air quality.Buoyancy forces due to density (or analogously potential temperature) variations play an active role in mixing influencing the overall flow structure in boundary layers.Mixing of momentum and density are quantified through the eddy (or turbulent) viscosity (K M ) and eddy (or turbulent) diffusivity (K H ) defined as where u w is the turbulent momentum flux, u is the mean streamwise velocity, z is the height, ρ w is the density flux, and ρ is the mean density.The turbulent Prandtl number a Corresponding author: wilsonjm@lipscomb.edub Current affiliation: Department of Mechanical Engineering, Lipscomb University, 1 University Park Drive, Nashville, TN 37204, USA (Pr t ) relates viscous to diffusive mixing by Pr t = K M /K H . Pr t is found to be on the order of unity for simple turbulent flows [1].While Pr t is generally assumed to be a constant in the range 0.5 − 1.0 in numerical modeling of stably stratified turbulent flows, observational studies [2][3][4], direct numerical simulation (DNS) [5][6][7], and analytical models [8][9][10] have shown that the value of Pr t is strongly linked to stability.The research presented in this paper seeks to investigate several Pr t propositions for efficacy in numerical simulations of stably stratified turbulent flows.

Stability Parameterization
In assessing the level of stratification an appropriate characterization or parameterization is required.The gradient Richardson number (Ri g ) is an important stability parameter given by where N = [(−g/ρ 0 )(dρ/dz)] 1/2 is the Brunt-Väisälä (or buoyancy) frequency, g is the gravitational acceleration, ρ 0 is the reference density, and S = du/dz is the mean shear.For stable stratification in a homogeneous shear flow, N > 0 resulting in positive values for the gradient Richardson number (Ri g > 0).The weakly stable regime (Ri g 1) represents a state of continuous turbulence.The strongly stable regime (Ri g 1) exhibits intermittent turbulence with suppressed scales of mixing and the formation of gravity waves or two-dimensional turbulence through the transfer of kinetic to potential energy [11].Recently, Flores & Riley [12] used DNS to suggest that this regime change can occur due to a lack of separation of turbulent and viscous scales.Notwithstanding, the delineation between these regimes, termed the critical gradient Richardson number (Ri gc ), is generally accepted to be near 0.20 − 0.25 (see e.g., [13]).
Another closely related parameter is the flux Richardson number (Ri f ) defined as the ratio of buoyancy production (B) to production of turbulent kinetic energy (P).
where B = −(g/ρ 0 )ρ w is the buoyancy production and P = −u w (du/dz) is the production of turbulent kinetic energy due to shear.The flux Richardson number (Ri f ) lends greater clarity on the state of flow with continuous turbulence for Ri f ≤ Ri f ∞ .While the exact value of Ri f ∞ remains an area of active research, Ri f ∞ = 0.20 − 0.25 is consistent with oceanic and atmospheric observations [14,15].For surface-bounded flows, dimensionless scaling through Monin-Obukhov similarity theory [16] or Nieuwstadt's local scaling [17] may also be used to quantify stability.The Obukhov length is used to measure the strength of stratification based on surface-level measurements.
where u * = [−u w ] 1/2 is the friction velocity defined at the surface, κ is the von Kármán constant, and (ρ w ) 0 is the surface density flux.If the variables in Eq. 5 are defined through local measurements the definition is known as the local Obukhov length (Λ).However, the flux Richardson number (Ri f ) is found to be a parameter of greater significance that can be used to assess surface-bounded and freeshear flows [18].The gradient Richardson number (Ri g ) is also valuable to our understanding of stably stratified turbulence because of its relevance to field observations and measurements.

The Turbulent Prandtl Number
In the context of Reynolds-averaged Navier Stokes (RANS) modeling, the turbulent Prandtl number (Pr t ) links the turbulent viscosity and diffusivity with Pr t = K M /K H . Turbulence modeling within the context of stably stratified flows generally takes two forms: (i) Defining a strong stably limit (Ri gc ) above which turbulence cannot exist [19,20]; or (ii) Assuming that turbulence may exist for all Ri g , albeit intermittent or decaying at large Ri g [11].
The former approach rooted on the concept of a definitive critical gradient Richardson number (Ri gc ).However, recent studies have shown evidence of turbulence at large Ri g maintained by strong shear [21].Taking the latter approach, two stability-dependent turbulent Prandtl number (Pr t ) propositions are selected for this study.
Firstly, Kim & Mahrt [4] (hereafter KM92) used atmospheric data from aircraft measurements of the SABL to propose a Pr t based on the theoretical formulation of Louis [22] as a function of Ri g given by where Pr t0 is the turbulent Prandtl number at neutral stability (Ri g = 0).For the weakly stable regime, Pr t ≈ Pr t0 + 4Ri g and transitions to scale as Pr t ≈ 5Ri g for strong stability with Ri f ∞ = 0.2 in the limit Ri g → ∞.A myriad of values for Pr t0 can be found in literature.The majority of these assignments are from observational measurements such as Businger's [23] value of Pr t = 0.74.More recently, Venayagamoorthy & Stretch [9] used relevant time and length scales from first principles to construct a Pr t0 for homogeneous turbulent flows.
where γ = (1/2)(T L /T ρ ) is the ratio of relevant time scales in turbulence models [1], T L = k/ε is the turbulent kinetic energy decay time scale, 2 is the scalar dissipation, κ is the molecular diffusion rate, L M = (q 2 1/2 )/S is a characteristic mixing length, q 2 = 2k is twice the turbulent kinetic energy (k), and L E = (ρ 2 ) 1/2 /|(dρ/dz)| is the Ellison length scale [24].The ratio L M /L E is found to be a strong function of stability that approaches unity in the limit Ri g → 0 with γ ≈ 0.7 for neutral stratification yielding Pr t0 = 0.7.
Additionally, the analytical Pr t model of Venayagamoorthy & Stretch [9] (hereafter VS10) is considered.The authors provided a theoretical modification to the empirical formulation of Schumann & Gerz [8] based on DNS data [5,7].This model provides insights for mixing of homogeneous stably stratified turbulence being applicable to stationary and nonstationary flows. where represents the the mixing efficiency in the limit Ri g → ∞.Venayagamoorthy & Stretch [9] found that Γ ∞ = 1/3 fits the data corresponding to Ri f ∞ = 0.25 [7].Pr t0 is taken to be 0.7 from the previous discussion.For the weakly stable regime, Pr t ≈ Pr t0 + Ri g and transitions to scale as Pr t ≈ 4Ri g for strong stability.Equation 8 provides a homogeneous Pr t framework that can be modified for any values of Γ ∞ , Ri f ∞ , and Pr t0 .For example, Osborn [25] suggested for oceanic flows and Weinstock [26] observed Γ ∞ ≈ 0.8 (Ri f ∞ = 4/9) for the upper atmosphere.This study maintains the values of Γ ∞ = 1/3, Ri f ∞ = 0.25, and Pr t0 = 0.7.

Paper Layout
The research presented in this paper investigates the performance of turbulent Prandtl number parameterizations for prediction of turbulent mixing and diffusion in stably stratified turbulent flows.The Reynolds-averaged framework is employed in conjunction with the standard k-ε turbulence model.The relevant governing and transport equations are presented in Section 2. Additionally, the details of the numerical simulations are included for the considered scenarios of plane Couette flow and pressuregradient driven channel flow.In Section 3, simulation results are evaluated with DNS data and thoroughly discussed.Conclusions are made in Section 4.

Theoretical Development
This section provides the theoretical background for numerical simulations.For this study, the RANS framework is implemented.One could question the relevance of RANS models when large-eddy simulation (LES) studies are increasingly attainable with modern computing power.
The basis of any turbulence model, RANS or LES, is to provide a description of the production, interaction, and destruction of turbulence.RANS simulations remain significantly less computationally expensive and have proven accurate for many engineering and geophysical applications.The k-ε turbulence closure scheme is selected for this initial study and provides information which can be used scaling towards SABL interactions with wind turbines, for example.Introducing the k-ω shear stress transport (SST) model is one possibility for extension to waked flows.The model performs well in the near-wall region (implicitly a low-Reynolds number model with sufficient grid resolution) and embeds the desirable free-shear performance of the k-ε model.Ultimately, the end goal of all turbulence models is to properly capture the physics essential to the problem at hand.If this objective can be achieved with a lower-order model (e.g., RANS) further sophistication is unnecessary.The Reynolds-averaged governing equations are presented alongside the two-equation k-ε turbulence closure scheme.Finally, the simulation details are presented along with a discussion on the numerical implementation in the open-source computational fluid dynamic (CFD) code OpenFOAM [27].

Governing Equations
The Reynolds-averaged equations for continuity and conservation of momentum assuming the validity of the Boussinesq approximation are given by Eqns.9 and 10, respectively, using Einstein summation convention.
where U j is the mean velocity, ρ 0 is a reference density, p is the mean pressure field, ν eff = ν + ν t is the effective viscosity, ν is the molecular viscosity, and ν t is the turbulent viscosity introduced by the turbulent viscosity hypothesis.For two-equation turbulence closure, ν t is closed for algebraically as the product of a relevant velocity and length scale (ν t = u * * ).The specification of these scales is dependent on the selection of turbulence model.For a stratified flow, a transport equation for density is required, where ρ is the mean density, κ eff = κ + κ t is the effective diffusivity, κ is the molecular diffusivity, and κ t is the turbulent diffusivity solved for using the gradient diffusion hypothesis and turbulent Prandtl number (Pr t ).It is noted that ν t and κ t are the modeled turbulent viscosity and diffusivity quantities analogous to K M and K H given by Eqns. 1 and 2, respectively.

The Standard k-ε Model
The two-equation k-ε turbulence closure scheme [28] is an energy transport model solving for the turbulent kinetic energy (k) and the dissipation rate of the turbulent kinetic energy (ε).The model is widely used for engineering and geophysical applications [11,29].The modeled transport equations for k and ε are given by where σ k is the turbulent Prandtl number for k, σ ε is the turbulent Prandtl number for ε, and C ε1 , C ε2 , and C ε3 are model constants.The velocity scale for this model is u * = k 1/2 and the length scale is * = C μ k 3/2 /ε.Combining these scales, the turbulent viscosity is given by ) 2 is a model variable typically defined with the engineering constant of 0.09.However, C μ is also speculated to vary with stability (see e.g., [30]) the focus of which is beyond the immediate scope of this study.As with most turbulence closure models, the k-ε model was developed for nonbuoyant flows and appropriate modifications are required for stratification.Assuming local equilibrium in a stratified flow, the dominant balance in Eq. 12 reduces to P + B = ε which can subsequently be shown to moderate C μ by (1−Ri f ) for the turbulent viscosity (ν t ) given in Eq. 14.
Table 1 presents the standard model constants used in this study [31].All coefficients except for C ε3 take on their standard values calibrated for neutral stability.The assignment of C ε3 has been the subject of much research.
Baumert & Peters [32] provided an excellent review of 2nd Symposium on OpenFOAM in Wind Energy values for C 3 .The general procedure for modeling stably stratified flows within the the k-ε model framework is to assign a constant turbulent Prandtl number, Pr t ≈ 0.85, and modify the 'buoyancy production of dissipation' term in Eq. 13.This modification is accomplished by varying C ε3 to calibrate the model.Simulation results were found to be very sensitive to the value of C 3 .A value of C 3 = −1.44,similar to the value of -1.4 found by Burchard & Baumart [33], provided the greatest accuracy compared with DNS data.The assignment of a negative coefficient is physically consistent with stable stratification yielding a negative buoyancy production (or destruction of turbulent kinetic energy) which is balanced by a positive value for the 'buoyant production of dissipation' term.Applying C 3 = −1.44 to all simulation runs, any variation in results can be directly attributed to Pr t .

Wall Boundary Conditions
Many stratified turbulent flows also occur in the vicinity of a boundary such as the SABL.Given the simulation cases, the effect of the wall boundary must be considered.For the k-ε model, standard smooth-wall functions are employed to calculate the values for k, ε, and u at the edge of the logarithmic layer (y + ≈ 30).This wall function methodology assumes the flux of k is zero through the wall surface resulting a zero-gradient wall boundary condition (dk/dz) w = 0, where the subscript w denotes a value at the wall surface.The shear stress over the wall-adjacent cell can be used to connect the wall shear stress (τ w ) and the mean streamwise velocity (u).
where u τ = C 1/4 μ k 1/2 is the shear velocity defined using the local equilibrium assumption, κ is the von Kármán constant, z + = zu τ /ν is the nondimensional height above the surface, the superscript + denotes a normalized quantity by the viscous scales u τ or ν, E = exp(κB) is an empirical constant related to the thickness of the viscous sublayer, and B ≈ 5.2.For the k-ε model, the von Kármán constant has an implied value, where κ ≈ 0.433 for the standard model constants (Table 1).Assuming that the shear production at the wall is given by P ≈ τ w (du/dz), the velocity derivative can be given by where the subscript P denotes the value at the first grid point.Finally, the dissipation rate of turbulent kinetic en-ergy at the first grid point (ε P ) is given by In general the equilibrium assumption used to develop these wall functions holds for continuously turbulent stratified flows.

Simulation Details
One-dimensional simulations of stably stratified turbulent plane Couette and channel flow are performed to assess the performance of Pr t propositions.Recently, García et al. [34] performed DNS studies for stably stratified turbulent Couette flow with a shear Reynolds number of Re τ ≈ 540.
The relevant flow and stratification specifying parameters are the shear (or friction) Reynolds number and Richardson number given by where u τ = (τ w /ρ 0 ) 1/2 is the shear velocity, τ w = μ(du/dz) is the wall shear stress, δ = h/2 is the half-channel depth, h is the channel depth, and Δρ = ρ bottom − ρ top is the density difference between the top and bottom surfaces.Figure 1(a) displays a schematic for this case.For plane Couette flow, the top wall is specified a velocity (U w ) and the bottom wall is stationary.The fluid viscosity in the nearwall region drives the flow velocity in the channel.The top and bottom surfaces are isothermal walls with Dirichlet boundary conditions for density.This study considers the continuously turbulent stratified cases of (i) Ri τ ≈ 83.5 and (ii) Ri τ ≈ 167.Additionally, García-Villalba & del Álamo [35] performed DNS studies for pressure-gradient driven wallbounded flow.This study considers channel flow simulations with Re τ = 550 and the continuously turbulent stably stratified states of (i) Ri τ = 60 and (ii) Ri τ = 120.Figure 1(b) displays a schematic for the channel flow scenario.The flow is driven by a pressure gradient in the x-direction calculated by −∂p/∂x = ρ 0 u 2 τ /δ.The top and bottom surfaces are isothermal walls with Dirichlet boundary conditions for density.In the channel core, the mean shear diminishes which leads to countergradient fluxes.
Table 2 presents the flow and stratification parameters for simulations.Simulations are preformed on the same computational grid.The top and bottom surfaces (z-direction) are defined as walls with no-slip conditions (u = 0) and Dirichlet conditions for mean density of ρ top and ρ bottom , respectively.Appropriate wall boundary conditions are specified based on the previous discussion in Section 2.3.The inlet and outlet surfaces (x-direction) are periodic boundary conditions.The front and back surfaces (y-direction) are treated as empty boundaries (i.e., no flow in the y-direction).The shear Reynolds number is monitored to match the prescribed values.
Numerical simulations were performed in the opensource computational fluid dynamic (CFD) software   OpenFOAM [27].OpenFOAM is a library of C++ utilities and solvers allowing for solutions to fluid flow problems using the finite-volume method.The standard k-ε turbulence model is implemented with the selected Pr t propositions to calculate turbulent diffusivity (κ t ) in the definition of buoyancy production (B) and the transport equation for mean density (Eq.11).Assuming steady-state homogeneous turbulence, the semi-implicit method for pressurelinked equations (SIMPLE) algorithm is selected.Spatial discretization is set to second-order accuracy.Solutions ran until appropriate convergence of all relevant flow and turbulent quantity residuals (R 10 −7 − 10 −8 ).Open-FOAM was selected for this study due to the open-source nature allowing for direct implementation and modification of Prandtl number formulations and turbulence models.OpenFOAM also provides thorough libraries which allow for the models tested on simple geometries to be scaled exponentially.

Results
This section presents the results of numerical simulations of stably stratified turbulent Couette and channel flow.The DNS data of Gracía et al. [34] (digitized) and García-Villalba & del Álamo [35], respectively, is used to evaluate model efficacy.

Couette Flow
Plane Couette flow simulations are performed for Re τ ≈ 540 with stratification levels of Ri τ ≈ 83.5 and 167.Mean velocity and density fields present an overall evaluation of Pr t model performance with all turbulence and numerical model parameters kept constant in simulations.
Couette flow gives a clearer picture of Pr t propositions because the mean shear does not dissipate in the channel core.Figure 2 presents nondimensionalized mean streamwise velocity (u + = u/u τ ) and mean density (ρ ) profiles compared with the digitized DNS data of García-Villalba et al. [34].The results of the Couette flow simulations reveal that all models compare very well with the DNS data.Marginal differences exist between the models at lower stratification (Ri τ ≈ 83.5).The mean velocity profile is slightly overpredicted and mean density profile is undermixed for all models except for KM92.This result can be attributed to the lower limit of Ri f ∞ = 0.2 in the model assumptions.For the higher stratification level (Ri τ ≈ 167), KM92 slightly overpredicts the velocity profile and undermixes the density profile.The remaining models all compare very well with the DNS data of García-Villalba et al. [34].These results further illustrate that a stability-dependent turbulent Prandtl number (Pr t ) is capable of turbulence parameterization and prediction of Couette flow.

Channel Flow
Pressure-driven channel flow simulations are performed with Re τ = 550 for stratification levels of Ri τ = 60 and 120. Figure 3 presents nondimensionalized mean streamwise velocity and mean density compared with the DNS data of García-Villalba & del Álamo [35].For Ri τ = 60, the constant Pr t predicts the mean velocity profile rather well, while VS10 slightly lags the DNS and KM92 varies significantly near the free surface.All of the models overmix the mean density field which is of no minor consequence.For higher stratification of Ri τ = 120, the mean velocity profile from KM92 falls well behind DNS profile.VS10 predicts the mean velocity with the most accuracy for both cases.As for the mean density, all simulations reveal overmixed profiles compared with DNS results.

Modified Turbulent Prandtl Number Formulation
Clearly, all Pr t formulations overpredict the mean density field for channel flow simulations.Karimpour & Venayagamoorthy [36] (hereafter, KV14) proposed a turbulent Prandtl number formulation for stably stratified wall bounded flow based on the observation that Pr t decays towards the neutral value near free stream (or core) for a channel flow.Their proposition is defined in Eq. 20.where Ri f = 0.25[1 − exp(−γRi g )] is a sole function of Ri g providing a simple empirical fit to the model of Mellor & Yamada [37], γ = 7.5 is an empirical constant, (1 − z/δ)Pr twd0 + Pr t0 ≈ 1.1 is the turbulent Prandtl number for a simple channel flow, and Pr t0 = 0.7 based on VS10. Figure 4 displays the results of channel flow simulations of KV14 compared with the other Pr t propositions.Clearly, KV14 significantly improves the prediction of the mean density field while accurately capturing the mean velocity field as well.
Results of channel flow simulations reveal that appropriate parameterization of Pr t can strongly influence the efficacy of model prediction of mean velocity and density.Mean velocity profiles that lag the DNS data indicate an overprediction of turbulent viscosity.For stably stratified flows, ν t is connected to the prediction of Ri f (Eq.14).Overmixing of the mean density profiles indicates the Pr t propositions underpredict the turbulent diffusivity (κ t ) due in part to the asymptotic behavior of gradient Richardson number (Ri g ) near the mid-channel depth (δ) with the minimization of shear.As stratification increases, this band of diminishing shear narrows in the channel core.

Conclusions
The research presented in this study investigate the performance of k-ε turbulence model and turbulent Prandtl number formulations for simulation of stably stratified boundary layer flows.For plane Couette flow, there is little differentiation between a constant Pr t = 0.85 and the formulations of KM92 or VS10 for mean velocity or density fields.In the case of channel flow, all three of these formulations failed to accurately predict the mean density field, E3S Web of Conferences 04003-p.6 significantly overpredicting mixing.This observation warranted a further investigation of the turbulent Prandtl number, including the proposition of KV14 for stably stratified wall-bounded flows.KV14 corrects of the inhomogeneities of flow close to a wall surface and significantly improves the prediction of mean density for channel flow simulations.
The research presented in this study indicates that the proper assignment of the turbulent Prandtl number in a RANS framework can significantly influence flow prediction.Regarding future research on the stably stratified atmospheric boundary layer, further investigations are needed to establish the behavior of the turbulent Prandtl number.Consideration must be given to variations in stability with the temporal evolution of the flow.Spatial variability must also be studied looking at the vertical profile of the SABL transitioning from the inner to outer layer dynamics.
Further implications of this study extend to modeling communities where stably stratified turbulence effects are significant in boundary layer flows such as wind energy.In the atmospheric boundary layer, diurnal variation leads to varying stratification conditions including stably stratified turbulent events during the nighttime hours.Improved wind speed predictions of stably stratified turbulent conditions could assist in the operation of wind turbines located in regions with substantial nighttime stratification effects (e.g., Great Plains of the United States).The framework can also be extended to additional RANS models and the subgrid-scale approximations in LES models in modeling stably stratified phenomena such as low-level jets.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 2015