An Explicit Algebraic Closure for Passive Scalar-Flux: Applications in Channel Flows at a Wide Range of Reynolds Numbers

: In this paper, we propose an algebraic model for turbulent scalar-ﬂux vector that stems from tensor representation theory. The resulting closure contains direct dependence on mean velocity gradients and quadratic products of the Reynolds stress tensor. Model coefﬁcients are determined from Direct Numerical Simulations (DNS) data of homogeneous shear ﬂows subjected to arbitrary mean scalar gradient orientations, while a correction function was applied at one model coefﬁcient based on a turbulent channel ﬂow case. Model performance is evaluated in Poiseuille and Couette ﬂows at several Reynolds numbers for Pr = 0.7, along with a case at a higher Prandtl number ( Pr = 7.0) that typically occurs in water–boundary interaction applications. Overall, the proposed model provides promising results for wide near-wall interaction applications. To put the performance of the proposed model into context, we compare with Younis algebraic model, which is known to provide reasonable predictions for several engineering ﬂows.


Introduction
Reynolds-Averaged Navier-Stokes (RANS) and Large Eddy Simulations (LES) approaches are today widely used in engineering applications for predicting the dynamics and scalar transfer of turbulent flows. Both approaches have existed for a long time and will remain useful for understanding several industrial configurations for several decades, playing different roles. In a LES, the larger scales in a flow are computed explicitly, but the effect of the smaller unresolved scales is modeled. The basic ingredient of LES modeling is a model for the sub-grid dynamics in terms of the resolved velocity field: a so-called eddy-viscosity model. An important example of an eddy-viscosity model for LES is the classical Smagorinsky model [1]. The Smagorinsky model gives good results for homogeneous isotropic turbulence, but removes too much kinetic energy from laminar and transitional flows. Minimum-dissipation eddy-viscosity models are a class of sub-filter models for large-eddy simulation that give the minimum eddy dissipation required to dissipate the energy of sub-filter scales. A class of minimum-dissipation that has many desirable properties cases, for example two-dimensional free shear flows and configurations involving curvature effects ( [16]). However, the linear nature of this specific closure is the main reason why it fails to capture the proper near-wall misalignment levels of scalar-flux vector, thus revealing the importance of incorporating nonlinear information regarding turbulence anisotropy. Therefore, in this study we propose an algebraic model that stems directly from the Younis formulation and involves products of the Reynolds stress tensor, a dependence missing from the multi-linear model. Special attention is given so that model complexity is kept minimal for ease implementation in existing industrial codes, while its performance ability is tested on several Couette and Poiseuille flows at a wide range of Reynolds numbers and a Prandtl number equal to 0.7 and 7.0. The former Pr number is encountered in aerodynamic applications, whereas the latter number typically occurs in marine applications and water-boundary interaction.
A description of the exact transport equations of the turbulent passive-scalar fluxes along with an extensive summary of the Younis approach are given in Section 2 with the motivation behind the use of a nonlinear term involving products of Reynolds stress given in Section 3. The proposed formulation is discussed in Section 4 and the performance of this closure is evaluated in Section 5, yielding encouraging results. Summary and conclusions are given in Section 6.

