Numerical Simulation of Turbulent Fluid Flow and Heat Transfer in a Ribbed Rotating Two-Pass Square Duct

The local turbulent fluid flow and heat transfer in a rotating two-pass square duct with 19 pairs of in-line 90◦ ribs have been investigated computationally. A Reynolds-averaged Navier-Stokes equation (RANS) with a two-layer k − ε turbulence model was solved. The in-line 90◦ ribs were arranged on the leading and trailing walls with rib height-to-hydraulic diameter ratio and pitchto-height ratio of 0.136 and 10, respectively. The Reynolds number, based on duct hydraulic diameter and bulk mean velocity, was fixed at 1.0 × 104 whereas the rotational number varied from 0 to 0.2. Results are validated with previous measured velocity field and heat transfer coefficient distributions. The validation shows that the effect of rotation on the passage-averaged Nusselt number ratio can be predicted reasonably well; nevertheless, the transverse mean velocity and, in turn, the distribution of regionalaveraged Nusselt number ratio are markedly underpredicted in the regions toward which the Coriolis force is directed. Further CFD studies are needed.


INTRODUCTION
Gas turbine blades are cooled internally and externally to raise turbine rotor inlet temperature for better thermal efficiency and power output of gas turbines.This paper concerns the fluid flow and heat transfer in turbine blade internal coolant passages.Among the relevant studies reported, computations or measurements of both turbulent fluid flow and heat transfer in a rotating two-pass square-sectioned duct with 90 • ribs are sparse.In particular, accurate predictions of the complex heat and fluid flow phenomena resulting from Coriolis force, centrifugal buoyancy, sharp turn, and shear layer over ribs have not yet been achieved.In this study, a preliminary attempt is taken to computationally simulate this challenging problem.
This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.
Data obtained from rotating smooth duct are often used as a base for comparison with those obtained from rotating ribbed duct.For rotating smooth ducts, Liou and Chen [13] presented a set of laser-doppler velocimetry (LDV) data for developing flows through a 180 • straight-corner turn at a rotating number Ro = 0.082 to investigate the effect of 180 • turn under rotation.Servouze [12] performed threedimensional LDV measurements in a rotating two-pass channel with a 180 • sharp turn.His results were mainly obtained at Re = 5000 and Ro = 0.33.Some results were obtained at Re = 25 000 and Ro = 0.033 and 0.066.Chen and Liou [8] investigated rotating effect on fluid flow in a smooth duct with a 180 • sharp turn for Re = 1 × 10 4 and Ro from 0 to 0.2.They found that both the skewness of streamwise mean velocity and magnitude of secondary flow velocity increased linearly, and turbulence intensity increased nonlinearly with increasing Ro.Their results also illustrated the reported spanwise heat transfer variation in the rotating channels, especially the high heat transfer enhancement on the leading wall in the turn.More importantly, they concluded that the direction and strength of the secondary flow with respect to the heat transfer wall were the most important fluid dynamic factors affecting the local heat transfer distributions inside a 180 • sharp turn under both rotating and nonrotating conditions.
In practical applications, ribs are often arranged on the channel wall to enhance the heat transfer.Tse and Steuber [9] investigated flow characteristics in a serpentine coolant passage with 45 • ribs using LDV.Liou et al. [1] used both LDV and TLCT (transient liquid crystal thermometry) to measure the local fluid flow and surface heat transfer distributions in an orthogonally rotating two-pass square-sectioned duct with in-line 90 • ribs mounted on the leading and trailing walls.Up to date, detailed flow field measurements using nonintrusive LDV or PIV in a rotating multipass ribbed channel with 180 • sharp turns are still lacking.For the studies of heat transfer only, Wagner et al. [10] made thermocouple measurements in serpentine square passages with staggered 90 • ribs on the leading and trailing surfaces.The ribs' H/D H and Pi /H were fixed at 0.1 and 10, respectively.Flow parameters were Re = 2.5 × 10 3 and Ro varying from 0 to 0.35.They showed that the maximum heat transfer coefficient increased up to 4.5 times (in contrast to 3.5 times for the rotating smooth-walled passages) the stationary fully developed circular tube values, whereas the minimum heat transfer coefficients decreased to 80% of the stationary 90 • -ribbed wall results.When the above 90 • ribs were replaced by 45 • ribs, Johnson et al. [2] found that the maximum heat transfer coefficient increased up to 5 times the stationary tube values but the minimum heat transfer coefficient decreased to 40% of the stationary 45 • -ribbed wall values.Parsons et al. [11] performed thermocouple measurements to investigate wall heating effect on heat transfer in a rotating two-pass square channel with in-line 90 • ribs for Re and Ro from 2.5 × 10 3 to 2.5 × 10 4 and 0 to 0.35, respectively.Zhang et al. [14] substituted in-line 60 • ribs for the in-line 90 • ribs of Parsons et al. [11] and examined the effect of rotation on heat transfer using thermocouples.Azad et al. [15] studied the effect of the channel orientation (β = 90 • and 135 • from the direction of rotation) on two-pass rectangular channels for both smooth and 45 • -ribbed walls.For the ribbed wall case, the channel orientation effect was less sensitive compared to the smooth-wall case.In all cases the rotating ribbed wall heat transfer coefficients were two to three times higher than their corresponding rotating smooth-walled channel values.
As to the numerical studies, a low-Re k − ω turbulence model had been implemented by Stephens and Shih [5] to investigate the effect of angled ribs on the heat transfer coefficients in a rotating two-pass duct.The main parameters examined were Reynolds number, rotation number, and buoyancy parameter.Their results indicated that the regionalaveraged heat transfer coefficient measured by Johnson et al. [2] was underpredicted for the stationary case.However, no conclusive comparison was made for the rotating case.Bonhoff et al. [6] computed the flow fields and heat transfer coefficients for rotating U-shaped coolant channels with 45 •angled ribs.A Reynolds stress turbulence model with wall functions for the FLUENT CFD code was used.Their results underpredicted the measured regional-averaged Nusselt number ratios of the first-pass leading and trailing surfaces.Rigby [7] predicted fluid flow and heat transfer in a 90 • -ribbed rotating coolant passage with a 180 • turn using a k−ω turbulence model.The results showed that the regionalaveraged heat transfer coefficient was overpredicted for the stationary case and underpredicted for the rotating case.Chen et al. [3,4] used two turbulence models: a two-layer k−ε model and a near-wall second-moment closure model to study heat transfer in rotating two-pass square channel with smooth walls.The near-wall second-moment closure model produced better regional-averaged heat transfer predictions in comparison with the data of Wagner et al. [10] whereas the two-layer k − ε model significantly underpredicted the experimental data.Using the same second-moment model, Al-Qahtani et al. [16] calculated flow and heat transfer in a rotating two-pass square channel with 45 • -angled ribs and channel aspect ratio of 2 : 1.The focus was to investigate the effect of channel aspect ratio and channel orientation on the flow and heat transfer.The results compared reasonably with regional-averaged heat transfer data of Azad et al. [15] for both stationary and rotating ribbed channels.
One of the reasons why extensive numerical studies of the fluid flow and heat transfer in rotating ribbed ducts with 180 • sharp turns have not been possible to date is the lack of local flow field and local (not regional-averaged) surface heat transfer data.This detailed local information can be beneficial to the validation and improvement of the numerical simulations.In addition, the boundary conditions including inlet and outlet flow condition and thermal condition of wall are crucial to the simulation.Most numerical works mentioned above fail in the validation of the fluid flow results and adoption of appropriate boundary conditions.In the present study, the data base of Liou et al. [1], including local fluid flow and surface heat transfer distributions, is used to validate the computed fluid flow and heat transfer in an orthogonally rotating two-pass square-sectioned duct with in-line 90 • ribs mounted on the leading and trailing walls.The thermal boundary condition adopted in previous numerical simulations was either constant wall temperate or constant heat flux.In order to implement a more realistic thermal boundary condition, the wall temperature distributions obtained from Liou et al. [1] are adopted in the present study (Note that a list of nomenclature presented in the paper are set in Table 1).In addition to local fluid flow validation, the investigation is also focused on the effect of rotation number on the local, regional averaged, and passage-averaged Nusselt number.

