Modeling of Laminar-Turbulent Transition in Boundary Layers and Rough Turbine Blades

Laminar-turbulent transition occurs in a broad range of industrial applications, e.g., turbine or rotorcraft rotor blades, which sometimes operate in harsh conditions, leading to surface roughness. Roughness could have great effect on transition. It is important to be able to estimate the influence of transition and roughness on their aerodynamic and/or heat transfer performance, through computational fluid dynamics (CFD) modeling.

Accurate prediction of laminar-turbulent transition remains a difficulty for general-purpose CFD engineering modeling. One reason is that several different mechanisms lead to transition in different situations. For aerodynamic flows with very low freestream turbulence intensities (⁠$≲0.2%$⁠), laminar to turbulent transition starts with linear instabilities, such as Tollmien–Schlichting waves, which develop downstream as nonlinear modes, and final break down to turbulence. This is called natural transition [1,2]. However, for freestream turbulence intensities $≳1%$⁠, the transition occurs via diffusion of turbulence into the laminar boundary layer, i.e., bypass transition [3–7]. There are also other transition mechanisms, such as separation-induced transition that may occur in the detached shear layer. It is quite difficult for one model to accommodate all such mechanisms.

Many transition models have been proposed over time. One focus has been to combine Reynolds-averaged Navier–Stokes (RANS) based turbulence models with a transition criterion. One criterion is the e^N method predicting transition based on boundary layer stability theory [8]. Another simpler criterion for transition prediction is to use a data correlation for critical Reynolds numbers. For instance, an empirical critical Reynolds number is compared to the flow Reynolds number and used to switch on the RANS model [9–12]. This approach requires boundary layer integral properties to be computed. To avoid boundary layer integrals, Langtry and Menter [2] proposed a method that expresses a correlation for the critical vorticity Reynolds number [13] in terms only of local boundary layer quantities. Also, it couples two additional transport equations (for the intermittency function γ and the transition vorticity Reynolds number Re_θt) with the Menter shear stress transport k–ω based RANS turbulence model. It is often called a “correlation-based” γ–Re_θt transition model. Most recently, this model was extended to include crossflow transitional effects [14].

A simpler approach is to use only one transport equation for the intermittency function γ and only local variables. Such a model was first proposed in Ref. [15], using the intermittency approach [16]. The model was further developed in Refs. [17] and [18]. Some further improvements are described herein.

Section 1.1 cited studies on transition modeling on smooth surfaces; there are only limited studies about transition modeling on rough surfaces. Roughness effects in wall-bounded flows have been referred to as an “Achilles Heel of CFD” [19,20]. George et al. [20] used Durbin's v² − f model [21] to study roughness effects on a turbulent boundary layer. Although there have been many experimental studies of turbulent flows over rough walls (see Ref. [22] for a review), there are very few benchmark datasets on transition on rough surfaces. Here, we use the turbine data of Stripf et al. [23].

With the development of high performance computing, direct numerical simulation (DNS) and large eddy simulation have recently been used on rough wall-bounded flows [24–27]. Since DNS resolves all scales and large eddy simulation only models small scales, both methods can give accurate predictions, but both methods require fine grid and high computing power.

Many researchers [28–33] used the discrete element method where a form drag term was added to the momentum equation, and the blockage effects of rough elements on the flow were considered. This method produced good results, but it was more suitable for known regular and periodic roughness geometries, and it needs detailed information about the roughness characteristics [34].

Roughness is commonly characterized by an equivalent sandgrain roughness height (r), as first proposed by Schlichting [35]. He defined r as the size of sandgrain in a certain pipe flow experiment producing the observed skin friction. Some researchers [36–42] proposed correlations to compute r, where the correlation [38] only uses statistical parameters of the surface roughness to calculate r.

A displacement of origin concept [43] was introduced through prediction of the shift of velocity profiles due to roughness, which provides a way to incorporate the equivalent sandgrain roughness into RANS turbulence models, through modified boundary conditions. A hydrodynamic roughness is obtained from the sandgrain roughness via a calibration procedure that fits the log-law displacement predicted by the model to an empirical formula. Then, the rough surface is represented by modified boundary conditions for model field variables. The method is applied to k–ω based RANS turbulence in Refs. [44] and [45]. Another objective of the current study is to apply the displacement of origin and equivalent sandgrain roughness concepts to bypass transition modeling.