Younis' Formulation
The motivation behind the work reported by [15] arose from the need to provide a better alternative to the existing gradient-transport closures for the turbulent scalar-fluxes. As a starting point, they considered the exact transport equation of the following fluxes, where the instantaneous flow variables are decomposed into a mean part, denoted by overbar, and a fluctuating part, denoted by prime symbol. Hereafter, we are using index notation whereby repeated indexes imply summation. The coefficient of fluid viscosity and the coefficient of scalar diffusivity are denoted by ν and γ, respectively, ρ is the density of the fluid and u i , φ are the instantaneous fluid velocity and passive scalar fields, respectively. The fluctuating pressure is denoted with p . The first and second terms on the RHS of Equation (1) represent the generation of scalar-flux as a consequence of turbulence-mean deformation interactions. The third term refers to the rate at which the scalar-flux is destructed, while the fourth term refers to the fluctuating pressure-scalar correlations and is responsible for the redistribution of the flux among the different components. The last term is interpreted as the turbulent-transport term. [15] used Equation (1) as a guide to provide a rationally assumed relationship between the scalar-flux vector and various tensor quantities, where R ij = u i u j is the Reynolds stress tensor, is the energy dissipation rate, φ is half the scalar dissipation rate, and φ 2 is the scalar-variance. S ij , Ω ij , and Λ i denote the mean strain rate tensor, the mean vorticity tensor, and the mean scalar gradient vector, respectively, defined as With the aid of tensor representation theory, Younis et al. constructed an explicit model for u i φ in accordance with Equation (2). Under this approach, u i φ can be expressed as a series of basis vectors Υ i , where α n can depend on all the tensor variables appearing in Equation (2). The basis vectors are formed from the products of the symmetric (S ij , R ij ), the skew-symmetric (Ω ij ) tensor and the vector (Λ i ), leading to the following algebraic expression.
To bring the above expression into a more compact form, [15] assumed sufficiently small anisotropies and turbulent time scales. Further assuming that the effects of S ij and Ω ij are in balance, results to the following multi-linear expression, where κ = R kk /2 is the turbulent kinetic energy, while the α i coefficients were replaced based on dimensional arguments. The first term on the RHS of Equation (6) represents the GDH, while the second term coincides with the production term that appears in Equation (1) and is associated with the mean scalar deformation, also known as the Generalized Gradient Diffusion Hypothesis (GGDH) model proposed in [9]. The remaining two terms involve products between the gradients of mean scalar and mean velocity a dependence proposed by the analyses of the works in [17,18]. The values of C i coefficients were determined based on the LES results of [19] for homogeneous flows subjected to a uniform shear with uniform scalar gradients Regarding inhomogeneous flows, detailed analysis of the model performance in relation to data from a large number of studies on wall-bounded flows showed a serious error in the normal flux component [16]. For example, the normal component of (6) takes the following form for flat-wall flows Using DNS data [20] for a fully-developed turbulent channel flow, it can be seen that near the wall (0 < y + < 24). Using the above expression into Equation (8) incorrectly predicts that u 2 φ and Λ 2 have the same sign. As a result, the value of C 1 coefficient was modified in the near-wall region by applying the following damping function, where α = −0.02 and β = 1.9 are the values proposed in [21], Pr = ν/γ is the Prandtl number, Re t = κ 2 ν is the turbulence Reynolds number, and A is the stress-flatness parameter, defined as where A 2 and A 3 are the second and third invariants of the normalized Reynolds-stress anisotropic tensor A ij . More details regarding the mathematical formulation can be found in [16]. As already mentioned in the introduction, the linear form of Younis closure is an important reason why this closure exhibits certain limitations. For example, it fails to capture the near-wall relative magnitudes between the scalar-flux components. In the following section we discuss the importance of incorporating a term that contains nonlinear information regarding the turbulent anisotropy, vital to predict the proper near-wall behavior.

Quadratic Dependence on Reynolds Stress
It is well known that models involving only the second term of Equation (6) cannot predict the streamwise scalar-flux component reasonably well [22]. In an attempt to extend their applicability, [23] performed a series of LES simulations for fully developed turbulent channel flows under different boundary conditions and for a wide range of Prandtl numbers. They relied on the findings of [24], who pointed out that scalar fluctuations are correlated more strongly with streamwise than transverse velocity fluctuations in the near-wall region, to propose an algebraic relation between the turbulent scalar-flux vector and quadratic products of Reynolds stress tensor: where τ is a turbulent time scale and C φ is a model coefficient. In wall-bounded channel flows where statistical quantities are functions only of the wall-distance, the above expression approaches the following near-wall limit, where subscripts 1, 2, and 3 are, respectively, the streamwise, wall-normal, and spanwise directions. [23] also observed that the correlation between φ and u 2 relatively increases with the decrease of Pr, reaching almost the same level as that between φ and u 1 for Pr = 0.025. As a result, they recommended that an effective algebraic scalar-flux model should depend both on linear and quadratic forms of Reynolds stress, written as where C φ,1 and C φ,2 are model coefficients. These coefficients should be determined so that the nonlinear term in Equation (14) becomes dominant in regions of high deformation rates, while the linear term becomes significant in the presence of weak strain rates, as well as in low Prandtl numbers. Hence, the present study aims to propose an algebraic model based on Younis general form (5) that involves the nonlinear term, and validate its performance in different channel flows.

