Large Eddy Simulation of Turbulent Flow and Heat Transfer in a Ribbed Coolant Passage

Numerical simulations of hydrodynamic and thermally fully developed turbulent flow are presented for flow through a stationary duct with periodic array of inline transverse rib turbulators. The rib height to hydraulic diameter ratio e/Dh is 0.1 and the rib pitch to rib height ratio P/e is 10. The effect of secondary flow due to presence of rib turbulators on heat and mass transfer has been investigated. The present work reviews the use of a large eddy simulation LES turbulence model, known as shear-improved Smagorinsky model SISM , for predicting flow and heat transfer characteristics in the fully developed periodic flow region. The computations are performed for Reynolds number of 2,053 and the working fluid chosen to be air, the Prandtl number of which is 0.7. Instantaneous flow field, time-mean, and turbulent quantities are reported together with heat transfer and a close match with experiments has been observed.


Introduction
The gas turbines are used for aircraft propulsion, land-based power generation, and industrial applications.The rotor inlet temperature is one of the important parameters being used for raising gas turbine engine performance.The increase in inlet temperatures results in increased temperature of turbine blades.Advanced gas turbines operate at temperatures well above the permissible melting point of the material of blades.This compels active cooling of turbine blades in order to prevent breakdown.Blades are typically cooled by air extracted from the compressor of the engine thereby incurring a penalty to thermal efficiency.This makes it utmost necessary to understand and optimize the cooling techniques.
Blades are cooled internally by passing the bleed-off air through ribbed internal passages.Laminar sublayer is found to offer thermal resistance to convective heat transfer from the thermally active surfaces.Ribs act as turbulent promoters by breaking the laminar sublayer and creating local wall turbulence because of flow separation from the ribs and reattachment, and they enhance momentum as well as heat transfer.The escalation in heat transfer is also accompanied by increase in resistance to the flow field.The presence of and turbulent heat transfer in a stationary square duct with transverse and angled rib turbulators.The Reynolds number in the study was low and no comparison with experimental data was provided.Effect of centrifugal buoyancy on turbulent heat transfer was also reported by Murata and Mochizuki 18 .The effect of centrifugal buoyancy, was found to increase and decrease the heat transfer for buoyancy induced aiding and opposing flows, respectively.
Ooi et al. 19 compared the performance of one equation Spalart-Allmaras and v 2 − f turbulence models with the standard k −ε model.They reported that heat transfer predictions obtained by using v 2 − f model closely with the experimental data than the other two models used in the study.Watanabe and Takahashi 20 performed LES for fully developed stationary ribbed duct flow.The duct aspect ratio studied was 0.5 and the parameters e/D h and P/e were 0.1 and 10, respectively.The Reynolds number based on bulk velocity was 117,000.Experimental results were also reported in their work.They observed a close match between the experiments and numerical results.Their results also reveal the unsteady mechanism that is responsible for enhancement of the heat transfer coefficient.Abdel-Wahab and Tafti 21 reported the effect of Coriolis and centrifugal buoyancy forces on heat transfer in a 90 • ribbed duct with rotation using LES.They analyzed that the reattachment length on leading wall is increased by the action of buoyancy forces, up to a certain value of rotation number where the separated shear layer extends all the way to the downstream rib resulting in degradation of heat transfer.However, buoyancy force does not affect reattachment length on trailing wall significantly.On the other hand, the turbulence levels on the trailing wall are amplified significantly by the action of buoyancy forces in addition to that caused by Coriolis forces.This results in heat transfer augmentation on the trailing wall.They also concluded that the buoyancy also results in strengthening of secondary flow in the duct cross-section, bringing about enhanced heat transfer along both the smooth walls.Tyagi and Acharya 22 emphasized the importance of large-scale coherent structures on mixing and consequently heat transfer.Role of subgrid stress modeling on the fluid flow and heat transfer in ribbed ducts for internal cooling was discussed in detail by Tafti 23 . Sewall et al. 24 comprehensively validated the use of LES for predicting flow and heat transfer in a stationary ribbed duct.
The objectives of the present study is to evaluate the use of shear-improved Smagorinsky model SISM , suggested by Lévêque et al. 25 , for predicting flow and heat transfer through a stationary ribbed duct of square cross-section, and elucidate the detailed physics encountered in the fully developed region.The flow and thermal field for the present flow configuration is statistically nonstationary.The shear-improved Smagorinsky model is proven to be effective for statistically stationary flow Toschi et al. 26 and Saha and Biswas 27 .There is a dearth of literatures available for nonstationary flow.This provides the prime motivation to perform the present study.The effect of secondary flow due to rib turbulators on heat and mass transfer in a stationary square ducts has been investigated.The capability of the shear-improved model in predicting turbulent flow field and heat transfer is evaluated.The computations are performed for Reynolds number of 2,053 in a square duct with inline normal ribs with e/D h 0.1 and P/e 10.These parameters give similar configuration to that studied by Murata and Mochizuki 17 .It is to be noted that the flow in a smooth square duct is found to be turbulent at a Reynolds of 2,060 Hanks and Ruo 28 .A lower value of critical Reynolds number is expected for flow through a rib-roughened square duct.Similarly, Nonino and Comini 29 reported a Reynolds number of around 1,000 for chaotic flow inside a square duct populated with periodic array of the ribs arranged in staggered fashion.Consequently, one can expect the flow at Reynolds number of 2,053 to be turbulent.Instantaneous flow field, time mean, and turbulent quantities are reported together with thermal field.