The rest of the paper is organized as follows: Section 2 is the introduction and formulation of the model for smooth walls, followed by roughness modification to the model to account for roughness effects on transition. Section 3 is mainly about validation tests on different type of applications: flat plate with zero, favorable, and adverse pressure gradients, periodic wakes impinging on a flat plate, rough low, and high pressure turbine blade (HPT). The final section is a summary of this work.

The objective of the development of the current model is to predict transitional flows with freestream turbulent intensities Tu at low (Tu < 1%), moderate(1% < Tu < 3%), and high Tu > 3% levels. It is a challenging task. Previous studies [46,47] have shown that transition onset momentum thickness Reynolds number Re_θt increases slightly, or remains nearly constant, i.e., in a region around Re_θt ∼ 200, with the decrease of Tu from high to moderate levels, however, Re_θt increases significantly to above 1000 when Tu decreases from moderate to low levels.

The relation between Re_θt and Tu was considered in the development of the current intermittency γ transition model, although the current model does not use Re_θt and Tu. Instead it uses similar local variables in order to make the model “local” for general CFD purpose by removing the integral variables, and simplify the model formulation by removing the transport equation for Re_θt.

where σ_l = 5.0 and σ_γ = 0.2 are empirical coefficients for laminar (kinematic) and turbulence (eddy) viscosities, respectively, which controls the diffusion of γ. σ_γ can affect the transition length, that is, the distance from the transition onset location to the full turbulence location. For instance, an increase of the eddy viscosity parameter σ_γ can reduce the diffusion of turbulence, delay the location of full turbulence, and increase transition length. This feature is used in the model formulation for roughness, presented in Sec. 2. On the other hand, decrease of σ_l, leads to an increase of the effect of laminar viscosity on diffusion, and thus delays transition onset. It seems to have a great effect on the laminar end of transition. Some related sensitivity analysis of σ_γ and σ_l can be found in Ge et al. [17].

Note that R_t is the ratio between eddy viscosity and kinematic viscosity, and Ω is the mean rotation rate. T_ω is related to turbulence levels, and R_c1 is the critical Reynolds number for transition to turbulence. In order to accommodate transition with low Tu, an inverse relation between R_c1 and T_ω is assumed, similar to the relation between Re_θt and Tu. R_ν is wall-distance-based, vorticity Reynolds number. It is similar to momentum thickness Reynolds number Re_θ in a Blasius boundary layer; that is, the maximum R_ν is equal to Re_θ [2].

The function F_γ turns on the production term when R_ν is between 1.08R_c1 and 90. It can be seen that the higher R_c1, the smaller the effective region of F_γ. γ_max > 1 is used to drive γ to 1. The first minimum in Eq. (11) ensures a γ ≤ 1.

where G_γ, similar to F_γ, is used to turn on/off the sink term. It becomes effective when 18 < R_ν < 100. The factor, F_turb, makes the sink negligible outside of the laminar boundary layer, where R_t is high.

Otherwise, this is the standard model, with $νT=k/ω$⁠, $Ck=0.09, Cω1=5/9, Cω2=3/40$⁠, $σω=σk=2$⁠.

where $|S|=SijSji, Sij=0.5((∂ui/∂xj)+(∂uj/∂xi))$⁠, and n_w is the unit wall normal vector.

Roughness affects both the onset and length of transition. Experimental studies [34,47,48] on a rough turbine blade have shown both earlier transition, and increased transition length in the “transitionally rough” turbulent regime.

In the current model, the rough surface is replaced by a smooth surface with modified boundary conditions to account for roughness effects. The model was modified by introducing dependence on equivalent sandgrain roughness r.

Note that σ_γ = 0.2 when r = 0, which is the values for smooth walls. R_r is similar to r⁺= ru_*/ν, but uses local variables.

The value C_r = 0.26 was selected to set how strong the influence of roughness is. C_rr is the effective displacement of the origin; some analysis of its influence on transition can be found in Ref. [18].

where Q is the difference between the rate of rotation and the rate of strain. Q > 0 for favorable pressure gradient, Q = 0 for zero pressure gradient, and Q < 0 for adverse pressure gradient. Hence, F_turb is reduced in regions of strong acceleration.

