International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
Contents lists available at SciVerse ScienceDirect
International Journal of Heat and Fluid Flow journal homepage: www.elsevier.com/locate/ijhff
Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry Julia Ling ⇑, Filippo Coletti, Sayuri D. Yapa, John K. Eaton Mechanical Engineering Department, Stanford University, Stanford, CA 94305, USA
a r t i c l e
i n f o
Article history: Received 25 October 2012 Received in revised form 10 June 2013 Accepted 8 July 2013 Available online xxxx Keywords: Turbulent diffusivity Turbulent Prandtl number Turbulent Schmidt number Film cooling Anisotropy
a b s t r a c t A process has been developed by which mean velocity and concentration measurements can be used to determine optimal turbulent diffusivity values for an angled jet in cross-flow configuration. This configuration has applications in film cooling for gas turbine blades. The measurements, obtained by magnetic resonance imaging techniques, provide 3D time-averaged velocity and concentration fields. The mean velocity field is fed into a Reynolds-Averaged Advection Diffusion solver, which uses a turbulent diffusivity model to solve for the mean coolant concentration distribution. This distribution can be compared to the experimentally-obtained concentration field by means of an error metric that quantifies the difference between the computational and experimental concentration fields. By minimizing this error, an optimal value of the turbulent diffusivity can be determined. This optimized distribution is then compared to a RANS simulation to evaluate the relative contribution to error of the turbulent momentum flux model versus the turbulent scalar flux model. Ó 2013 Elsevier Inc. All rights reserved.
1. Introduction Film cooling is a widely used active cooling scheme for gas turbine blades in which coolant air is bled through holes that are distributed over the surface of the blade. In order to predict and improve the film cooling performance, it is important to be able to model the turbulent mixing of the coolant with the main flow. However, traditional Reynolds Averaged Navier Stokes (RANS) simulations often fail to accurately capture this mixing process because of the complex and unsteady nature of the flow and the simplistic nature of standard mixing models (Acharya et al., 2001; Boudier et al., 2007; Hoda et al., 2000). Most commercial computational fluid dynamics (CFD) codes offer a variety of RANS methods for modeling the turbulent momentum flux, including models such as k- and k-x, which employ the Boussinesq turbulent viscosity hypothesis, and Reynolds stress transport models. However, these codes often model the turbulent thermal flux simply through the gradient diffusion hypothesis which stipulates an isotropic turbulent diffusivity. This diffusivity is usually determined by a fixed turbulent Prandtl number, Prt, in accordance with the Reynolds analogy; often Prt = 0.85 is used, based on the empirical value for a flat plate boundary layer (Kays, 1994). However, there is no guarantee that this value for the turbulent Prandtl number will be close to the correct value in more complex flows.
⇑ Corresponding author. Tel.: +1 (650)725 3898. E-mail address:
[email protected] (J. Ling).
A study by Tominaga et al. (2007) investigated several geometries, including plume dispersion in a boundary layer, and showed that the dispersion predictions were sensitive to the turbulent Schmidt number prescribed. Liu et al. (2008) used a k- model to study a film cooling geometry of a single angled jet in cross flow and showed that changing the turbulent Prandtl number had a strong effect on the predicted surface effectiveness. Recently, Schneider et al. (2012) used a Large Eddy Simulation (LES) of a trailing edge film cooling configuration to generate Reynolds stresses and a mean velocity field that were then fed into a RANS solver to calculate the temperature distribution at a fixed turbulent Prandtl number of 0.7. This study showed that even with an accurate mean velocity field and Reynolds stresses, RANS failed to correctly predict the temperature distribution because of the errors incurred by the turbulent mixing model. These studies suggest that it is important to choose the correct turbulent Prandtl number in order to accurately predict the mixing of the coolant flow with the main flow and the resultant surface effectiveness. Otherwise, even if the flow and the Reynolds stresses are correctly predicted, the modeling of the turbulent mixing will be inaccurate. Furthermore, in many flows, the gradient diffusion hypothesis with a fixed Prt is not appropriate. Kaszeta et al. (2000) reported measurements of the eddy viscosity in a film cooling configuration, and showed that the eddy viscosity was significantly anisotropic. Kohli et al. (2005) used cold wire measurements in conjunction with Laser Doppler Velocimetry (LDV) to measure an isotropic but spatially varying turbulent Prandtl number for a configuration with a row of film cooling holes. Their results showed that the
0142-727X/$ - see front matter Ó 2013 Elsevier Inc. All rights reserved. http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
2
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
Nomenclature Prt Sct
vjet vb D u
turbulent Prandtl number turbulent Schmidt number mean jet velocity mean bulk velocity of main flow hole diameter mean velocity
turbulent Prandtl number varied by a factor of 4 spatially. These results highlight the importance of rethinking the uniform isotropic Prt = 0.85 assumption in film cooling configurations. It should be recognized that several models for turbulent mixing have been proposed which are both more complicated and more accurate than the gradient diffusion hypothesis with a set isotropic turbulent Prandtl number. The generalized gradient diffusion hypothesis (GGDH) proposed by Daly and Harlow (1970) is one level up in complexity, but is still often inaccurate in predicting the direction of the scalar flux (Abe and Suga, 2001; Pope, 2000; Tavoularis and Corrsin, 1985). Abe and Suga (2001) proposed an algebraic model for the turbulent scalar flux that uses quadratic products of the Reynolds stress tensor and gives some improvement in the predictions in a wall-shear region of a channel flow. LES have been shown to yield improved accuracy in film cooling modeling (Boudier et al., 2007; Hoda et al., 2000; Tyagi and Acharya, 2003; Peet and Lele, 2008; Renze et al., 2008). However, because of its stability-enhancing properties and its ease of implementation, the gradient diffusion hypothesis is incorporated into many CFD packages and remains in widespread use. Because of the computational cost of LES, RANS is still widely used. It is therefore desirable to work within the framework of the gradient diffusion hypothesis and RANS to optimize the turbulent diffusivity modeling. Experimental measurement of the turbulent diffusivity, especially anisotropic diffusivity, is challenging. There are very few experiments that have measured turbulent diffusivities or turbulent Prandtl or Schmidt numbers in film cooling flows. Kohli et al. (2005) made simultaneous velocity and temperature measurements, enabling them to directly calculate the turbulent heat fluxes in their film cooling geometry. They used the turbulent momentum fluxes and turbulent heat fluxes to calculate the wall-normal turbulent viscosity and thermal diffusivity on the symmetry plane. This method relies on time-resolved, synchronized measurements of the velocity and temperature field at the same point in space. This paper presents a methodology for calculating spatiallyvarying turbulent diffusivity and anisotropic uniform turbulent diffusivity, and determining the sensitivity of the model to the diffusivity value. Furthermore, the paper demonstrates the feasibility of using this technique for validating and optimizing other more complicated turbulent mixing models. This methodology relies on full field mean velocity and concentration measurements obtained using Magnetic Resonance Imaging (MRI) techniques in water. The measurements are carried out in water instead of air because MRI measures the magnetization of the hydrogen spins in water. Because of the fully turbulent and highly 3D nature of this flow, the turbulent mixing dominates over the molecular mixing, so the molecular properties of the fluid will not have a significant effect on the concentration distribution, except at the smallest scales. Neither our experiment nor a typical RANS can resolve these smallest scales. Unlike other techniques that can give the velocity and concentration data in a plane or on a surface, MRI provides fully 3D mean velocity and concentration data. Because MRI measures concentration, not temperature, the passive scalar analogy is
c
at cRAAD cMRC
mean concentration turbulent diffusivity concentration distribution as calculated using RAAD solver concentration distribution as measured using MRC
employed, such that the coolant concentration is measured instead of the fluid temperature. The passive scalar analogy has been leveraged for several previous film cooling studies, including those by Holloway et al. (2002) that measured CO2 concentration at discrete surface taps, and those by Goldstein and Cho (1995) that used the naphthalene sublimation technique. The methodology uses the experimentally measured mean velocity field as an input to a Reynolds-Averaged Advection Diffusion (RAAD) solver, which yields a coolant concentration distribution given an assumed turbulent diffusivity model. This calculated coolant distribution is compared to the experimentally measured mean concentration distribution, and the difference is quantified in an error metric cost function that is minimized to determine the optimal turbulent diffusivity. A critical component of this methodology is identifying a suitable error metric; the process by which the error metric was chosen will be explained and justified. This method has the potential not only to improve RANS modeling of this angled jet in crossflow geometry, but also to serve as a testbed to validate and optimize other turbulent mixing models in different flow configurations. The specific objectives of this research are to (i) develop this new technique using data from an experiment by Coletti et al. (2013) on a single inclined jet in cross flow, (ii) calculate an optimized spatially linearly-varying turbulent diffusivity and an optimized uniform anisotropic turbulent diffusivity from this experimental data, (iii) determine the sensitivity of the calculated coolant concentration distribution to the diffusivity, and (iv) demonstrate the viability of using this technique to investigate other, more complicated turbulent mixing models. The purpose of this paper is not to present a new turbulent diffusivity model but to demonstrate a novel method for calculating the turbulent diffusivity using experimental mean field data. This method represents a significant advancement because it does not rely on simultaneous velocity and concentration data, which are difficult and time-consuming to acquire over a plane, much less over a 3D volume. Instead, it uses mean field data from MRI techniques, which can be acquired in complex flows over a short period of time.
Fig. 1. Schematic of the single inclined jet in crossflow geometry with reference axes. The domain over which the objective function is evaluated is highlighted.
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
2. Optimization procedure 2.1. Configuration The geometry investigated is a single inclined jet in cross-flow shown schematically in Fig. 1. The hole is angled 30° from horizontal, has a diameter D of 5.8 mm and is 24 mm long, and is fed from a rectangular plenum with a volume of 34.5 mL. The test section has a square cross section with side-length 50 mm, and has a boundary layer trip 250 mm upstream of the hole. The nominal main flow rate was 75 liters per minute and the nominal coolant jet flowrate was 0.79 liters per minute. This yields a bulk velocity (i.e. average main flow velocity) of vb = 0.5 m/s and a blowing ratio, vjet/vb, of 1. The Reynolds number based on the bulk velocity and channel hydraulic diameter was 25,000. 2.2. Experimental results The optimization procedure relies on experimental velocity and concentration data obtained using MRI techniques by Coletti et al. (2013). A complete description of the experimental apparatus and measurement procedure is given in Coletti et al. (2013). All measurements were taken in water. Mean velocity data were obtained using Magnetic Resonance Velocimetry (MRV), which gives a full 3D 3-component time-averaged velocity field. The resolution ^ 0:6 mm ^z ð0:11 D ^ ^ was 0.66 mm ^ x 0:6 mm y x 0:10 D y ^, and ^z are the streamwise, wall-normal, and 0:10 D ^zÞ, where ^ x; y spanwise directions respectively. The data matrix of sampled points was 384 162 110. The velocity uncertainty from MRIbased noise was 3% of the mean jet velocity at one standard deviation. Magnetic Resonance Concentration (MRC) was used to obtain the full 3D time-averaged concentration field. The spatial resolution was identical to that of the velocity data with a data matrix of 384 162 104. The uncertainty associated with the MRC measurements is 2.5% at one standard deviation. 2.3. Optimization methodology The optimization relies on experimental mean velocity data as an input to a RAAD solver, which solves for the coolant concentration distribution on the MRV data grid. The RAAD equation, shown in Eq. (1), requires a closure of the turbulent flux term. This closure was accomplished by means of a turbulent diffusivity assumption based on the gradient diffusion hypothesis, as shown in Eq. (2). Because the turbulent diffusivity is much greater than the molecular diffusivity in this flow, the molecular diffusivity can be ignored, giving Eq. (3).
c ¼ ar2 c r u0 c0 ru u0 c0 ¼ at rc c ¼ r ðat rcÞ ru
ð1Þ
3
turbulent thermal diffusivity. Once the mean coolant concentration field is calculated, it is compared to the experimental mean concentration field, and the difference between these two fields is quantified in a cost function. Different turbulent diffusivity values would yield varying cost function values, and the optimal turbulent diffusivity value corresponds to the minimum cost function value. This optimization process is shown schematically in Fig. 2. This procedure thus chooses the values of the parameters in a turbulent diffusivity model that would give the most accurate possible concentration distribution when used with the correct velocity field in the RAAD equation.
2.4. RAAD solver The RAAD solver employed in the optimization process is an inhouse code that reads in the MRV data and defines a mesh with 2.0 million cells with each of the data points at the center of a cell. The cells are hexahedrons spaced 0.6 mm apart in the wall-normal and spanwise directions and 0.66 mm apart in the streamwise direction, corresponding to approximately 9 cells across the hole diameter. A section of this grid near the hole inlet is shown in Fig. 3(a). The experimental MRC data are interpolated onto the same grid as the MRV data, so that the concentration distribution calculated by RAAD can be compared on a cell-by-cell basis to the MRC data. The underlying software is detailed in (Iaccarino, 2011; Pecnik et al., 2012). The solver uses an implicit finite volume approach and treats the diffusive terms with second order central differencing, and the convective terms with first-order upwinding. The upwind scheme was employed because of its good convergence properties. However, on coarse grids, first order upwind schemes can add artificial diffusivity. This effect was quantified by comparing the results of a second order central difference scheme and the first order upwind scheme on a model problem with similar grid refinement and velocity magnitude. The cost function was evaluated for a range of isotropic uniform turbulent diffusivity values. The optimal values of diffusivity were found to vary by less than 15% between the two discretization schemes, and the shape of the cost function was unchanged. Since this variation is on par with the optimization sampling resolution, the upwind scheme was adopted because of its much more favorable convergence properties. For the RAAD solver, the main flow inlet was defined 35 mm (6 hole diameters) upstream of the hole exit, and the concentration at this inlet was prescribed to be 0, since no coolant is expected to flow backwards to the plane of the inlet. The coolant inlet was set 2 mm below the exit plane of the hole, where the concentration was set to be 1. This condition is reasonable because at a blowing ratio of 1, significant main flow ingestion into the hole is not expected.
ð2Þ ð3Þ
In these equations, c is the coolant concentration, u is the velocity vector, and at is the turbulent diffusivity. Overbars indicate a time average and bold indicates a vector. In the most general turbulent diffusivity model, at is a spatially-varying tensor. In the isotropic case, at is scalar; in the spatially uniform case, at is a constant. In this paper, two different optimizations are presented: a uniform anisotropic case in which at is a constant diagonal tensor and a spatially varying isotropic case in which at is a scalar function in space. With the experimental mean velocity field as an input, only the turbulent diffusivity needs to be prescribed in Eq. (3). Because the passive scalar analogy for temperature is used, the optimized quantity is the turbulent scalar diffusivity rather than the
Fig. 2. Schematic of the optimization process.
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
4
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
Fig. 3. Sections of the meshes used for the (a) RAAD and (b) RANS solvers.
The outlet was 152 mm (26 hole diameters) downstream of the hole exit. Because the concentration distribution was no longer changing significantly in the streamwise direction by this plane, the outlet was given a Neumann boundary condition of zero concentration gradient in the streamwise direction. A Neumann condition was also imposed at the walls of the channel, since there is zero diffusion into the wall. This is analogous to an adiabatic condition in a heat transfer setting. All computations were run on a single node that had 24 hyperthreaded cores, Intel Westmere-EP processors, and 36 GB of RAM. 2.5. Cost function The cost function is the means of quantifying the difference between the experimental and computational concentration distributions; in the optimal case, the cost function will be minimized. In order to perform a meaningful optimization, it is critical to choose the cost function wisely. In the present case, there are three main conditions that the cost function must satisfy. First, the cost function has to be sensitive to changes in the optimization variable, turbulent diffusivity. Second, it has to accurately reflect the difference between the concentration field predicted by RAAD and the experimental MRC concentration field. Third, it has to be robust within the experimental uncertainty, and relatively insensitive to experimental outliers. There are two main aspects of the cost function: the error metric and the domain. The error metric is the means of integrating the differences between the RAAD and experimental concentration fields to arrive at a single scalar error value for each test. The domain is the volume over which the error metric is applied. Several error metrics were considered, including the mean difference Emean, the L2 norm EL2, the L1 norm EL1, and the Huber penalty EH. These are defined below, where the sums are over all the cells included in the domain, and n is the total number of cells in the domain:
1X Emean ¼ ðcRAAD cMRC Þ n X 1=2 1 ðcRAAD cMRC Þ2 EL2 ¼ n 1X jcRAAD cMRC j EL1 ¼ n X 1 EH ¼ /ðcRAAD cMRC Þ; n( x2 ; jxj <¼ 0:2 /ðxÞ ¼ 0:2ð2jxj 0:2Þ; otherwise
ð4Þ ð5Þ ð6Þ ð7Þ ð8Þ
The disadvantage of using the mean difference Emean is that error cancelation can occur, where regions of positive error can cancel regions of negative error, thereby giving an artificially low cost function value. The other three error metrics all prevent such
cancelation by means of either squaring or taking the absolute value of the error. They differ in how they weight errors of different magnitudes. The L2 norm, which would give a least-squares optimization scheme, is weighted towards larger errors because it squares the errors before summing. Therefore, the L2 norm is strongly affected by outliers. The L1 norm, on the other hand, weights the larger magnitude errors with much less severity, and more strongly penalizes smaller errors. The Huber penalty is designed to be a balance between the L1 and L2 norms; it does not penalize the outliers as heavily as the L2 norm or the small errors as heavily as the L1 norm because it squares the smaller errors and takes the absolute value of the larger errors (Boyd and Vandenberghe, 2004). This property has made the Huber penalty popular for optimizations that potentially have large outliers as well as Gaussian noise. The disadvantage of the Huber penalty is its added complexity as compared to the L1 or L2 norms. The L1 norm was selected for this study because of its robustness against large errors caused by experimental outliers and its simplicity. It also exhibited good agreement with the L2 and Huber metrics in determining a minimum of the cost function. There are also several possible choices of domains over which to evaluate the error metrics. One obvious choice would be the whole region over which we have MRC data, which is a 26 cm long section of the test channel. The disadvantage of using the whole region is that the expected coolant concentration is zero over a large region of the domain. Small errors in the MRC measurements accumulate over the zero concentration regions, effectively overwhelming the error due to turbulent mixing modeling. Alternatively, the domain could be the region relevant to calculating surface effectiveness: the cells nearest the bottom wall, downstream of the hole and within 1.5 hole diameters on either side of the centerline. However, the near-wall region has the highest uncertainty for MRC measurements because of partial volume effects, in which a volumetric pixel, or ‘‘voxel’’, falls partially within a solid region and partially within a region of flow, thereby giving an artificially low concentration measurement. For this reason, the near-wall region is not an ideal domain. The domain that was chosen was the ‘‘near-jet’’ region, as shown in Fig. 1. It extends in ^ x from immediately upstream of ^ from the hole to 20 hole diameters downstream of the hole, in y 1 voxel above the bottom wall to 2 hole diameters above the bottom wall, and in ^z 1.5 hole diameters on either side of the centerline. This region was chosen because, based on the experimental results, it is the region in which significant coolant concentration is expected and in which the relative uncertainty of the experimental MRC data is expected to be low. In order to evaluate the chosen cost function, its results were compared against the eight other possible pairings of the three error metrics (Emean was excluded because it suffers from error cancelation) and three domains suggested above in an isotropic uniform turbulent diffusivity optimization. The concentration
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
5
Fig. 4. Comparison of potential cost functions values for an isotropic turbulent diffusivity optimization. Cost function values have been normalized to lie between 0 and 1 for each domain. The L1 norm in the near jet domain was 108% higher at at = 1E5 m2/s than at the minimum at at = 3E5 m2/s.
distribution was calculated for a constant scalar at varying between 1E5 m2/s and 5E5 m2/s, and all nine possible cost functions were evaluated at each value of at. The results are shown in Fig. 4, with each proposed cost function normalized over each domain so that it ranges between 0 and 1. All of the cost functions tested are minimized at the same value of isotropic diffusivity: at = 3E5 m2/s and the L1 norm performs similarly to the Huber penalty in all three domains. This agreement between all nine cost functions tested confirms that optimizations carried out with the chosen cost function should be robust with respect to small changes in the error metric and domain. It is also important to quantify the effect of the uncertainty in the experimental velocity field on the cost function evaluation. Therefore, a Monte Carlo uncertainty analysis was employed: the velocity at each point was randomly perturbed within its uncertainty for 1000 trials and was then fed into the RAAD solver with a fixed isotropic uniform turbulent diffusivity of 3E5 m2/s. The cost function was evaluated for each trial and the standard deviation of the cost function value was found to be only 0.6% of its unperturbed value. Therefore the cost function is also robust with respect to the uncertainty in the velocity field. 2.6. RANS method It is instructive to compare the results from the optimized RAAD simulation to a RANS model of this configuration in order to differentiate between errors incurred by modeling the velocity field versus modeling the turbulent scalar fluxes. The main distinction between the RANS solver and the RAAD solver is that the RANS solver must also calculate the mean velocity field, whereas the RAAD solver was given the experimentally measured velocity field as an input. Therefore, whereas the RANS solver will have errors due to both an incorrect velocity field and turbulent diffusivity modeling, the RAAD solver has errors only due to the latter effect. The RANS simulation was run using FLUENT 13.0 standard k- with enhanced wall treatment and a coupled second-order solver. This turbulence model was chosen because of its prevalent usage in both academic and industrial applications for film cooling geometries. Furthermore, Ferguson et al. (1998) compared the performance of the standard k- model without wall functions to the performance of the renormalization group (RNG) k- model and a Reynolds Stress Model (RSM) in a film cooling configuration and showed that the standard k- model was able to more closely match experimental data. Karvinen and Ahlstedt (2005) compared the performance of various turbulence models in a jet-in-crossflow configuration, and showed that standard k- and k-x SST gave the best agreement with experimental results in terms of their prediction of the velocity profiles and turbulent kinetic energy profiles.
However, both of these reports agreed that none of the turbulence models employed matched experimental data in all regions of the flow. The inlet of the main channel was set 87 mm (15 D) upstream of the hole leading edge, which was 117 mm (20.2 D) downstream of the boundary layer trip. The outlet was 203 mm (35 D) downstream of the hole. The grid consisted of 2.5 million hexahedral cells that were meshed using a Cooper scheme, with 140 cells around the circumference of the hole. A section of this mesh is shown in Fig. 3(b). Boundary layer meshes were applied to both the bottom wall of the channel and the inside of the hole with the first row located at y+ < 1. When the simulation was converged, all residuals were less than 106. Grid-independent convergence was confirmed by running the same simulation on a grid with 4.6 million cells. The turbulent viscosity varied on average by less than 1% between the results on the two grids in the ‘‘near jet’’ domain. The velocity inlet condition for the main channel was prescribed using experimental MRV data. The main inlet turbulence boundary conditions were prescribed based on empirical relations from Moin and Kim (1982) and Pope (2000) and an estimated boundary layer thickness of 6 mm based on experimental velocity profiles. Profiles of the prescribed turbulent kinetic energy and turbulent dissipation are shown in Fig. 5. These prescribed inlet conditions let to turbulence intensities of approximately 4.6% at the main inlet mid-channel. The velocity inlet for the coolant plenum was set at a uniform velocity of 0.079 m/s (0.16 vjet) based on the
Fig. 5. Turbulence inlet profile for the main channel for RANS.
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
6
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
experimentally measured coolant flowrate, and the turbulence inlet conditions were also prescribed uniformly as 1% turbulent intensity and 1 mm turbulent length scale. Because the flow inside the coolant plenum is mostly stagnant, the velocity field was found to be much less sensitive to the coolant flow inlet conditions than the main flow inlet conditions. 3. Optimization of two-parameter turbulent diffusivity models In order to validate our optimization process, two different simple two-parameter models were chosen. The first was anisotropic uniform turbulent diffusivity, where ay and az were the two parameters and ax was set to be their average. The second was an isotropic turbulent diffusivity that varies linearly in the streamwise direction, and is constant in the other two directions: at = ax + b, where a and b are the two parameters to be optimized. In both cases, a brute force exhaustive search approach was employed, in which a full matrix of parameters was tested. This approach was chosen to determine the sensitivity of the objective function to the two parameters and to better understand the shape of the cost function over the range of values. However, because the brute force method was used, we were limited to two-parameter models. It would have also been possible to sacrifice information on the shape of the cost function and use a more efficient optimization algorithm for a model with more parameters. These two-parameter models are not meant to be suggestions for realistic turbulent scalar flux models, but rather serve as case studies illustrating the versatility of the optimization process and how the sensitivity of the concentration distribution to the turbulent diffusivity model can be analyzed by examining the cost function. 3.1. Anisotropic uniform turbulent diffusivity The first turbulent diffusivity model that was optimized was anisotropic, uniform diffusivity. In this model, the turbulent diffusivity is a constant diagonal tensor. Because the convective term in ^ x is much greater than the diffusive term, the primary parameters ^ and ^z: ay and az. of interest were the turbulent diffusivities in y These two parameters were optimized over a 7 by 7 matrix of diffusivity values, ranging from 1E5 m2/s to 5E5 m2/s for ay and az. At each combination, ax was set to be the average between ay and az. The 49 cases in the 7 by 7 matrix took approximately 1 week of run time on a single node. The results of this anisotropic diffusivity optimization are shown in Fig. 6. The minimum value of the L1-norm ‘‘near-jet’’
Fig. 6. Results of the anisotropic uniform diffusivity optimization. Contours of the cost function for varying ay and az.
cost function is 0.01589, and occurs at ay = 2.5E5 m2/s and az = 3.5E5 m2/s. This cost function minimum value means that in the optimum case, the average difference magnitude between the calculated and experimental concentration distributions is 1.589% concentration in the near-jet domain. It should also be noted that even a ‘‘perfect’’ concentration distribution would not achieve a cost function evaluation of zero because of noise in the experimental data. At this optimum point for ay and az, a 1-parameter optimization was then run in ax. In the range of 1E6 m2/s to 5E5 m2/s, the cost function was minimized at ax = 1E6 m2/s, which is the molecular viscosity. The optimization in ay and az was then re-run at this optimal value of ax to see if the minimum point moved in ay or az. It was determined that neither the location of the minimum nor the shape of the cost function changed significantly for ax = 1E6 m2/s. The new minimum value of the cost function was 0.01578, a change of less than 1% from the previous minimum, indicating that the cost function was relatively insensitive to ax, as would be expected. This result is interesting in that it indicates a significant anisotropy, with 40% higher values of turbulent diffusivity in the spanwise direction than in the wall-normal direction. It should be noted that the anisotropy of the flow is most likely not uniform. Nevertheless, these results demonstrate that the turbulent scalar fluxes are anisotropic in this flow, with on average greater transport in the spanwise direction than in the wall-normal direction. Kaszeta et al. (2000) measured the anisotropic turbulent viscosity in a single hole film cooling geometry and also reported generally larger values in the spanwise direction than in the wall-normal direction. The shape of the cost function is important for two key reasons. The first is that it indicates that within this range of turbulent diffusivity values, there is a single minimum – there are no local minima. This is both important and encouraging, as it greatly reduces the computational effort that would be required to optimize the cost function using a more efficient method to sample points in ay and az. The second crucial piece of information provided by the cost function surface is the sensitivity of the concentration distribution to the turbulent diffusivity. For example, it can be seen that an increase by 40% of ay or az from their optimal values yields an increase in the cost function by 11.1% or 21.9% respectively. A decrease by 40% of ay or az causes an increase in the cost function by 12.3% or 21.9% respectively. The concentration distribution is therefore most sensitive to variation in az. Furthermore, while the cost function is approximately equally sensitive to positive and negative perturbations of this magnitude in az, it is more sensitive to negative perturbations in ay than positive ones. This indicates that in this flow it is ‘‘safer’’ to over-estimate ay than under-estimate it in order to recover the correct concentration distribution. It is instructive to compare the concentration distribution from RAAD with the optimized turbulent diffusivity to the experimental results from MRC to determine in which regions RAAD is failing to predict the correct distribution. Fig. 7(a and b) shows the experimentally-measured concentration distribution in two planes: the center plane of the channel in the spanwise direction and a wall-parallel slice 1 mm (0.2 D) above the bottom wall. Fig. 7(c and d) shows the difference between the computed and experimental concentrations in the same two planes. As can be seen from this figure, the match is not entirely satisfactory, though the difference is less than 20% of the absolute concentration over most of the test volume. It should be remembered that very near the walls, and in the hole, the experimental MRC data are prone to error due to partial volume effects. Furthermore, because of the finite resolution of the grid, the exact shape of the hole could not be captured, leading to increased error near the hole exit.
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
7
Fig. 7. Subfigures (a and b) show contours of the experimentally-measured MRC concentration distribution. Subfigures (c and d) show contours of the difference between the RAAD concentration field for the optimized anisotropic uniform diffusivity and the MRC concentration field.
There are several features in Fig. 7 that relate to the inadequacy of the uniform, anisotropic model. Fig. 7(b and d) shows two parallel lobes just downstream of the hole where the experimental results show high concentration but where RAAD underpredicts the concentration. These lobes correspond to the region occupied by the counter-rotating vortex pair where the concentration would be 100% if there were no diffusivity; a deficit of coolant in this region corresponds to too high of a turbulent diffusivity. Conversely, in the region between these two lobes, RAAD overpredicts the coolant concentration. In essence, too high of a turbulent diffusivity is bleeding coolant from the two lobes and into the center region, causing the observed error pattern. This error arises because of local deviation away from the globally optimized turbulent diffusivity values, indicating the insufficiency of a uniform turbulent diffusivity model. Either a spatially-varying turbulent diffusivity or a more complex model for the turbulent scalar fluxes is required. 3.2. Linearly-varying isotropic turbulent diffusivity The second turbulent diffusivity model that was optimized was isotropic turbulent diffusivity that varies linearly in the streamwise ð^ xÞ direction. In this case, the turbulent diffusivity is a scalar function of streamwise location. This model was chosen because of its simplicity (it depends only on two parameters), and because it seemed likely that the turbulent diffusivity would vary in the streamwise direction. The two parameters in the optimization are the slope, a, and the offset, b, such that at = ax + b where x = 0 at the leading edge of the hole. The optimization was run over a 6 by 5 matrix in which a was varied in the range [1E4 m/s, 4E4 m/s], and b was varied in the range [1E5 m2/s, 5E5 m2/s]. A slope of a = 1E4 m/s corresponds to a variation in the diffusivity of 1.5E5 m2/s between the hole exit and the end of the domain. These ranges were chosen after sampling a broader range of values to determine the approximate ranges in which the minimum would occur. Fig. 8(a) shows the results of this optimization. Although the turbulent diffusivity only varied in one direction, it should be noted that the cost function was still evaluated over the full 3D region of interest in the near jet domain, since variation of the turbulent diffusivity in any direction affects the full concentration distribution. The cost function obtained its minimum value of 0.01613 at a = 2E4 m/s and b = 2E5 m2/s. This represents only a 1% decrease in the cost function from the optimal spatially-uniform case of a = 0 and b = 3E5 m2/ s. The case of a = 1E4 m/s and b = 1E5 m2/s did not converge because it gave negative
Fig. 8. Results of the linearly-varying isotropic diffusivity optimization.
diffusivity over a large portion of the flow. The cost function is therefore minimized for a turbulent diffusivity that increases in
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
8
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
the streamwise direction. This rise in turbulent diffusivity in the streamwise direction could be due to growing vortex length scales. What is most interesting about the shape of the cost function is that there is a compensating effect between the slope and offset: a lower slope with a higher offset performs comparably to a higher slope with a lower offset. This behavior suggests that there could be a critical region where it is most important to have the correct diffusivity. There were 7 combinations of a and b that had cost function values within 10% of the minimum. The profiles of turbulent diffusivity for these 7 combinations are shown in Fig. 8(b). As can be seen from this plot, all 7 combinations have similar values in the region 5–10 D downstream of injection, suggesting that this is the region in which it is most critical to have the correct turbulent diffusivity. Mahesh (2013) characterized a jet-in-crossflow as having two distinct regimes: a near field that is pressure-driven and a far-field that is momentum-driven. Perhaps the region 5–10 D downstream of injection is critical because it marks the transition to the far field regime. 4. Comparison to RANS model It is interesting to compare the results from the optimized RAAD case to the results from the RANS calculation running standard k-. Unlike the RAAD solver, the RANS solver must also calculate the velocity field and model the turbulent momentum fluxes. Since RANS calculates the turbulent viscosity, it can implement either a fixed Sct model or a fixed at model, allowing for a comparison between these two models. Therefore, 7 different fixed values of Sct and at were employed, ranging from 0.5 to 2.0 and 1E5 m2/s to 7E5 m2/s respectively. In all cases, the turbulent diffusivity was isotropic. Fig. 9 shows the turbulent diffusivity calculated by RANS on the channel midplane for Sct = 0.85, along with the corresponding 10% concentration isocontour. The turbulent diffusivity ranges between 2E5 m2/s and 6E5 m2/s in the region occupied by the jet. In order to compare the results of the RANS to that of RAAD, the RANS results were linearly interpolated onto the MRV experimental data grid, and then the cost function was evaluated for the RANS results. Fig. 10 shows the cost function values for the varying values of Sct and at, with a horizontal line indicating the minimum error achieved using RAAD with an isotropic uniform turbulent diffusivity. The minimum values of the cost function occur at Sct = 0.85 and at = 5E5 m2/s. Interestingly, the minimum value of the cost function for a fixed at is only 6% higher than that with a fixed Sct, indicating that in this case, a fixed Sct is not a much more accurate model than a completely uniform at. Perhaps the similarity in the performance of these two models is due to the relatively uniform nature of the turbulent viscosity predicted by RANS in the near-jet region, as seen in Fig. 9. Fig. 11 shows contours of the difference between the RANS concentration field for Sct = 0.85 and the MRC concentration field.
Fig. 9. Contours of turbulent diffusivity in m2/s on the channel centerplane as predicted by RANS for Sct = 0.85. The 10% concentration isocontour is shown as a black line.
Fig. 10. Results of RANS optimization for varying fixed turbulent diffusivity (solid —) or Schmidt number (dashed - - -). Optimal value for isotropic uniform RAAD optimization shown for comparison (dotted . . .).
Although globally, the optimal uniform anisotropic diffusivity in RAAD out-performs the optimal turbulent Schmidt number in RANS, there are still local regions (e.g. in the two lobes of high concentration immediately downstream of injection) where RANS is more accurate. These are regions where the local turbulent diffusivity is highly non-uniform, and where the uniform anisotropic assumption is therefore least appropriate. However, overall, the regions of high error are similar between the best anisotropic uniform diffusivity case for RAAD and the best Sct case for RANS. Both simulations over-predict the centerline concentration immediately downstream of the hole in the separation region and underpredict the concentration on either side. Both codes also over-predict the concentration just upstream of the hole. These regions can therefore be identified as the zones where the gradient diffusion hypothesis is invalid. It is possible that the breakdown of the gradient diffusion hypothesis in these regions is based on features of the turbulence structure, a discussion of which is beyond the scope of this study. By comparing the RANS results with the RAAD results, it is possible to examine the effects of an incorrect velocity field on the optimal turbulent diffusivity values. The optimal Sct and at from RANS gave cost function values that were 18% and 25% higher respectively than the cost function value for the optimal isotropic uniform turbulent diffusivity value from RAAD. Therefore, even with a much higher grid resolution near the hole, RANS was not able to match the performance of a solver with the correct velocity field (within experimental uncertainty). Furthermore, the location of the cost function minimum for an isotropic uniform diffusivity was different between RANS and RAAD: RANS minimized the cost function at at = 5E m2/s and RAAD minimized the cost function at at = 3E5 m2/s. This difference in optimal turbulent diffusivities was due to RANS slightly over-estimating the streamwise velocity in the near-wall region; as a result, higher turbulent diffusivity was required to compensate to match the experimental concentration field. It is also possible to roughly estimate the relative contributions of turbulent momentum flux modeling and turbulent scalar flux modeling to errors in the concentration distribution. A precise evaluation of the relative contributions to the error is not available at this time because the RANS and RAAD models were run on different meshes and because of the contribution of experimental noise to the cost function. However, the contribution of experimental noise to the cost function can be roughly estimated by
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
9
Fig. 11. Contours of the difference between the RANS concentration field for Sct = 0.85 and the MRC concentration field.
looking at the mean value of the experimentally measured concentration magnitude in a region where zero concentration is expected (e.g. upstream of the jet). This mean value was approximately 0.0076, which suggests that a correct concentration distribution would give a cost function value of 0.0076, and that the cost function evaluations are therefore offset by this amount. If 0.0076 is subtracted from the optimal values of the cost function for an isotropic uniform turbulent diffusivity for RAAD and RANS, then the cost function values become, respectively, 8.65E3 and 1.27E2. Because RAAD only has to model the turbulent scalar flux, but RANS must model both the turbulent scalar flux and the turbulent momentum flux, it can be deduced that the ratio of these two errors represents a rough estimate of the contribution of the turbulent scalar flux model to the total error in the concentration distribution. This ratio, 68%, would suggest that the turbulent scalar flux model contributes a significant, and perhaps a majority, of the error to the concentration distribution in this case.
5. Conclusions A new methodology was presented for optimizing turbulent diffusivity models for an inclined jet-in-crossflow using mean fullfield velocity and concentration data. These optimizations yielded information on in which regions it is most critical to have the correct diffusivity value, what the sensitivity of the concentration distribution is to these diffusivity values, and how robust the model is to perturbations in the diffusivity. By comparing to RANS, it was possible to assess the relative merits of a fixed at model versus a fixed Sct model and to evaluate the performance of RANS versus RAAD. By examining in which regions both RANS and RAAD incorrectly predicted the concentration distribution, it could be concluded in which regions a gradient diffusion assumption was inappropriate. Lastly, it was possible to give a rough evaluation of the relative contribution of the turbulent scalar flux model to the total error in the concentration distribution. This estimation showed that the turbulent scalar flux model is a major contributor to the total error in the concentration distribution. Several extensions to this methodology are possible. This methodology can be applied to more complex flow configurations and more realistic turbulent diffusivity models. The turbulent diffusivity could be allowed to be a function of the mean velocity gradient tensor or other higher parameter spaces, and the optimization could be performed using more efficient algorithms (e.g. genetic algorithms, simulated annealing, etc.) with the trade-off of losing information on the shape and sensitivity of the cost function. It would also be valuable to compare the cost function shapes of several different film cooling configurations to determine how sensitive the cost function is to geometric variation. The versatility of this technique stems from the fact that it only requires volumetric mean field measurements, which can be provided efficiently by MRI techniques.
Acknowledgments J. Ling was supported by an NSF Graduate Research Fellowship. The authors acknowledge the following award for providing computing resources that have contributed to the research results reported within this paper: MRI-R2: Acquisition of a Hybrid CPU/ GPU and Visualization Cluster for Multidisciplinary Studies in Transport Physics with Uncertainty Quantification, Award Number 0960306. This award is funded under the American Recovery and Reinvestment Act of 2009 (Public Law 111-5). The authors also gratefully acknowledge the guidance and help of G. Iaccarino and C. Gorle. References Abe, K., Suga, K., 2001. Towards the development of a Reynolds-averaged algebraic turbulent scalar-flux model. Journal of Heat and Fluid Flow 22 (February), 19– 29. Acharya, S., Tyagi, M., Hoda, A., 2001. Flow and heat transfer predictions for film cooling. Annals of the New York Academy of Sciences 934 (May), 110–125. Boudier, G., Gicquel, L., Poinsot, T., Bissieres, D., Berat, C., 2007. Comparison of LES, RANS and experiments in an aeronautical gas turbine combustion chamber. Proceedings of the Combustion Institute 31 (January), 3075–3082. Boyd, S., Vandenberghe, L., 2004. Convex Optimization. Cambridge University Press, pp. 294–301. Coletti, F., Benson, M., Ling, J., Elkins, C., Eaton, J., 2013. Turbulent transport in an inclined jet in crossflow. The International Journal of Heat and Fluid Flow, http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.06.001. Daly, B., Harlow, F., 1970. Transport equations in turbulence. Physics of Fluids 13, 2634–2649. Ferguson, J., Walters, D., Leylek, J., June 1998. Performance of turbulence models and near-wall treatments in discrete jet film cooling simulations. In: ASME International Gas Turbine and Aeroengine Congress and Exhibition 43rd. Goldstein, R., Cho, H., 1995. A review of mass transfer measurements using naphthalene sublimation. Experimental Methods in Thermal and Fluid Science 10 (May). Hoda, A., Acharya, S., Tyagi, M., 2000. Reynolds stress transport model predictions and large eddy simulations for film coolant jet in crossflow. ASME Turbo Expo 2000 (May). Holloway, D., Leylek, J., Buck, F., 2002. Pressure-side bleed film cooling: Part I steady framework for experimental and computational results. ASME Turbo Expo 3 (June). Iaccarino, G., 2011. Stanford University ME469 Course Notes. Karvinen, A., Ahlstedt, H., 2005. Comparison of turbulence models in case of jet in crossflow using commericial CFD code. Engineering Turbulence Modelling and Experiments 6, 399–408. Kaszeta, R., Simon, T., 2000. Measurement of eddy diffusivity of momentum in film cooling flows with streamwise injection. Journal of Turbomachinery 122, 178– 183. Kays, W., 1994. Turbulent Prandtl number – where are we? Journal of Heat Transfer 116 (May), 284–295. Kohli, A., Bogard, D., 2005. Turbulent transport in film cooling flows. Journal of Heat Transfer 127, 513–520. Liu, C.-L., Zhu, H.-R., Bai, J.-T., 2008. Effect of turbulent Prandtl number on the computation of film-cooling effectiveness. International Journal of Heat and Mass Transfer 51, 6208–6218. Mahesh, K., 2013. The interaction of jets with crossflow. Annual Review of Fluid Mechanics 45, 379–407. Moin, P., Kim, J., 1982. Numerical investigation of turbulent channel flow. Journal of Fluid Mechanics 118, 341–377. Pecnik, R., Terrapon, V., Ham, F., Iaccarino, G., Pitsch, H., 2012. Reynolds-averaged Navier–Stokes simulations of the HyShot II scramjet. AIAA Journal 50, 1717– 1732. Peet, Y., Lele, S., 2008. Near field of film cooling jet issued into a flat plate boundary layer: LES study. ASME Turbo Expo 2008 (June).
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005
10
J. Ling et al. / International Journal of Heat and Fluid Flow xxx (2013) xxx–xxx
Pope, S., 2000. Turbulent Flows. Cambridge University Press, New York, USA. Renze, P., Schroder, W., Meinke, M., 2008. Large-eddy simulation of film cooling flows at density gradients. International Journal of Heat and Fluid Flow 29 (February), 18–34. Schneider, H., Bauer, H., von Terzi, D., Rodi, W., 2012. Coherent structures in trailing-edge cooling and the challenge for turbulent heat transfer modelling. ASME Turbo Expo 2012 (June).
Tavoularis, S., Corrsin, S., 1985. Effects of shear on the turbulent diffusivity tensor. Journal of Heat and Mass Transfer 28 (January), 265–276. Tominaga, Y., Stathopoulos, T., 2007. Turbulent Schmidt number for CFD analysis with various types of flowfield. Atmospheric Environment 41 (December), 8091–8099. Tyagi, M., Acharya, S., 2003. Large eddy simulation of film cooling flow from an inclined cylindrical jet. Journal of Turbomachinery 125, 734–742.
Please cite this article in press as: Ling, J., et al. Experimentally informed optimization of turbulent diffusivity for a discrete hole film cooling geometry. Int. J. Heat Fluid Flow (2013), http://dx.doi.org/10.1016/j.ijheatfluidflow.2013.07.005