Computational Model and Governing Equations
The computational model used for the present study is shown in Figure 1.The model assumes fully developed fluid flow and heat transfer in periodically repeating spatial structures in a square duct.The working fluid chosen is air Pr 0.7 .
In the literature, as per the author's knowledge, there are two approaches to simulate the fully developed flow using periodic boundary conditions in a domain having geometric periodicity.In the first approach, the simulation is carried out taking the characteristic velocity as the frictional velocity which gives the mean pressure gradient to be unity.The second method is one in which the average velocity of the flow inside the duct is used as the characteristic velocity.In the former approach, since the mean pressure gradient is fixed by the time-averaged frictional velocity which is constant , the inlet average velocity varies with time and for a given frictional Reynolds number, the average Reynolds number based on inlet average velocity fluctuates in an unsteady flow not known a priori and can only be ascertained once the time-averaged field is calculated after the flow reaches the dynamic steady state.On the other hand, in the latter case, since the average Reynolds number is fixed, which is more practical as the mass flow rate is generally given as one of the input parameters, the mean pressure gradient fluctuates to keep the inlet velocity i.e., Reynolds number constant.It is to be noted that for a physically steady flow, the mean pressure gradient comes out to be constant for a particular Reynolds number.The algorithm of the present study is such that the mean pressure gradient β t is iterated to a value that ensures the inlet velocity equal to the constant average velocity corresponding to the Reynolds number.As a consequence, the mean pressure gradient is expected to vary as the flow in the present study is an unsteady one and results in unequal frictional forces from one time to the other.
In the present study, incompressible Navier-Stokes equations along with the energy equation in primitive variable forms are used to simulate fluid flow and heat transfer.As the domain is periodic in streamwise directions, it necessitates the means gradient of pressure be isolated from the fluctuating components which can be expressed as follows: where β t is the linear component of nondimensional pressure and is adjusted at each time step to achieve the desired mass flow rate.On substituting into Navier-Stokes and energy equations, the governing equations in dimensionless form, Saha and Acharya 30 , are expressed as where τ ij u i u j − u i u j is the subgrid scale stress SGS tensor and the subgrid scale energy flux q j −u θ .The unknown function Λ in the energy equation is given by The coupling between θ and Λ is solved iteratively as described by Wang and Vanka 31 .The governing equations are nondimensionalized by using length scale D h , that is, the duct hydraulic diameter, velocity scale u avg , pressure by ρu 2  avg , and time by D h /u avg .Temperature is nondimensionalized by using the scale θ T − T w / T m1 − T w , where T w is the wall temperature and T m1 is the bulk temperature of the fluid evaluated at inlet.Finite-volume method on collocated grid is used to solve 2.2 , 2.3 , and 2.4 .To prevent the pressure velocity decoupling, momentum interpolation method MIM originally proposed by Rhie and Chow 32 has been used.Both convective and diffusion terms are discretized using central difference scheme.Time advancement is achieved semi-implicitly.
The SGS stress tensor in 2.3 is parameterized by eddy viscosity model as The eddy viscosity ν T , as defined by Smagorinsky model, is given by C S Δ 2 |S|, where C S is the Smagorinsky constant with value varying in the range of 0.065-0.25 and . The major drawback of the Smagorinsky model is the introduction of spurious dissipation near the wall damping the small perturbations and thereby transition to turbulence.The subgrid scale energy flux is given by where Pr SGS is the turbulent Prandtl number for subgrid-scale component and its value is set to 0.5, Moin et al. 33 .
No-slip boundary condition u i 0 is imposed on all confining surfaces.The periodic boundary condition generalized by 2.8 is used along the streamwise direction: where φ u i or p, "L x " is the length of domain in streamwise direction, and "n" represents the number of modules included in the computational domain, one for the present study.The corresponding periodic condition for energy equation is where, θ m1 and θ m2 are the bulk temperatures calculated at the inlet and outlet of the computational domain, respectively.Constant wall temperature θ w 0 is imposed on all the confining surfaces.
The geometry of the present study is a square duct with ribs attached to top and bottom surfaces at periodic intervals implying the geometry to be periodic in the streamwise direction.It is generally expected that the flow is periodic when the fully developed state is reached.It may so happen that the flow periodicity and geometric periodicity do not match even in fully developed condition, but, it has been established by Saha and Acharya 34 that the geometric periodicity guarantees flow periodicity at least in nonrotating conditions.Consequently, the simulations are carried out using the streamwise length having one geometric periodicity.This helps in scale-resolved inertial scales computations and leads to reduction in computational cost.For almost similar geometry pin-fin 35 , the study with two different periodic modules has also been carried out but no change in the time-averaged field has been observed.
One can model the periodic module according to easiness in generating grid.In the present study, we have put the ribs in the middle of the domain.One may choose a geometry in which two halves of the ribs are placed at inlet and outlet of the domain.
The heat transfer augmentation ratio Nu/Nu 0 is calculated using the correlation of Nusselt number for smooth duct: Nu 0 0.023 Re 0.8 Pr 0.4 .