A favorable pressure gradient, i.e., a positive Q, leads to a positive R_c3, and thus, increases the upper bound of G_γ. Note that R_c2 was introduced to increase the lower bound of G_γ so that the roughness-induced R_ν will not spuriously turn on the sink term very near the wall, in the fully turbulence region.

The current model is validated on both smooth and rough wall cases in Sec. 3.

The new model was first tested on transitional flows over a smooth, flat plate, zero-pressure gradient boundary layer, with different freestream turbulence intensities (Tu_in) and freestream velocities (U_∞). With a plate length of 1.5 m and a kinematic viscosity of $ν=1.5×10−5m2s−1$ ERCOFTAC cases T3A-, T3A, and T3B [49] correspond to $U∞=25 m/s, 5.2 m/s,9.4 m/s$⁠, and $Tuin=0.9%, 3.5%, 6%$⁠, respectively.

The inlet turbulent kinetic energy $kin=1.5(TuinU∞)2$⁠, specific dissipation rate $ωin=kin/(Rtν)$ and R_t were determined by matching the computed Tu(x) to data. The freestream turbulence intensity profile is plotted in Fig. 1. It agrees excellently with the experimental data. We found R_t = ν_t/ν = 8.7, 14, 100 for T3A-, T3A, and T3B, respectively.

The predicted skin friction coefficient C_f for the current model is compared with the model of Ge et al. [17] and the experiment data for T3A-, T3A, and T3B in Figs. 2–4, respectively. It can be seen from Fig. 2, that for the low Tu case T3A-, the current model significantly improves the prediction.

It is might be noticed that the predicted C_f before transition for T3A- is higher than the experimental values. This is a spurious influence of $FRs$ on γ_eff (Eq. (11)). Non-zero R_s around the leading-edge of the plate seems to have an influence on skin friction for transitional flows with low freestream turbulence.

The current model improves the prediction for flows with very low freestream Tu, but at the same time, it maintains the good prediction for higher Tu. Figures 3 and 4 show that the current model gives similar predictions as the model of Ge et al. for cases T3A and T3B.

The new model was also tested on the T3C1, T3C2, T3C3, T3C4, and T3C5 cases; that is, transitional flows over a smooth, flat plate with favorable/adverse pressure gradients, with different inflow freestream turbulence intensities (Tu_in) and freestream velocities (U_∞). The computational domain for T3C cases has a bottom flat plate with a length of 1.65 m. The variable pressure gradients are created by a top curved, slip wall, as shown in Fig. 5. Some inlet parameters for 5 T3C cases are given in Table 1. It can be seen that T3C1 has a high inlet freestream turbulent intensity, while other cases have moderate turbulence levels.

Skin friction C_f profiles for T3C cases are compared with experimental data in Fig. 6. Turbulence intensity (Tu(x)) decay agrees with the experiment data, it is not plotted here, for brevity. Both C_f and Re_x are scaled by the local freestream velocity U_fs (taken from y = 0.05 m, similar to Ref. [50]).

For T3C1, with high inlet Tu, the transition location and length are predicted well by the current model. For T3C2, with moderate Tu, transition occurs a bit later than T3C1. The predicted of skin friction is in better agreement with the experimental data than the result of Ge et al. [17]. The T3C3 case has smaller U_in and Tu_in than T3C2, which delays transition to near the end of the plate. Again, the current model agrees well with the experiment. T3C4 further decreases U_in while maintaining similar level of Tu_in as T3C3. While T3C4 was considered as a benchmark for separation-induced transition in Ref. [17], it leaves very little room near the end of the plate for the transition to turbulence, i.e., just one experimental data point. The current model has similar level of agreement with the experimental data as Ref. [17]. T3C5 has the highest U_in and second highest Tu_in, which makes its transition location move upstream, to the favorable pressure gradient region. Compared with the experiment, the current prediction seems to be good.

Transition occurs on compressor or turbine blades that are subjected to impinging turbulent wakes. In order to remove the influence of pressure gradients and surface curvature, periodic wakes were swept across a flat plate in the study by Wu et al. [51,52]. The computational domain for their DNS data is x = [0.1, 3.5], y = [0, 0.8] in nondimensional units. The inlet velocity is U_ref = 1 plus wakes. The Reynolds number is $Re=UrefL/ν=1.5×105$⁠. The boundary layer starts from a Blasius profile at the inlet. Self-similar wakes with half width of b = 0.1, velocity deficit of U_def = 0.14, and translational velocity of 0.7 are periodically swept across the inlet.