Proposed Formulation
In this section, we introduce the proposed model, which is essentially a combination between the multi-linear model of Younis (6) and the functional form proposed by Abe and Suga (14), thus incorporating nonlinear information regarding turbulence anisotropy. Special attention is given to keep the model as simple as possible so that it can be easily implemented in existing industrial codes. This is achieved by investigating the contribution of each term appearing in the above equations in an attempt to propose a minimal combination of these terms that is able to provide reasonable predictions. We chose to neglect the C 1 -related term by assuming that information regarding this term is partly included in the isotropic part of the second term of Equation (6). We have also excluded the term associated with C 4 , even though its inclusion could provide additional flexibility, for the following reasons. First, this term cannot induce the proper near-wall anisotropy levels of the scalar-flux vector, as given in Equation (13). Second, this term involves products between the mean velocity gradients and the Reynolds stress, thus being more complex than the C 3 -related term, which also contains gradients of the mean velocity field. Third, its absence facilitated our effort to apply minimal corrections to the model coefficients, as will be shown in this Section. Consequently, the proposed closure takes the following compact form, where we keep the same indexing for the model coefficients as in Equation (6). The above equation essentially differs from Equation (14) in involving the C 3 -related term, which accounts for the products between mean scalar and mean velocity gradients.

Redistributive Term
In order to bring closure to the transport equations of the turbulent scalar-fluxes and stresses, it is required to model correlations that involve the fluctuating pressure p . As the action of the pressure interactions is known to redistribute the turbulent kinetic energy among the diagonal components of the Reynolds stresses, it is a common practice to re-express the pressure term that exist in the stress equations as the sum of a traceless distributive term, called the pressure-strain correlation, and the divergence of the correlation between the fluctuating velocity and the pressure field (added to the turbulent diffusion). As pressure is also present in the scalar equations, an analogous treatment has been employed that splits the pressure-related term in two parts: a correlation term p ∂φ ∂x i that redistributes the scalar-flux among its different components, and its divergence that is absorbed in the turbulent-transport term which explains the terms appearing at the RHS of (1). Given the level of approximation and empiricism that, inevitably, is needed in order to derive an algebraic model form, it is not possible to associate one particular term in the resulting algebraic expression (15) with the pressure-correlation term. Several models have been developed to model this term [25][26][27], suggesting that no consensus has been achieved regarding the modeling of this term. In an attempt to associate the terms appearing in the proposed closure (15) with the pressure-scalar correlation, we consider the elegant analysis conducted in [17], who represented the fluctuating parts of scalar and pressure as Fouries series in order to derive the following expression that relates pressure-correlation term with turbulent variables that appear in the transport equation where C θ,i are model coefficients. Comparing the proposed closure with the above equation, we observe that only the C 2 term is explicitly present, while the other two terms are implicitly associated through the remaining terms.