2.10
The grid size used to perform the simulation over the domain is 86 × 86 × 86.To resolve the near wall viscous effects, the grid density is kept high near the walls as well as on the rib surfaces.the corresponding values are 0.35 in each direction and are well below 1.0 generally needed for turbulent flow calculations without any wall treatment.The temporal variation of streamwise velocity component and the corresponding power spectra plotted downstream of the rib at x/e 8, y/e 0.5, and z 0.5 and at the center of the computational domain x y z 0.5 have been shown in Figure 2. The power spectra without any energy builtup at higher frequencies at both locations Figure 2

Results and Discussion
The results obtained for the geometry studied, shown in Figure 1, are reported in the present section.The ratios e/D h and P/e are kept fixed at 0.1 and 10, respectively.The computations are performed at Reynolds number of 2,053.Instantaneous flow fields, time mean, and turbulent quantities are presented.

Instantaneous Field
The secondary velocity vectors superimposed on the temperature contours, at x 0.5 and 1.0, are presented in Figure 4.The secondary flow in a smooth square duct is characterized by the presence of eight recirculation regions arranged symmetrically about the diagonals of the square cross-section.The flow downstream of the ribs Figure 4 b shows stronger secondary flow structures than the location at the mid-section of the rib Figure 4 a .Close examination of the secondary flow structures at the two locations shows that more number of secondary vortices are present near the top and bottom wall at x 1.0, which corresponds to the outflow plane.The temperature gradients near sidewalls for both locations are seen to be almost equal.However, the gradient near the top region of the two ribs at x 0.5 is found to be higher compared to corresponding near wall region top and bottom walls at x 1.0.This is because of the fact that the streamwise bulk velocity is quite high near the top region of two ribs at x 0.5 in comparison to the location x 1.0.Figure 5 presents the instantaneous contours of normalized Nusselt number and skin friction coefficient scaled with Re on one of the sidewalls and one of the ribbed-roughened walls.The regions of higher heat transfer are visible on either walls but the highest heat transfer has been observed on the rib surface.The flow reattachment at the isolated locations may be the reason for higher heat transfer on both the rib surface and the flat bottom wall.As expected, the strong recirculation region downstream of the rib results in poor heat transfer and is obvious from the contours of the Nusselt number.Similarly, skin friction distribution shows higher wall shear stress on rib surface and the sidewall where flow gradient is higher.
Figure 6 presents the temporal variation of space-averaged Nusselt number normalized with Nu 0 and skin friction coefficient scaled with Reynolds number for one of the sidewalls and one of the ribbed-roughened walls.Strong unsteadiness in the Nusselt number and skin friction are clearly brought out in the two plots.Close examination between the Nusselt number and skin friction at both walls reveals a close correlation between the Nusselt number and skin friction.The matching of phases between the two set of data in both walls is obvious from the matching of peaks and valleys of the two signals.Therefore, it can be said that the space-averaged velocity and temperature fields are well correlated in the present flow configuration.The Nusselt number is observed to vary between 1.4 and 2.1 on the ribroughened wall and as expected it is more compared to the sidewall value.
Instantaneous moduli of vorticity and temperature distribution on cross-sectional plane at z 0.5 are shown in Figure 7.The vortex field near the leading edge of rib, indicated by circle, is considerably strong as compared to that on trailing edge.This leads to destabilization of flow field near the leading edge resulting in shear layer separation.Hightemperature gradient near the leading edge is also observed in this region along with a cold   fluid region in the bulk stream.The reason for existence of the high-temperature gradient is anticipated due to formation of additional shear layer near the tip.The hot fluid is dragged along by these separating shear layers into the core thereby increasing the heat transfer.Iso-surface of Q-criteria coloured by ω z , top and front views, is portrayed in Figure 8. Q-criteria, which deal with the second invariant of velocity gradient tensor, are defined by