Similar to the flat plate test cases in Secs. 3.1 and 3.2, the computational domain of the current study includes the leading edge of the plate; that is, x = [–0.05, 3.5], y = [0, 0.8]. The inlet reference velocity and Reynolds number, as well as the inlet wake period $T=1.67$ are the same as Ref. [52]. The inlet wake parameters such as U, k, and ω are adjusted so that the respective parameters at x = 0.1 agree with those of Wu and Durbin [52]. U_def = 0.48 and b = 0.03 were chosen for the current inlet (x = –0.05) wake U, k, and ω profiles. Wake profiles at x = –0.05 and x = 0.1 are compared with Ref. [52] in Fig. 7. The general agreement is quite acceptable.

Space-time contours of phase-averaged skin friction coefficient C_f(x, t), predicted by the current model and the DNS data [52], are shown in Fig. 8. C_f profiles along the plate at a fixed time instant are horizontal sections of the contour plot. It can be seen that the current model results agree well with the DNS data. The blue tongue-shaped regions around x = 1 represent low C_f in the laminar region before transition. The surrounding yellow, wavy contours, around x = 1.5, represent the end stage of the transition, right before turbulence begins. Due to the strong influence of the period wakes, the laminar skin friction and the transition location fluctuate with time; and, at a certain time instant (e.g., at phase = 0.75), the passing wake leads to bell-shaped, dark blue, protrusions of the laminar skin friction.

A curve of the time-averaged skin friction coefficient, which is obtained by averaging over all phases, is plotted in Fig. 9. It is compared to the DNS data of Wu an Durbin [52]. Both curves agree on the transition location and on its development. Compared with the DNS data, C_f in the turbulence region seems slightly under predicted for the current model, which is not unexpected: the RANS model is not accurate at the low Reynolds number just after transition.

The Stripf's [34,47,48] low pressure turbine (LPT) was selected for validation tests. Note that Boyle and Stripf [34] proposed a correlation formula to compute equivalent sandgrain roughness from the roughness geometries statistics, e.g., root-mean-square and skewness, to validate their model. The current model uses the equivalent sandgrain roughness computed using their method, for comparison purposes.

The computational domain is shown in Fig. 10. The true chord c = 113.34 mm; U_in = 33.086 m/s; Tu_in = 3% and Re_in = 250,000. The roughness height ranges from 0 to 395 μm (fully rough). The blade temperature is set as 300°K and the fluid 400°K. ρ = 1.2 kg/m³, C_p = 1000 m²/s²K. The laminar and turbulent Prandtl numbers are set as Pr = 0.72 and Pr_T = 0.86, respectively. An effective thermal diffusivity κ_eff = ν/Pr + ν_T/Pr_T is used in the additional energy equation, to represent heat transport.

The inlet turbulence kinetic energy k_in is mainly based on the inlet turbulence intensity (Tu_in), and the inlet specific dissipation rate ω_in is adjusted based on the Tu decay profile. As indicated in Ref. [18], however, there is some uncertainty regarding the freestream turbulence around the turbine blade. The k_in and ω_in in this study were adjusted based on the Nu distribution for the smooth surface. Then, the influence of the roughness is modeled by changing the equivalent sandgrain roughness (r). The current Tu_in is 1.5%, which is similar to Ref. [18].

The Nusselt number is defined as Nu = hc/κ, where h is the heat transfer coefficient, and κ is the thermal conductivity. Nu distributions along the suction side of the blade are plotted in Fig. 11, compared with the experimental data. The x coordinate is the surface distance normalized by the true chord (s/c). It can be seen that the profiles are in good agreement to data on both transition onset location and transition length, for the various roughness. Compared with the smooth surface, heat transfer on the rough surfaces rises earlier, and is significantly enhanced. Note that the profiles for two high roughness cases are a little over-predicted right after transition. It may be due to the turbulence model with Knopp boundary conditions, which is suitable for fully developed turbulence flows, while the current turbulent flow right after transition is not fully developed yet.

The computational domain of the two-dimensional incompressible flow through the passage of the HPT blade is shown in Fig. 12. The true chord is c = 93.95 mm. The top and bottom boundaries of the computation are chosen in such a way that a cyclic condition can be imposed.