Model Calibration
One way to determine model parameters is to use benchmark inhomogeneous cases, such as channel flows, to assign a set of functional forms to the model coefficients through a fitting procedure. However, this process might lead to a complex set of functions that do not necessarily perform well for cases that deviate from the ones used to calibrate the model. Another method is to determine the model parameters through a two-step process: As a first step, initial values are assigned to the coefficients by fitting the proposed closure to existing homogeneous data, which are subsequently modified through the use of correction functions so that they include inhomogeneous effects. It is common to use these expressions so that they correct the near-wall behavior of the closure, thus called near-wall functions, while recovering the homogeneous values far away from the wall. Preliminary efforts to adopt the above approach in this study did not lead to satisfactory improvements regarding the streamwise flux. Eventually, we have proposed a correction function only for the C 3 term, denoted as f C3 , that exhibits the proper analytical behavior very close to the wall boundary (Appendix A.1), provides fair predictions for the near-wall peak value of the streamwise component for a wide range of Reynolds numbers (see Section 5), and does not necessarily recover the homogeneous value of the coefficient far away from the wall region. As in [15], homogeneous coefficients are determined through a linear regression fitting using data derived from the LES non-buoyant results of the work in [19], who considered homogeneous shear flows under different orientations of the mean scalar gradient. Next, we have considered a fully-developed channel flow at friction Reynolds number Re τ = u τ δ ν = 395, where u τ is the friction velocity and δ is the half-channel height, to modify the coefficient associated with the mean velocity gradients (C 3 ). The flow is driven by a uniform mean pressure gradient along the streamwise direction, with uniform scalar flux applied at both walls. The Prandtl number is set equal to 0.71, whereas the presented data is expressed in wall units. Figure 1 shows model predictions for scalar fluxes, obtained by importing DNS results of the works in [28,29] into Equation (15). For this case, only a streamwise component of the mean velocity field exists that varies along the normal component x 2 . Figure 1a shows model predictions for the streamwise component either with (corrected) or without (uncorrected) the inclusion of a correction function for the homogeneous C 3 coefficient, defined in Equation (19). We observe that the proposed model based on Equation (18) overestimates the near-wall peak value of the streamwise component, while very good predictions are provided for the normal component ( Figure 1b). As a result, we focus on improving only the streamwise component, without affecting the normal component. A simple way to do that is by modifying C 3 , as the corresponding term does not appear in the algebraic expression for the normal flux. We have concluded to the following expression, where S = S ij S ji /2 is the magnitude of the mean strain rate and τ = κ/ is the turbulent time scale. Details regarding the construction of the correction function are given in Appendix A.2. For y + < 15, the uncorrected C 3 -related term plays a dominant role (Figure 1c), with its maximum occurring at about y + = 10, a location at which mean strain rate is also maximized (not shown here). The C 5 -related term contributes the most for y + > 15, while the impact of C 2 -related term becomes non-trivial in regions far away from the wall that are characterized by weak deformations. Applying the correction function on C 3 coefficient yields a reduction of the associated term, which remains important in the high-shear region, while resulting in the quadratic term being dominant in the entire region outside the viscous sublayer (y + > 7). By contrast, the C 3 -related term does not contribute to the normal component, which is prevailed by the nonlinear term for y + > 15. Consequently, the near-wall behavior of the model is driven by two mechanisms: one that involves products between mean velocity and scalar gradients, and a mechanism arising from the turbulence-turbulence interactions. Combining these two terms yields an alternative form for the proposed model: The above equation suggests that the non-linear interactions provide a gradient, in addition to the actual mean gradient, thus yielding an effective gradient G e ij : This idea, called the "effective gradients" hypothesis, has been originally proposed by [30] to construct a Reynolds-stress transport closure, and has been recently extended by the authors of [31][32][33] for passive scalar transport. Thus, Equation (20) can be thought as an extension of Abe and Suga proposal (12), as it contains a linear term that becomes important in regions of weak deformation rate, while replacing the quadratic term by a term involving the effective mean gradient that dominates the region of high and moderate deformations rates.

Model Assessment
In this section we investigate the estimation ability of the proposed closure for channel flows under different boundary conditions and Reynolds numbers. For the examined cases considered, the Prandtl number equals to 0.71 and 7.0 and the flow variables are expressed in wall units. The criterion to choose those two Pr values for applying the developed algebraic model is to provide results that prove the efficient applicability of the proposed model in fluid flows related to aerodynamic and hydrodynamic conditions. In order to reduce uncertainties associated with the numerical solution of model equations, we import results from DNS into the proposed algebraic expressions for the scalar-fluxes. To facilitate the discussion during the validation procedure, we have performed additional computations to account for the performance of Younis' model, as described in Section 2. For all cases considered, only a streamwise component of the mean flow exists that varies along the wall-normal direction x 2 . Under these conditions, the mean velocity gradient G ij and the expressions for the flux components according to Equation (15) become and respectively. The flow geometry and coordinate system are shown in Figure 2.

Poiseuille Flows
We evaluate model performance at fully-developed turbulent channel flows for a wide range of friction Reynolds numbers, particularly Re τ = 150, 395, 640, 995, 2013 and 4088. For all cases, non-slip boundary condition is adopted in the wall-normal direction (i.e., at the top and bottom walls) while the scalar boundary condition is uniform scalar-flux on both walls.