Time-Averaged and rms Field
Figure 9 shows the time-averaged streamline at the duct center plane z 0.5 .The major characteristics of flow over a rib, namely, a large recirculation region behind the rib is clearly brought out.Additionally, two small secondary eddies, one between the primary recirculation zone and the rib, while the other in front of rib, are clearly resolved by the shearimproved model used in the study.The estimation of reattachment length, as reported in literature, is highly sensitive to the accuracy of predicting the level of turbulence in the shear layer.Overprediction of the reattachment length indicates that turbulence is underpredicted by the model and vice versa.The reattachment length is observed at around 5e which is in agreement with the results of Murata and Mochizuki 17 thereby indicating that the SISM model captures the turbulence level to a good accuracy near the wall.Figure 10 presents the time-averaged streamlines superimposed on temperature contours at three different streamwise locations.The overall secondary flow pattern at three locations look to be similar with only difference at the location of the ribs where the secondary vortices near the two sidewalls are found be weaker.This may be due to the fact that the highly accelerated flow at this location thus weakening the secondary flow structures.Figures 11 a -11 d show the distribution of resolved turbulent quantities at the midtransverse plane.Close examination of the contours of all quantities, as expected, reveals the highest normal as well as Reynolds shear stress within the separating shear layers.The corresponding line plot at x/e 1.5, just downstream of reattachment point, is shown in    The spanwise fluctuations w rms become the maximum in the stagnation region in front of the rib having a magnitude of 34.This high magnitude may be due to the flow impingement in front of the rib.The distribution of turbulent shear stress u v shows magnitude of −2.3 in the separated shear layer of the rib.The wall normal gradient of the three normal stresses is such that the streamwise u rms and spanwise fluctuations w rms show almost identical slope while the distribution of transverse fluctuations v rms reveals a lower value of the gradient.However, Reynolds shear stress −u v presented in Figure 12 shows a wall normal variation having the lowest slope.The time-averaged streamwise velocity and temperature contours at three different streamwise locations are compared and illustrated in Figures 13 and 14.The velocity contours are plotted for u * u × u avg /u τ and compliment well with the results of Murata and Mochizuki 17 .For the case of stationary ribbed duct, the similarity between the streamwise velocity and the temperature contours confirms the correlation between the streamwise velocity field and the corresponding thermal field.As discussed by Murata and Mochizuki 17 , the velocity boundary layer on the ribbed surface is found to be much thicker than temperature boundary layer at location, between two consecutive ribs.At the rib location, the two profiles remain more or less similar to each other.
Precise calculation of flow physics is one of important requirements for predicting the heat transfer accurately.Figure 15 shows the Nusselt number Nu/Nu 0 distribution on one of the ribbed walls and one of the sidewalls.As expected, the highest Nusselt number on rib wall occurs on the front face of the rib due to unsteady impingement of the shear layers.These eddies enhance the mixing of fluid near wall bringing cool fluid from core to wall and remove hot fluid from the wall, resulting in higher heat transfer.The augmentation ratio in comparison to a smooth duct in this region is 3.5.Immediately downstream of the rib, in the trapped recirculation zone, the heat transfer augmentation is very low.Further downstream of the rib, the heat transfer increases due to the presence of primary recirculation zone which improves fluid mixing.The maximum heat transfer occurs just near the reattachment region where the Nusselt number ratio is found to be 3.35.The heat transfer decreases as we move toward the smooth sidewalls.This is because of the fact that the formation of thick boundary layer on the either side bottom wall and the side wall.At the junction of rib with the smooth wall, high heat transfer is observed because of the presence of both the strong primary as well as secondary flow.The maximum Nusselt number ratio on top surface of rib not shown here is found to be 5.0 which decreases to unity as sidewall is approached.The maximum Nusselt number 8.0 occurs on the leading/front edge of the rib on which the flow impinges.The higher heat transfer Journal of Applied Mathematics   16 reveals the predicted augmentation ratio at the center plane of the ribbed wall, z 0.5.Two peaks, one in front of the rib and other near the reattachment length, as described earlier, can be clearly identified from this figure.Table 1 shows the comparison of wall-averaged Nusselt number with the data of Wagner et al. 4 .The friction factor value for the geometric configuration studied is evaluated equal to 0.26.