INTERNAL COOLANT CHANNEL SIMULATED
A sketch of the configuration and dimensions of the internal coolant channel simulated is depicted in Figure 1.The hydraulic diameter of the square cross-sectional flow path was D H = 22 mm and the divider-wall thickness was 0.5 D H .At the turn, the clearance between the tip of the divider wall and the duct outer wall was fixed at 1 D H .The 90 • transverse ribs were mounted on the leading and trailing walls, directly opposite one another (not staggered).One pair of ribs was installed at the divider-wall tip in the turn and 9 pairs of ribs in each pass.The rib height-tohydraulic diameter ratio (H/D H ) and the rib pitch-to-height ratio (Pi /H) were 0.136 and 10, respectively, in each pass.
The Reynolds number, based on the bulk mean velocity of 7.58 m/s and hydraulic diameter, was fixed at 1.0 × 10 4 .The range of the rotation number examined was from 0 to 0.20 corresponding to the range of rotational speed from 0 to 600 rpm.These conditions were selected to compare with the experimental data reported by Liou et al. [1].

Governing equations
Figure 1 also depicts the coordinate system to be used for formulating the governing equations and describing the results.The viscosity and thermal conductivity are assumed constant and the viscous dissipation and compression work in the energy equation are neglected.To capture the various forces in a rotating system, the centrifugal-buoyancy is not neglected.
It is worth noting that centrifugal-buoyancy may not be neglected as the rotational rate is large.The time-dependent Navier-Stokes equations and energy equation in a rotating frame of reference with the axis of rotation aligned with the spanwise direction take the form  The last two terms in (2) are the Coriolis forces and centrifugal-buoyancy forces experienced by a fluid element resulting from system rotation.L x and L y are the distances from the center of rotation in the x and y direction, respectively.