Low Reynolds Number Cases
We initially consider three different low Reynolds number cases, namely, Re τ = 150, 395, and 640. The quality of the predictions is compared with the corresponding DNS results [29,34,35]. Figure 3 compares the turbulent scalar-fluxes as obtained from the present model and the Younis model for the three cases. Solid line denotes the proposed model and dash-dotted line denotes Younis' model. Symbols denote the DNS results. For all cases, both models provide reasonable agreement with DNS regarding the near-wall peak value of the streamwise component. In addition, both closures underestimate the value of this component further away from the wall, while the predictions of the proposed closure are closer to the DNS results. Contrary to algebraic closures, DNS predict non-zero values for the streamwise component at the channel centerline. This happens because the turbulent scalar field is being transported by the mean flow or the turbulence. These transport processes constitute a non-local mechanism that cannot be captured by a rational algebraic closure [15]. Except from the lowest Re τ case, the proposed model achieves a considerably better agreement than Younis' model regarding the normal flux at the entire region. In particular, the current model accurately captures the near-wall behavior of the fluxes, while it exhibits a mild overestimation of the magnitude beyond y + ≈ 80 for Re τ = 640, still being able to produce satisfactory results.

High Reynolds Number Cases
Next (figure 4), we investigate model performance at higher Reynolds numbers than before, particularly Re τ = 995, 2013 and 4088, against available DNS results [36,37].
Regarding the streamwise flux component, the current model captures well the near-wall peak magnitude for all cases, being able to properly adjust to the presence of high-shear rates. Furthermore, it achieves very good agreement with the DNS predictions for the normal component at the near-wall region (y + < 40) for all cases, suggesting that the use of the terms involving Reynolds stress is enough to properly capture the highly anisotropic correlations between u 2 and φ that exist at the near-wall region. As in the Re τ = 640 case, we observe its tendency to overestimate, although much milder than Younis' model, the magnitude of the normal flux beyond the slope change.

Inclination Angle of the Scalar-Flux Vector Near the Wall Boundary
Here, we further investigate the near-wall performance of the proposed closure. A useful parameter for that purpose is the scalar-flux ratio R u φ , defined in Equation (13). Figure 5 shows a comparison of the scalar-flux ratio predictions of both models with the corresponding DNS data at three different Reynolds number cases, namely, Re τ = 395, 995 and 4088. For comparison purposes, we have also included the DNS predictions for R 12 /R 11 as this ratio is expected to vary as scalar-flux ratio in the near-wall region, as already mentioned in Section 3. We observe that the proposed model achieves better agreement with the DNS results for all cases. Especially for Re τ = 395, it captures accurately the near-wall limit up to y + ≈ 30, in accordance with the previous discussion. This is attributed to the dominant role that C 5 -related term plays in the buffer layer for both flux-components, as already shown in Figure 1, which exhibits the proper near-wall physical behavior (see discussion in Section 3). In contrast, Younis' model fails to capture this limiting behavior for y + > 10, partly because no term appearing in its model Equation (6)

Couette Flows
We now ascertain the performance of the models in a fully-developed Couette flow with the top wall moving at constant speed u w relative to the bottom wall. As a result, the mean flow is driven by the shear stress due to the relative movement between top and bottom walls, while the scalar difference between the top and bottom walls is kept constant. We consider two cases at relatively low Reynolds numbers, particularly Re w = 2 u w δ ν = 8600 and 12,800, for which detailed DNS data are provided in [28]. Figures 6 and 7 show a comparison between the two models for the scalar-flux components with the corresponding DNS results. Compared to Younis closure, the proposed model tends to overestimate the near-wall peak value of the streamwise component, while it provides better estimations beyond y + > 20. Regarding the normal component, predictions agree quite well with the DNS data for both cases, especially for the higher Re w -case, while Younis' model underestimates the wall-normal component outside the viscous sublayer (y + > 10).  Figure 6 but for Re w = 12,800.

