Pressure scrambling effects and the quantification of turbulent scalar flux model uncertainties
PPressure scrambling effects and the quantification ofturbulent scalar flux model uncertainties
Zengrong Hao, Catherine Gorl´e
Wind Engineering Laboratory, Department of Civil and Environmental Engineering,Stanford University, Stanford, CA 94305, USA
Abstract
Closure models for the turbulent scalar flux are an important source of un-certainty in Reynolds-averaged-Navier-Stokes (RANS) simulations of scalartransport. This paper presents an approach to quantify this uncertainty insimulations of complex engineering flows. The approach addresses the uncer-tainty in modeling the pressure scrambling (PS) effect, which is the primarymechanism balancing the productions in scalar flux dynamics. Inspired bythe two classical phenomenological theories of return-to-isotropy (RI) andisotropization-of-production (IP), we assume that the most likely directionsof the PS term are around a fan-shaped region bounded by the RI and IPdirections. Subsequently, we propose a strategy that requires two additionalsimulations, defining perturbations of the PS directions towards the RI andIP limits. The approach is applied to simulations of forced heat convection ina complex pin-fin array configuration, and shows favorable monotonic proper-ties and bounding behaviors for various quantities of interest (QoIs) relevantto heat transfer. To conclude, the results are analyzed from the perspectiveof transverse scalar transport in a shear flow; the analysis indicates that theproposed approach is likely to exhibit monotonic behaviors in a wide rangeof scalar transport problems.
Keywords:
Turbulent scalar flux model, Model form uncertaintyquantification, Pressure scrambling effects, Return-to-isotropy,Isotropization-of-production
Preprint submitted to Physical Review Fluids April 16, 2020 a r X i v : . [ phy s i c s . f l u - dyn ] A p r . Introduction Reynolds-averaged-Navier-Stokes (RANS) simulations of turbulent flowswith passive scalar transport require closure models for the Reynolds stressesand the turbulent scalar fluxes. The epistemic uncertainty in both modelforms can affect the accuracy of the results, and tools to assess the corre-sponding uncertainty in predicted quantities of interest (QoIs) would providevaluable information when using the simulations for engineering design andanalysis. Epistemic uncertainty quantification (UQ) of Reynolds stress mod-els has been explored in several recent studies [1, 2, 3, 4, 5, 6], but quantifyingthe uncertainty introduced by scalar flux models has received comparativelyless attention.The most common approach to scalar flux modeling for engineering appli-cations is based on the standard gradient-diffusion hypothesis (SGDH). Thishypothesis relates the scalar flux to the product of a scalar diffusion coeffi-cient and the mean scalar gradient. The scalar diffusion coefficient is definedas the ratio of the turbulent viscosity and a turbulent Prandtl number, Pr t (or Schmidt number, Sc t ). When investigating the performance of the SGDH,the primary focus has been on the influence of the selected values for Pr t orSc t . Their considerable effect on the predicted scalar field has been demon-strated in a variety of flow problems, including pollutant dispersion in urbanareas [7, 8], reacting flows in propulsion systems [9, 10], particle burning inpacked beds [11], and film cooling in turbomachines [12, 13]. Calibration ofPr t or Sc t has also indicated that it is far from a universal constant; it ishighly dependent on the flow configurations and regimes, with a reportedrange as large as 0 . ∼ . t primarily demon-strate the effect of increasing or decreasing the magnitude of the scalar fluxvector. They do not investigate the effect of the assumption that the direc-tion of the vector remains identical to the local mean scalar gradient, whichhas been demonstrated to be incorrect in a priori analyses of scalar trans-port problems. A formal method for UQ of scalar flux models should alsoaddress this assumption. A framework to achieve this was first proposed byGorl´e and Iaccarino [10]. The method used a generalized gradient diffusionhypothesis (GGDH [16]), which defines a tensorial diffusion coefficient thatis a function of the Reynolds stresses, and it propagates perturbations to theReynolds stresses through the GGDH model to modify the predicted scalarflux. When defining perturbations based on a high-fidelity data set it was2hown that this framework could provide a prediction with a confidence inter-val that encompasses high-fidelity data; this could not be achieved when onlymodifying Pr t . In [17], this approach was extended to consider the propaga-tion of physics-based perturbations to the Reynolds stresses through differentscalar flux models, including a second-moment model that solves three addi-tional transport equations for the scalar fluxes. The results showed that theperturbations provide an interval prediction for the average Nusselt numberin a pin-fin heat exchanger that encompasses the value obtained from high-fidelity simulations; however, the intervals predicted for the local distributionof the Nusselt number did not fully encompass these data. It was suggestedthat this shortcoming is related to uncertainty in the modeled scalar flux.In this paper, we propose a UQ method that addresses the inherent inade-quacies of models for the dynamics (i.e. the transport equation) of scalar flux.We particularly focus on the pressure scrambling (PS) term, which is the pri-mary mechanism acting to balance the productions. Since the mechanismof PS is not fully understood to date, this term is also a major uncertaintysource of second-moment models in practice. Therefore, the basic idea of thiswork is to ‘bound’ the (vectorial) PS term through a reasonable approachbased on our limited understanding of its effects.
2. Scalar flux model UQ method
For incompressible turbulent flow, with constant density ρ , molecular vis-cosity ν , and molecular scalar diffusion coefficient κ , the transport equationfor a mean passive scalar Θ is given by:¯ D ¯ Dt Θ = ∂∂x i (cid:18) κ ∂ Θ ∂x i − u (cid:48) i θ (cid:48) (cid:19) , (1)where U i and Θ are the mean velocity and scalar fields, u (cid:48) i and θ (cid:48) the corre-sponding turbulent fluctuations, and ¯ D/ ¯ Dt ≡ ∂/∂t + U i ∂/∂ i is the materialderivative. The term u (cid:48) i θ (cid:48) in Eq. (1) is the turbulent scalar flux to be closed.The exact transport equation for u (cid:48) i θ (cid:48) reads:¯ D ¯ Dt u (cid:48) i θ (cid:48) = − u (cid:48) i u (cid:48) j ∂ Θ ∂x j − u (cid:48) j θ (cid:48) ∂U i ∂x j + p (cid:48) ρ ∂θ (cid:48) ∂x i − ( ν + κ ) ∂u (cid:48) i ∂x j ∂θ (cid:48) ∂x j + ∂∂x j (cid:32) νθ (cid:48) ∂u (cid:48) i ∂x j + κu (cid:48) i ∂θ (cid:48) ∂x j − u (cid:48) i u (cid:48) j θ (cid:48) − p (cid:48) ρ θ (cid:48) δ ij (cid:33) . (2)3n the right-hand-side of Eq. (2), the first two terms, denoted as G I i and G II i hereafter, represent the direct interactions between the mean gradientfields and the large scale turbulent eddies. These interactions are immedi-ately responsible for the production of u (cid:48) i θ (cid:48) ; they are in closed form given theReynolds stress field u (cid:48) i u (cid:48) j . The third term, denoted hereafter as Π i , repre-sents the pressure scrambling (PS) effects. This is the primary mechanismbalancing the two production terms, but it is not well understood to date,primarily due to the highly non-local nature of pressure. The fourth term,denoted hereafter as ε θi , represents the molecular dissipation; it is usuallynegligible when the Reynolds and P´eclet numbers are both high enough, asis the case for the flow considered in this paper. The last four terms, denotedtogether as D θi , are all in divergence form and represent the diffusive trans-port. Despite their complex form, it seems that employing a certain typeof GDH model (to be introduced in § Π i ; this motivates our focus on quantifying theuncertainty in this term. Existent understanding of the PS effects is informed by the wave natureof the pressure fluctuation p (cid:48) , which tends to disorganize the turbulent eddystructures and thus decorrelate two variables, in this case u (cid:48) i and θ (cid:48) . Further-more, the governing Poisson equation for p (cid:48) − ρ ∂ p (cid:48) ∂x i ∂x i = ∂ ∂x i ∂x j (cid:0) u (cid:48) i u (cid:48) j − u (cid:48) i u (cid:48) j (cid:1) + 2 ∂U i ∂x j ∂u (cid:48) j ∂x i (3)indicates that three separate mechanisms contribute to p (cid:48) : the slow pressure p ( s ) induced by the interaction between turbulent eddies; the rapid pres-sure p ( r ) induced by the interaction between turbulent eddies and meanflow distortion; and the harmonic pressure p ( h ) induced by the boundaryconditions, e.g. wall-blocking and free-surface reflection. Correspondingly,the wave effects of p (cid:48) = p ( s ) + p ( r ) + p ( h ) can also be discussed in terms ofthese separate mechanisms, and the PS term in Eq. (2) can be written as Π i = Π ( s ) i + Π ( r ) i + Π ( h ) i . In this analysis, we focus on the uncertainty inthe PS term resulting from the turbulent nature of the flow, i.e. on Π ( s ) i and Π ( r ) i . 4irst, a theoretical basis for modeling the slow part Π ( s ) i , is provided bythe return-to-isotropy (RI) hypothesis (Rotta [19], Monin [20]), which arguesthat the vector Π ( s ) i tends to diminish the scalar flux. This implies that itsdirection should be opposite the scalar flux direction: Π ( s ) i // − u (cid:48) i θ (cid:48) . (4)Second, the theoretical basis for modeling the rapid part Π ( r ) i , is providedby the isotropization-of-production (IP) hypothesis (Naot [21], Owen [22]),which argues that the vector Π ( r ) i tends to counteract the scalar flux produc-tion by mean flow distortion, i.e. G II i . This implies that its direction shouldbe opposite the direction of G II i : Π ( r ) i // − G II i . (5)Eqs. (4) and (5) have been used to establish PS models through calibra-tion in canonical flows. For example, a basic model has been proposed basedon calibration in homogeneous shear flow with a transverse scalar gradient[23]. The drawback of this approach is that the accuracy of these models cannot be guaranteed when extrapolating to more complex engineering flows:the lack of information on the magnitudes of the slow and rapid pressureterms in Eqs. (4) and (5) makes the model highly dependent on the cal-ibration process. However, from the perspective of model UQ, the limitedinformation provided by the RI and IP hypotheses supports identifying plau-sible directions of the combined PS term Π ( s ) i + Π ( r ) i : directions close to theRI-IP plane are more likely to occur than directions considerably deviatingfrom the plane. This provides the basis for the proposed UQ approach. As a starting point, we adopt the basic second-moment model [23, 24] asthe baseline (BSL) model for the UQ framework. In this model, the pressurescrambling terms are represented as: Π ( s ) i + Π ( r ) i = − c θ εk u (cid:48) i θ (cid:48) − c θ G II i , (6)with c θ = 2 . c θ = 0 .
4. The dissipation ε θi is assumed to be negligible,and the last four terms in Eq. (2) are modeled as: D θi = ∂∂x j (cid:32) c θd kε u (cid:48) j u (cid:48) k ∂u (cid:48) i θ (cid:48) ∂x k (cid:33) (7)5ith c θd = 0 . k ≡ u (cid:48) i u (cid:48) i / ε ≡ ν ( ∂ j u (cid:48) i ) ( ∂ j u (cid:48) i )the energy dissipation rate.Eq. 6 implies that the baseline PS term lies on the RI-IP plane. As dis-cussed in Section 2.2, its exact direction on this plane is an important sourceof model uncertainty, since it is determined by the relative magnitudes of theslow and rapid terms. Hence, we propose a simple perturbation strategy toquantify this uncertainty: for a local PS vector predicted by the BSL model,we perturb its direction within the RI-IP plane as shown in Fig. 1. The per-turbation requires the definition of the uncertain parameter r RI or r IP in theinterval [0,1], where 0 corresponds to no perturbation, while 1 correspondsto a perturbation to the RI or IP limit. In this initial analysis, we prescribea single value of r RI or r IP and apply it uniformly in the entire computationaldomain. Figure 1: The strategy of perturbing the PS direction in the UQ approach
3. Results
The test case considers heat transfer in a pin-fin array consisting of eightrows of staggered cylinders (pins) between two parallel flat plates (fins). Itwas studied experimentally by [25, 26, 27], and we have reported a high-fidelity LES [28] for the operating conditions considered in this study. Fig. 2shows the computational domain, which spans half of the lateral ( y ) pinspacing and half of the vertical ( z ) fin spacing due to symmetry in bothdirections. All pin and fin surfaces have a fixed and uniform temperature Θ w .At the inlet, fluid with a higher temperature Θ in enters the channel, inducing6 igure 2: Computational domain and mesh for the pin-fin array configuration forced heat convection between the fluid and the pin and fin surfaces. Thebulk Reynolds number is given by Re D = V m D/ν = 10 , where D is thediameter of the cylindrical pins and V m is the average velocity on the throatcross-section between two laterally adjacent pins. The Prandtl number Pr is0.71. Only steady solutions (i.e. ∂/∂t = 0) of Eqs. (1) and (2) are considered.Details on the numerical set-up can be found in [29].The main QoIs are the heat transfer rates on the pin and fin surfaces, char-acterized by the Nusselt number Nu ≡ ( ∂ n Θ ) w D/ ( Θ b − Θ w ), where ( ∂ n Θ ) w is the wall-normal temperature gradient and Θ b is the bulk temperature ofthe surrounding fluid. The PS term ( p (cid:48) /ρ ) ∂ i θ (cid:48) cannot be directly obtained from the LES sincethe contribution from unresolved fluctuations cannot be neglected. Hence,we estimate the term through the budget of Eq. (2), with D θi and ε θi modeledas in the BSL model: p (cid:48) ρ ∂θ (cid:48) ∂x i ≈ U j ∂u (cid:48) i θ (cid:48) ∂x j + u (cid:48) i u (cid:48) j ∂ Θ ∂x j + u (cid:48) j θ (cid:48) ∂U i ∂x j − ∂∂x j (cid:32) c θd kε u (cid:48) j u (cid:48) k ∂u (cid:48) i θ (cid:48) ∂x k (cid:33) (8)The fields U i , Θ , k , u (cid:48) i u (cid:48) j , and u (cid:48) i θ (cid:48) are all directly obtained from the LESdata, while ε is partially modeled as in [28].7o analyze the direction of the resulting PS vectors, we define two metrics m and m as shown in Fig. 3; these metrics quantify the discrepancy betweenthe estimated PS direction and the RI and IP directions. Figure 3: Definition of the metrics m and m , indicating the relation between the es-timated PS direction and the RI ( m = 0, m = − /
2) and IP ( m = 0, m = +1 / Fig. 4 shows the joint distribution of m and m , weighted by the PSterm magnitude, over the entire domain. The most likely PS directions con-centrate around the RI-IP plane, and 57% of the values are inside the ellipse (cid:112) m + m = 1 /
2. This result confirms the validity of the hypothesis thatdirections close to the RI-IP plane are more likely to occur than directionsconsiderably deviating from the plane.
Figure 4: Joint distribution of ( m , m ) to visualize the estimated PS directions, weightedby their magnitude, relative to the RI and IP directions over the entire test section of thepin-fin array. .3. Practical performance of the UQ method In this section we present results obtained from a series of simulationswith different r RI or r IP values. To evaluate the performance of the approachin absence of Reynolds stress model errors, the simulations solve Eqs. (1)and (2) using the BSL model with input for U i , k , ε and u (cid:48) i u (cid:48) j obtained fromthe LES data.First, Fig. 5 shows the globally averaged Nusselt numbers on the fins, Nu f for the full range of the uncertain parameter r IP or r RI . With r IP changingfrom 0 to 1, Nu f monotonically increases by up to 38%; in contrast, with r RI increasing, Nu f monotonically decreases, to a far lesser extent, by upto 3%. A similar trend is reflected by the temperature contours in Fig. 6.From the bottom to the top, the fluid cooling is gradually enhanced, asindicated by decreasing outlet temperatures. Given the monotonicity, and Figure 5: Globally averaged Nusselt numbers on the fins changed with the uncertainparameter the observation that the interval obtained from simulations with r RI , IP = 0 . f averaged over individual1 . D × . D small fin patches on the left, and Nu p averaged over the surfacesof individual pins on the right. Compared to the LES results, the BSL modelunderestimates the local values by 14 ∼
21% for Nu f , and by 11 ∼
14% forNu p . The perturbation with r IP = 0 . f , reaching alevel 8 ∼
17% higher than the LES results. Conversely, the perturbation with9 igure 6: Temperature counters on the middle plane of the pin-fin array r RI = 0 . f prediction to a level 16 ∼
23% lower thanLES. The resulting interval encompasses the LES predictions for Nu f in theentire test section. The results for Nu p exhibit a similar trend, altough the r IP = 0 . ∼ Figure 7: Locally averaged Nusselt numbers on the fins (left) and the pins (right)
Fig. 8 depicts the local distribution of the Nusselt number along the pinon row p in10ig. 7 is due to the small uncertainty predicted around the stagnation pointand the downstream wake region. The narrow uncertainty in the predictionfor the upwind values (0 ◦ ≤ β pin < ◦ ) results from small angles betweenthe RI and IP directions in these regions: due to symmetry the angle iszero at β pin = 0 ◦ , and there is no perturbation to the PS directions. Inthe downstream region (135 ◦ < β pin ≤ ◦ ), the result indicates that theperturbations do not strongly affect the heat flux from the bulk flow to thewake region; this deficiency will be a focus of future work. Figure 8: Nusselt number distributions on the pin surface (row
4. Further discussions
The results in Section 3 clearly demonstrate the monotonicity of QoIs(Nusselt numbers) with respect to r RI or r IP , a favorable property for aninterval UQ approach. In this section, we preliminarily analyze the reasonfor this characteristic from one perspective.Figure 9 depicts the case of a parallel shear flow: U x > ∂U x /∂y > ∂ Θ /∂y > y ) di-rection. If the turbulence is close to a state of local equilibrium, we have u (cid:48) x u (cid:48) y < u (cid:48) y θ (cid:48) < u (cid:48) x θ (cid:48) >
0; consequently, the vectors of the two pro-duction terms G I , II i and the BSL predicted PS term Π BSL i can be illustratedas shown in Fig. 9. Given these relative directions, a perturbation of the PS11irection towards the IP limit will increase the magnitude of u (cid:48) y θ (cid:48) , while aperturbation towards the RI limit will decrease the magnitude of u (cid:48) y θ (cid:48) . Sincethis is the flux component directly responsible for the transverse scalar trans-port, and the most relevant component to the QoIs in general, this explainsthe monotonic behavior of the predictions as a function of r RI or r IP . Figure 9: Sketch of a typical scalar flux budget in a parallel shear flow with ∂U x /∂y > ∂ Θ /∂y > Additionally, we note that the flux directions predicted by Π BSL i are usu-ally closer to the RI than the IP direction, which causes the higher sensitivityof the QoIs to r IP than to r RI . This asymmetry is consistent with results ina parallel shear flow: available experimental data for homogeneous shearturbulence [30] indicates that the ratio of ∆ φ RI and ∆ φ IP is approximately1 / ∼ /
5. Conclusions
This paper proposes an approach to quantify scalar flux model uncer-tainties in RANS simulations of turbulent scalar transport. The approachaddresses the inherent inadequacy of the model for the pressure scrambling(PS) term in the transport equations of scalar flux. Specifically, it perturbsthe PS directions predicted by the basic second-moment model within afan-shaped region bounded by the directions corresponding to the return-to-isotropy (RI) and isotropization-of-production (IP) theories.The UQ method is applied to the prediction of heat transfer in a pin-finarray. First, high-fidelity LES data are analyzed to support the proposed12erturbation approach: PS directions estimated through the scalar flux bud-get are shown to concentrate around the RI-IP plane. Next, the resultsdemonstrate that the approach exhibits favorable monotonic behavior: per-turbations of the PS directions towards the IP (or RI) limit consistentlyenhance (or suppress) heat transfer. The resulting plausible intervals for theQoIs, provided by the results of r IP , RI = 0 .
9, encompass the LES predictionsfor the globally averaged and the locally averaged fin Nusselt numbers. Fi-nally, we explain the monotonic behavior of the proposed approach from theperspective of transverse scalar transport in shear flows, which indicates thatthe proposed UQ approach is likely to exhibit monotonic behavior in a widerange of scalar transport problems.To further improve the method and promote its use in engineering appli-cations we identify two future areas of research. First, the analysis of the LESdata indicates that some PS directions are out of the RI-IP plane. As shownin Section 4, perturbations out of this plane seem unnecessary for estimatingplausible intervals for the QoIs in simple cases, but this should be furtherinvestigated for more complex flows. Second, the current method does notintroduce uncertainty in the magnitude of the PS vector. The need to intro-duce this additional uncertainty should be further investigated using analysisof high-fidelity data, or testing of other established PS models (see [31] forreview). To conclude, it is worth noting that the proposed method providesa basis for exploring data-driven approaches where high-fidelity data couldbe used to more accurately inform the PS vector perturbations and reducethe uncertainty in the predictions.
6. Acknowledgements
The research was funded by EUFORIA (grant number IWT-140068).