Turbulence models
The Reynolds stress −u i u j and heat flux −u i T shown in ( 5) and ( 6) must be modeled to close the governing equations ( 4), ( 5), and (6).A two-layer k − ε turbulence model [17] was used.In the two-layer approach the standard k − ε model is used for the core region, whereas the one-equation model is used for the near-wall region which includes the sublayer, the buffer layer, and a part of the fully turbulent layer.Unlike the low-Reynolds-number model that requires the solution of the transport equations of both k and ε all the way to the wall, the one-equation model only requires the solution of the turbulent kinetic energy equation in the near-wall region.
The rate of energy dissipation and eddy viscosity in the near-wall region are obtained as The turbulent length scales that account for the necessary damping effects in the near-wall region in terms of the turbulence Reynolds number R y are adopted.The length scales proposed by Wolfshtein [18], modified by Chen and Patel [17], can be expressed as where y is the normal distance from the wall and According to the suggestion of Chen and Patel [17], C 1 is given by where κ is the von Karman constant, to ensure a smooth eddy-viscosity distribution at the junction of inner and outer layer.A ε = 2C 1 is assigned so as to recover the proper asymptotic behavior where ν molecular kinematic viscosity, in the sublayer.The third parameter, A µ = 70, is determined from numerical tests to recover the additive constant in the logarithmic law in the case of a flat plate boundary layer.

Boundary conditions
The measured mean velocity and turbulence intensity profiles for various values of Ro at X/H = 73 (or X * = 10) station [1] located in between the 8th and 9th rib pairs (Figure 1) upstream of the turn in the first pass were used as the inlet flow conditions.This station corresponds to the inlet reference station reported by Liou and Chen [13] in their study of rotating effect on fluid flow in a smooth duct with a 180 • sharp turn.The dissipation rate of turbulent kinetic energy was calculated from the prescribed inlet turbulent kinetic energy and a length scale (0.01 times channel height).At the outlet, a convective boundary condition was prescribed.No slip condition and measured wall temperature distribution were applied to the walls.