Conclusion
Large eddy simulation of flow through stationary ribbed duct with normal ribs has been performed using the shear improved Smagorinsky SISM model.The instantaneous flow and temperature fields reveal the complex nature of the turbulent flow resulting due to the presence of periodic rib turbulators.It is also found that the time-averaged flow and heat transfer data obtained is concurrent with the results obtained by dynamic LES as well as experimental data.The instantaneous flow and temperature field show all the nuances of the wall-bounded separated shear flow.Secondary flow and the turbulent fluctuations are also captured with good accuracy which results in precise prediction of surface heat transfer.This shows the shear improved Smagorinsky SISM can be used even for nonstationary flow as an alternative to the dynamic subgrid stress model, which is found to be computationally more expensive.Prandtl number, ν/α q j : Turbulent thermal subgrid stress Re:

Nomenclature
Reynolds number, u avg D h /ν t: T ime f: Friction factor T : Temperature u τ : Frictional velocity, τ w /ρ y : yu τ /ν u, v, w: Nondimensional velocity in x, y and z directions also u i and u j u i , u j : Velocity fluctuations x, y, z: Nondimensional Cartesian coordinates also x i and u j .
Lévêque et al. 25 modified the Smagorinsky model as ν T C S Δ 2 |S| − | S | where the magnitude of the mean shear | S | is subtracted from the magnitude of instantaneous resolved rate of strain tensor |S|.

3 a 3 b 5 Figure 2 :
Figure 2: Time variation and power spectra of streamwise velocity at a downstream of the rib; and b center of the computational domain.

Figure 3 :
Figure 3: Comparison of streamwise and cross-stream velocity with the reference data all results are at Re 12, 000 .

Figure 4 :
Figure 4: Instantaneous velocity vector superimposed on temperature contours on cross-planes a x 0.5 and b x 1.0.

Figure 5 :
Figure 5: Instantaneous contours of a normalized Nusselt number Nu/Nu 0 and b skin friction coefficient scaled with Reynolds number C f × Re on one of the sidewalls and bottom ribbed walls.

Figure 6 :
Figure 6: Time variation of normalized Nusselt number Nu/Nu 0 and skin friction coefficient scaled with Reynolds number C f × Re on one of the a sidewalls b bottom ribbed walls at an arbitrary time instant.

2 Figure 7 :Figure 8 :
Figure 7: Instantaneous moduli of vorticity and temperature distribution on cross-sectional plane at z 0.5.a modulus of vorticity, and b temperature field.

Figure 10 :
Figure 10: streamline superimposed on temperature contours at three streamwise locations a center of rib x/e 5 , b at x/e 7.5 , and c downstream of reattachment point x/e 1.5 .

Figure 11 :Figure 12 :
Figure 11: Distribution of resolved turbulent quantities a u rms ×100 , b v rms ×100 , c w rms ×100 , and d u v ×100 at center plane z 0.5 .

Figure 13 :
Figure 13: Time-averaged contours of streamwise velocity, u * u×u avg /u τ locations same as in Figure 10 .

Figure 12 .Figure 14 :
Figure 12.It is to be noted that the values of different stresses are scaled up by 100.The maximum value of u rms is 41.2 which occurs in the separated shear layer downstream of the rib and is the lowest 3.3 in the secondary recirculation region trapped between the rib and the primary recirculation region formed behind the rib.On the rib-roughened floor the magnitude remains around 19.The boundary layer Figure12is still recovering at this location.The boundary layer profile, however, is not revealed by the cross-stream fluctuation v rms ; instead, it monotonically reaches maximum value of 21 in the separated shear layer.At the front stagnation region, the magnitude of v rms is 14, whereas in the separated shear layer, the value ranges to 16.

Figure 15 :
Figure 15: Time-averaged Nusselt number on one of the a ribbed walls b sidewalls.

Figure 16 :
Figure 16: Time-averaged Nusselt number variation at the central transverse mid-plane ribbed wall z 0.5 .

Table 1 :
Comparison of averaged Nusselt number with Wagner et al. 4 .
augmentation on rib-roughened surface than that on the smooth wall is quite obvious from Figures 15 a -15 b .Figure