The inlet velocity is U_in = 39.915 m/s and the inlet Reynolds number Re_in = 250,000 based on the true chord. The blade surface roughness ranges from 0 to 320 μm (fully rough). The thermal properties such as temperature, density, Prandtl numbers are the same as the LPT case.

The heat transfer coefficient (h = q/ΔT where q is the heat flux, ΔT = 400 − 300 = 100) on the blade suction side predicted by the current model is compared with the reference experimental and model data in Fig. 13. It was found that the transition location moves from near the trailing edge, to the leading edge of the blade, as roughness increases from r = 0 μm to r = 320 μm. For all the roughness considered, the transition location predicted by the current model is in reasonable accord with the reference model and experimental data. After transition, however, the model indicates that the increasing roughness causes an increase of h in the fully turbulence regime, while the measured level of h, after transition, is insensitive to roughness. This might be a discrepancy of the k–ω predictions. The RANS prediction seems to make sense—that higher roughness increases mixing in a turbulent boundary layer, leading to higher heat transfer.

A new, intermittency function-based transition model was introduced. It is based on local flow variables, and was developed for low, moderate, and high freestream turbulence intensity, over smooth or rough surfaces. It uses a single transport equation for the intermittency function. In order to satisfy requirements for general-purpose CFD, nonlocal operations, such as boundary layer integral thickness, were avoided in the current model. The roughness model was formulated as a modified, displacement of origin, version of the k–ω RANS model.

Validation tests were performed for various types of flows over smooth/rough surfaces, including the ERCOFTAC, experimental, zero-pressure-gradient and favorable/adverse-pressure-gradient boundary layer cases, periodic impinging wakes, and turbine blades with regular roughness. The predicted skin friction agreed reasonable well with the experimental data. However, under low free-stream turbulence levels the model incorrectly produced a small increase in the level of the laminar skin friction. This is an inadvertent effect of a modification for separated flow.

The model was tested on periodic wakes impinging on a flat plate boundary layer, to study wake-induced transition, which often occurs on compressors and turbine blades. Both phase-averaged skin friction contours and time-averaged skin friction profiles, predicted by the current model, agreed well with the reference DNS data.

The roughness model was validated on Stripf's low and high pressure turbine blades, with a range of roughness heights. The Nusselt number and heat transfer coefficient along the blade surface, predicted by the current model, are in agreement with experimental data. After transition, the suction side heat transfer is over predicted. This may be because of the turbulence model with Knopp boundary conditions, which is suitable for fully developed turbulence flows, while the current turbulent flow right after transition is not fully developed yet. In summary, the model prediction for both the transition onset location and transition length is decent.

• National Science Foundation (Award No. CBET 1228195).

• Small Business Innovation Research, U.S. Army SBIR Phase II research program (Contract No. W911W6-14-C-0003).

C_k, C_ω1, C_ω2, σ_ω, σ_k =

k–ω model parameters

C_r =

roughness effects coefficient

d =

wall distance

D_γ =

diffusion term

$d0+$ =

hydrodynamic roughness function

E_γ =

sink term

F_Q =

pressure gradient parameter

$FRt,FRs$ =

effective intermittency function limiters

F_γ =

source term parameter

F_turb =

sink term parameter

G_γ =

sink term parameter

k_w, ω_w =

k, ω at the wall

n_w =

unit wall normal vector

P_γ =

source term

Q =

Ω_ijΩ_ij − S_ijS_ij

r =

equivalent sandgrain roughness

R_r =

local variable-scaled equivalent sandgrain roughness

R_s =

adverse pressure gradient detector

R_t =

ν_T/ν

Re_c1 =

critical Reynolds number

Re_c2 =

roughness-induced Reynolds number

Re_c3 =

roughness/pressure gradient-induced Re

Re_θt =

momentum thickness-based critical Reynolds number

Re_ν =

wall distance-based vorticity Reynolds number

Re_νnew =

vorticity Reynolds number

S =

mean strain rate $|S|=SijSji$

T_ω =

turbulence level parameter

Tu =

freestream turbulence intensity

γ =

intermittency function

γ_eff =

effective intermittency function in k equation

γ_max =

parameter to strengthen source term

ν =

kinematic viscosity

ν_t =

turbulence viscosity k/ω

σ_l =

empirical coefficient for laminar viscosity

σ_γ =

empirical coefficient for turbulence viscosity

ω =

vorticity

Ω =

mean rotation rate $2ΩijΩij$