NUMERICAL METHOD OF SOLUTION
By using the control-volume-based finite difference method of Patankar [19], governing equations were integrated over a finite number of control volumes covering the entire solution domain with each control volume containing a grid point.
The convection terms of all equations, including the equations of scalar variables, were discretized by a third-order upwind scheme QUICK [20] and the diffusion and source terms were discretized by the second-order central difference.The time derivative was discretized by a first-order accurate fully implicit method with a small time step for timeaccurate solutions and a larger time step for a steady solution.
As for handling the pressure-velocity coupling arising in incompressible flow equations, the PISO (pressure implicit with splitting of operators) [21] method was utilized in this study for a stable and faster convergence at each time step.The grid system was nonstaggered and the momentum interpolation method of Rhie and Chow [22] was adopted.The resulting set of algebraic equations was solved with the lineby-line TDMA (tridiagonal-matrix algorithm).The calculation results at each time step were declared convergent when the maximum of individual sum of absolute cell residues in the entire domain for each dependent variable normalized by the respective inlet flux was less than a prescribed convergence criterion, typically a value of 0.05%.
The computations to be reported in the following sections were performed on a 312×48×84 grid.Additional runs for the coarser (212 × 32 × 56) and finer (422 × 64 × 112) meshes were undertaken for a check of grid independence.The distance from the nearest nodes to the wall was kept constant (y + 0.1) for the three runs.A comparison of the results showed that the maximum changes of 1.2% in the axial velocity profiles and 1.8% in Nusselt number distribution between 312 × 48 × 84 and 422 × 64 × 112 grid sizes were smaller than 3.5% and 4.7% found between 312 × 48 × 84 and 212 × 32 × 56 mesh sizes.Consequently, the accuracy of the solutions on a 312 × 48 × 84 grid size is deemed satisfactory.

Flow development around the turn in Y * planes
The mean velocity vector distributions in the central Y * = 0 and near-wall Y * = −0.9planes parallel to the leading and the trailing surfaces are shown in Figure 2 for Ro = 0.15.For clearness, only part of the velocity vectors is shown.Fluid flow in these planes is perpendicular to the direction of the Coriolis force and, therefore, more directly affected by the presence of ribs and 180 • sharp turn.In the first pass (X/H > 0 and −1 < Z * * < 0), the fluid flows from right to left.The mean flow basically parallels the divider wall and no reverse flow can be found in the Y * = 0 plane which is 2.7 rib heights away from the rib tops of each rib pair.The flow just gradually adjusts itself for making a 180 • turn.In the Y * = −0.9plane, about midway between the rib top and leading surface, the flow exhibits a divergent nature and is rather symmetric to the centerline.An examination of the first two rows of velocity vectors immediately behind each rib indicates the presence of rib-top separated recirculating flow zones and flow reattachments.Note that due to the presence of boundary layers on Z * * = −1 and Z * * = 0, the reattachment lengths adjacent to are longer than those away from the side wall (Z * * = −1) and the divider wall (Z * * = 0), resulting in a curved mean reattachment line.The envelope of the curved mean reattachment line, rib rear edge, side wall, and divider forms a concave-planar lens-like shape.One may expect a similar shape of low heat transfer coefficient contour, as shown in Figure 3, since the rib-top separated recirculating flow zone often deteriorates heat transfer [1].Also note that the above-mentioned divergent flow pattern leads to two small corner vortices immediately in front of each rib, as indicated by the enlarged view b in Figure 2. Thus, two small heat transfer deterioration zones are also expected to exist to reflect the presence of these two corner vortices, as evidenced from Figure 3. Inside the turn, the flow pattern in the central plane Z * * = 0 is similar to that which occurred in a rotating smooth-wall duct with a corner vortex around the junction of blade tip and Z * * = −1 (enlarged view a in Figure 2).The combination of the 180 • turn-induced secondary flow, rib's disturbance, and duct rotation, however, results in a very complex flow pattern in the plane Y * = −0.9.In particular, the Coriolis force in the rear half of the turn points toward the leading wall (Y * = −1) making the disappearance of rib-top separated separation bubble behind the rib located at Z * * = 0.The larger mean velocity vectors in this part of the turn provide the rationale for the corresponding heat transfer enhancement shown in Figure 3. Attention is now focused on the fluid flow immediately downstream of the 180 • sharp turn and near the divider wall.The well-known turn-induced separation bubble normally taking place near the divider tip in the turn and reattaching on the divider wall in the second pass for the smooth duct case is practically absent in the Y * = 0 plane of the 90 • -ribbed duct, but still can be found between 0 < X/H < 10 in the Y * = −0.9plane, a result in agreement with that reported from the experiment of Liou et al. [1].The high heat transfer enhancement on the leading wall of the second pass (Figure 3) is attributed to the direction of the Coriolis force being directed toward the leading wall.Due to the centrifugal force associated with the 180 • turn, the axial mean velocity profiles in the Y * = 0 plane have their peaks toward the outer surface just after the turn.As the flow proceeds toward the outlet, the velocity vector distributions gradually recover their symmetry as in the first pass.