Higher Pr Case
As mentioned in Section 3, Kim and Moin [24] studied the near-wall behavior of the scalar-flux vector in turbulent channel flows for a wide range of Prandtl numbers. As a result, they observed mild changes in scalar ratio R u φ with increasing Prandtl number in the range of Pr ≥ 0.71. The above observation, along with the fact that the proposed model does not explicitly depend on this parameter, has motivated us to evaluate model performance at higher Pr than 0.7. In particular, we consider a Poiseulle flow at Re τ = 395 and Pr = 7, for which detailed DNS data are provided by Kawamura [28]. Boundary conditions are identical to the ones described in Section 5.1. Comparison is made again with Younis' model, which directly depends on Pr through damping function f C1 as indicated in Equation (10). Figure 8 shows a comparison between both models with the corresponding DNS results for the scalar-flux components. Compared to Younis' model, the proposed model tends to underpredict the near-wall peak value of the streamwise component, whereas it provides better predictions for y + > 15. As in all previous cases, both models exhibit a sharp reduction of the profile while moving away from the walls. Regarding the normal component, the present model achieves a much better agreement with the DNS results outside the viscous sublayer (y + > 10).  Figure 9 shows the corresponding predictions for the scalar ratio R u φ . The proposed model agrees reasonably well with the DNS results for y + < 15. Compared to Younis' model, the proposed model provides much better predictions for this ratio in the entire domain.

Summary and Conclusions
In this study, we have proposed an explicit algebraic closure for the turbulent scalar-flux vector based on the general formulation of Younis, which is essentially an extended version of the model proposed by Abe and Suga (14) through the inclusion of an extra term involving products between the gradients of the mean velocity and scalar field. This closure consists of three terms, making it simpler than Younis' model and an elegant choice for use in general-purpose computational codes. It is also motivated by the "effective gradients" hypothesis, which postulates that turbulence-turbulence interactions provide a gradient that acts supplementary to the mean shear, thus giving an alternative physical interpretation of the proposal compared to other closures. To minimize model bias to inhomogeneous applications, the values of the model coefficients are determined based on existing LES predictions of homogeneous shear flows in the presence of arbitrary mean scalar gradients, while a simple correction function was applied to the C 3 -term to properly capture the near-wall behavior. In order to test the quality of the proposed model, we have considered different types of channel flows, particularly Poiseuille and Couette flows, for a broad range of Reynolds numbers and Pr = 0.7. The proposed model performs satisfactorily in Couette flows at two different Reynolds numbers, especially for the normal component, thus showing its sensitivity to the near-wall inclination of the scalar-flux field. Regarding Poiseuille flows, fair predictions are obtained for both flux components, with the model being able to well capture the near-wall peak magnitude of the streamwise component regardless the Reynolds number choice. Furthermore, the proposed model accurately captures the near-wall behavior for the scalar ratio R u φ for distances up to y + ≈ 30 (Re τ = 395), significantly further than Younis' model (y + ≈ 10). Model performance is further evaluated in a Poiseuille flow at a higher Prandtl number (Pr = 7) that typically occurs in marine applications and water-boundary interaction, again providing reasonable predictions. The above cases served as a guidance regarding the role of the different terms appearing in the algebraic expressions, showing the ability of the proposal to capture the levels of misalignment between the different flux components. As part of future work, we intend to propose an extension of the proposed RANS-based model so that it properly accounts for Prandtl effects under different scenarios. Special attention will be given to Pr numbers that typically occur in aerated flows. To achieve that, we intend to consider several cases that involve Prandtl numbers in the range of 1 to 10, and perform computations to examine the sensitivity of the closure to the predicted stress fields. Furthermore, relevant experimental studies will be performed in the future to provide data for a comparative study between different turbulence models for the case of water flow around a monopile.  For a fully-developed turbulent channel flow, flow statistics vary only along the wall-normal (x 2 ) direction. For fixed walls, the analytical limiting behavior of turbulent quantities is whereas the stress flatness parameter A becomes thus varying as O(y 2 ) near the wall boundary. Regarding the scalar field, assuming constant wall values yields ∂φ ∂x 2 ∝ O(y 0 ) .
The above expressions suggest the following near-wall behavior for the scalar-flux components, Substituting Equations (A1)-(A4) into modeling expressions (23) leads to the following asymptotic expressions for the model coefficients, The term involving C 3 coefficient is expected to be important close to the wall boundary, due to the presence of high deformation rates. To meet this requirement, we choose this term to vary as the streamwise scalar-flux component, which leads to the following condition,