Flow development around the turn in Z * planes
As the flow proceeds toward the turn in the first pass, the flow patterns in terms of the mean velocity vector plots in the Z * = −0.5 plane predicted in the present study and measured by Liou et al. [1] are depicted in Figure 4 for Ro = 0.15.Both prediction and measurement show that in the first pass, the velocity vectors are directed toward the trailing wall by the Coriolis force.With the maximum streamwise mean velocity at each X/H skewed toward the trailing wall.However, the computation underpredicts the Coriolis effect resulting in a smaller transverse mean velocity component near the trailing wall.Inside the 180 • sharp turn, the combined effects of curvature and rotation make the flow pattern in the region −1 < X * < 0 (or −7.3 < X/H < 0) of Z * = −0.5 plane assume a skewed Dean-type vortex pair.Computational results are able to reveal the presence of the squeezed Dean vortex located at left top corner (−1 < X * < −0.5, 0.5 < Y * < 1) adjacent to the blade tip and trailing wall.Note that due to technical difficulty, LDV measurements could not be performed in the corresponding region.
In the Z * * = 0 cross-sectional plane, Figure 5 depicts that the dominant vortex occupies more than 80 percent area of the mid-turn cross-section and directs the major secondary flow to impinge and sweep the leading wall.Some differences between the prediction and experiment can be noticed for the flow patterns in the right top half of Figure 4.In the rear part of the turn, the Dean-type flow patterns (Figure 6) are no longer recognized because the flow has turned about 180 • and is ready for flowing radially inward.The direction of Coriolis force is reversed.More fluid flow is now gathered near the leading wall in the region −1 < X * < 0 of the Z * = 0.5 plane (Figure 6), opposite to that in the region −1 < X * < 0 of the Z * = −0.5 plane (Figure 4).Further downstream, the mean flow profiles are skewed toward the leading wall, as shown by the predicted and measured results in Figure 6.

Regional-averaged Nu distribution
Figure 7 depicts the computed and measured [1] regional averaged Nusselt number ratio (Nu rg / Nu 0 ) distributions in the ribbed rotating two-pass duct with a sharp 180 • turn for Ro = 0.15 and Re = 10 000.Both computations and experiments give the same trend of Nu rg / Nu 0 distribution, although computations underpredict quantitatively.The Nu rg / Nu 0 on the leading wall in the first pass (region indices −4 to −2) is lower than that on the trailing wall (Figure 7a) since the Coriolis force is directed toward the trailing wall for the radially outward duct flows.The Nu rg / Nu 0 on the leading wall in the second pass (region indices 2 to 4), however, is higher than that on the trailing wall since the Coriolis force is now also directed toward the leading wall for the radially inward duct flows.Inside the 180 • turn, Nu rg / Nu 0 is the highest as a result of enhanced convective and impinging heat transfer by the flow curvature, turn-induced secondary flow, and presence of ribs.
It should be pointed out that almost all the previous CFD studies underpredicted Nu rg / Nu 0 in rotating smooth [3,4] and ribbed (90 • or 45 • ) [5,6,8] ducts, as has been mentioned in the introduction.One of the reasons for the Leading, experiment [1] Trailing, experiment [1] Leading, number Trailing, number  underprediction of Nu rg / Nu 0 is apparent from the computed and measured flow pattern comparisons given in Figures 4, 5, and 6.As addressed before, computations tend to underpredict the transverse mean velocities induced by the Coriolis force, particularly in the turn (Figure 5), and near the trailing wall of the first pass (Figure 4) and the leading wall of the second pass (Figure 6).As a result, Nu rg / Nu 0 in the corresponding regions shown in Figures 7b and 7c are significantly underpredicted.Figure 7d cites the heat transfer in a stationary 90 • two-pass ribbed square channel reported by Johnson et al. [2] for Re = 3.0 × 10 4 .It is seen that even for the stationary case, the two-layer isotropic eddy-viscosity turbulence model fails to capture the steep increase of Nu rg / Nu 0 .Although the effect of the Reynolds stress anisotropy is accounted for in the second-moment closure model and the comparison shown in Figure 7d is improved, some large differences at several streamwise stations still can be found.Consequently, more extensive assessment of various turbulence models is recommended for the future study.Large eddy simulations are also encouraged.
It is also interesting to know the effect of rotating number on the passage-averaged heat transfer (Nu p ) rr for a ribbed rotating two-pass duct.Figure 8 depicts such a plot where (Nu p ) rr is normalized by its stationary counterpart (Nu p ) rs .As one can see from Figure 8, approximately linear function dependences of (Nu p ) rr /(Nu p ) rs on Ro are found for all passages on both the leading and trailing walls.More specifically, for the Ro range examined, (Nu p ) rr /(Nu p ) rs linearly increases and decreases with increasing Ro on the trailing and leading walls of the first pass whereas the trend is reversed in the second pass.In general, the present simulations and previous experiments show similar (Nu p ) rr /(Nu p ) rs versus Ro trends on both leading and trailing walls of the first and second passes.Quantitative agreement between the computed and measured (Nu p ) rr /(Nu p ) rs is also good.However, caution should be brought to one's attention that the comparison made for Nu rg / Nu 0 versus region index in Figure 7 is rather poor.

CONCLUSIONS
Turbulent fluid flow and heat transfer simulations have been performed in a rotating two-pass square-sectioned duct with 90 • in-line ribs mounted on the leading and trailing walls.The previous LDV and TLCT measured data are used for imposing velocity and thermal boundary conditions, and validating the computed results.The predicted flow patterns successfully explain the measured local heat transfer distributions.As the computed local Nusselt number distributions averaged over the passage (Nu p ) rr and normalized by its stationary counterpart (Nu p ) rs are plotted versus rotating number Ro examined, reasonably good agreement, both qualitatively and quantitatively, is attained.The (Nu p ) rr /(Nu p ) rs on the leading and trailing walls of the first and second passes can be correlated as linear functions of Ro.Although the streamwise evolution of the regional-averaged Nusselt number ratios for Nu rg / Nu 0 along the flow pass is captured qualitatively, the values are significantly underpredicted inside the 180 • sharp turn, on the trailing wall of the first pass, and on the leading wall of the second pass.Similar underpredictions were also found in the previous CFD studies using even Reynolds stress turbulence models.Two causes can be identified for the above underpredictions.First, the two-layer turbulence model adopted is known for the failure in accounting for the effect of turbulence anisotropy, particularly inside the turn.Second, the computed flow field shows the underprediction of the transverse mean velocity directed by the Coriolis force around these regions.The reason for the second cause needs to be further investigated in a future study using, say, large eddy simulations.

ACKNOWLEDGMENT
Support for this work was partially provided by the National Science Council of R.O.C under contract NSC 89-2212-E-007-135.

Figure 1 :
Figure 1: Sketch of configuration, coordinate system, and dimensions of test section.

Figure 7 :
Figure 7: (a) Region index.Computed regional-averaged Nusselt number ratios compared with the measured data [1] for Ro = 0.15 and Re = 10 000 on (b) leading and (c) trailing walls.(d) Presently computed and previously computed and measured [23] local Nusselt number ratios.

Figure 8 :
Figure 8: Comparison between presently computed and previously measured passage-averaged Nusselt number ratio versus rotation number in (a) first pass and (b) second pass.