Abstract
Spanwise oscillation applied on the wall under a spatially developing turbulent boundary layer flow is investigated using direct numerical simulation. The temporal wall forcing produces a considerable drag reduction over the region where oscillation occurs. Downstream development of drag reduction is investigated from Reynolds number dependency perspective. An alternative to the previously suggested power-law relation between Reynolds number and peak drag reduction values, which is valid for channel flow as well, is proposed. Considerable deviation in the variation of drag reduction with Reynolds number between different previous investigations of channel flow is found. The shift in velocity profile, which has been used in the past for explaining the diminishing drag reduction at higher Reynolds number for riblets, is investigated. A new predictive formula is derived, replacing the ones found in the literature. Furthermore, unlike for the case of riblets, the shift is varying downstream in the case of wall oscillations, which is a manifestation of the fact that the boundary layer has not reached a new equilibrium over the limited downstream distance in the simulations. Taking this into account, the predictive model agrees well with DNS data. On the other hand, the growth of the boundary layer does not influence the drag reduction prediction.
1. Introduction
Turbulence in an ubiquitous phenomena in aerospace applications is, more often than not, detrimental for the performance of, for example, aircrafts. In particular, the skin friction, which produces about half of the total drag, continues to be a challenge to decrease by means of flow control. Strategies for achieving drag reduction are broadly classified as either active or passive. Active mechanisms are further classified as closed loop or open loop type, depending on whether feedback from the flow is used for calculating the control actuation. Considerable progress both in the understanding of drag reduction mechanisms and in the development in techniques for obtaining drag reduction has been made in the past ten years. This progress has been made possible by the advancement of computer capabilities and more or less imaginative control strategies. Spanwise wall oscillations are an active mechanism for achieving high drag reductions which has the attractive property of being open loop and thus lends itself suitable for practical implementation without the need for complicated sensors to detect the flow state.
However, two outstanding issues remain before any active control mechanism can be implemented in real engineering applications. Firstly, the power to drive the control is still much higher than the power gained from drag reduction. This may be addressed by novel techniques replacing the actual motion of the wall, for example, plasma actuators [1, 2], Lorentz force [3, 4], or dielectric electroactive polymers [5]. Secondly, when moving from low Reynolds number in the computer or laboratory setting to the high Reynolds number in real applications, the effect on the drag reduction efficiency is not known. This second issue is in the context of aerospace applications equally important as the first one and is the topic of the present investigation.
Ever since the first studies of turbulent flow in a channel with wall oscillations using direct numerical simulation (DNS) by Jung et al. [6], the topic has been slowly progressing towards higher Reynolds numbers [7ā9]. The latest work on this topic has been completely focused on channel flow, including the more analytical work by [10]. For a list of references to work done between the pioneering study and the later ones mentioned above, please refer to, for example, [11] which utilizes data from many different investigations.
The influence of Reynolds number on DR is an important question for practical implementation aspects. This is still an open question and a clear understanding of the Reynolds number dependence of DR is presently lacking. Results from previous channel flow investigations with near-optimal oscillation parameters suggest a relation of the formwhere is the Reynolds number based on channel half width (), friction velocity (), and kinematic viscosity (). This relation was first suggested by Choi et al. [12]. Channel flow simulations with optimal oscillation parameters have subsequently suggested that . However, as indicated by Gatti and Quadrio [13], the value of the exponent in relation (1) is most probably dependent on the oscillation parameters; thus a universal relation will be unattainable. In addition, in the study of the Reynolds number averaged Navier-Stokes equations using a perturbation analysis by Belan and Quadrio [14], it is found that a relation of the typematches the data better. Hence, the DR approaches a constant value as .
Most of the previous DNS studies have focused on channel flow with the benefits of periodic boundary conditions which reduce the computational costs. However, not much work has been done to explore the effects on a turbulent boundary layer using numerical simulations. Note that the local drag reduction decreases downstream for a spatially growing boundary layer flow in contrast to the channel flow where drag reduction remains constant in both time and space after attaining an equilibrium value. The same numerical code used for the present investigation has previously been used for both temporal and spatial oscillations of the wall under a turbulent boundary layer [15ā17]. Recently, Lardeau and Leschziner [18] reported their findings for boundary layer flows using DNS where they focused on the transition of the skin friction to a steady low-drag state and perform phase-averaged analysis of turbulent statistics.
In passing we also note that the wall oscillation technique has been utilized for transition control in the boundary layer; see [19ā21].
Regarding practical applications, the oscillating wall technique may be still far from being matured. With the existing technological capabilities available, it is hard to imagine large portions of the aircraft wing or fuselage surface oscillating. However, since other methodologies, based on passive control such as riblets, have yet to be commonplace, the prospect of other methods based on active control, possibly refined by feedback control in order to increase their energy saving capacity, should be considered for future applications. In passing, one may furthermore note that the field of active flow control is open to new innovative ideas, although based on oscillating wall motion. One example is the rotating disc idea pursued by Ricco and his group [22, 23].
Although laminar flow is prevalent in many applications, for example, biochemistry [24] or biology [25], which allows for numerical simulations using the Navier-Stokes equations, turbulent flow is always encountered in aerospace applications. If realistic Reynolds numbers are considered when performing numerical work, the Navier-Stokes equations will not be possible to solve due to the wide range of scales (as measured by the Reynolds number) present in the flow. Hence, turbulence modelling will be required and either the Reynolds Average Navier-Stokes equations (RANS) or Large Eddy Simulations (LES) need to be employed. However, for the phenomena of drag reduction due to wall forcing, as investigated in this work, there is no reliable turbulence model available. Hence, we focus on low Reynolds number flow which, although turbulent, can be numerically simulated using DNS, that is, by obtaining the solution to the Navier-Stokes equations without any turbulence modelling. Only after understanding how the turbulence is affected by the wall forcing, a model which can be used for higher Reynolds number engineering application may be developed. For the purpose of DNS, highly efficient and accurate numerical schemes based on spectral methods have been employed in this work.
The rest of the paper is organized as follows. In Section 2, the numerical scheme and simulation parameters are presented. The results are discussed in Section 3, which first concerns older data found in the literature in Section 3.1. Based on previous simulations and experiments, a predictive relation for DR at high Reynolds numbers is suggested, which differs from the general conclusion of a Reynolds number dependency according to (1). Then, in Section 3.2, a predictive formula based on the shift of the velocity profile in drag-reduced flows is derived and compared with the new simulation results. In Section 3.3 the velocity profile shift is further investigated and the implications of the results are related back to the findings in Section 3.2. Lastly, conclusions are mentioned in Section 4.
2. Numerical Method and Simulation Parameters
The numerical code used, SIMSON (A Pseudo-Spectral Solver for Incompressible Boundary Layer Flows), was originally developed at KTH, Stockholm [31].
An outline of the numerical scheme is presented in Section 2.1; the various parameters used and the resolution are presented next in Section 2.2, where also the implementation of the wall motion is discussed.
2.1. Numerical Scheme
A pseudospectral method with Fourier discretization in the streamwise and spanwise directions and Chebyshev polynomials in the wall-normal direction has been used. The simulations are started with a laminar boundary layer at the inflow. A random volume force near the wall at the beginning of the computational domain is used to trigger the flow to transition.
At the end of the domain, a fringe region is added to enable simulations of spatially developing flows. The flow in this region is forced from the outflow of the physical domain to the inflow. In this way the physical domain and the fringe region together satisfy periodic boundary conditions. The implementation is done by adding a volume force, to the Navier-Stokes equations:
The force applies to the fringe region only where is the strength of the forcing and is the laminar inflow velocity profile to which the solution is forced. The fringe function is required to have minimum upstream influence and is designed as with
Here is the maximum strength of the fringe, and denote the spatial extent of the region where the fringe is nonzero, and and are the rise and fall distance of the fringe function, respectively. Figure 1 shows how the fringe function varies with with and . is a continuous step function that varies from zero for to unity for and is given byThe advantage of this expression is that has continuous derivatives of all orders.
The function is also used in enforcing the wall oscillation boundary condition as described in Section 2.2.
The time integration is performed using a third-order Runge-Kutta scheme for the nonlinear terms and a second-order Crank-Nicolson method for the linear terms. A 3/2-rule is applied to remove aliasing errors from the evaluation of the nonlinear terms when calculating FFTs in the wall parallel plane.
2.2. Numerical Parameters
All quantities are nondimensionalized by the free stream velocity () and the displacement thickness () at the starting position of the simulation (), where the flow is laminar. The Reynolds number is set by specifying at the laminar inlet ().
Wall oscillation is imposed over a long segment of the wall to examine the influence of (, being the momentum thickness) on DR. Table 1 gives the different parameters chosen for the computational box and grid. The length of the computational box is given in simulation length units (). The resolution used for the simulations is shown in + units. Note that unless otherwise stated the + superscript indicates that the quantity is made to be nondimensional with the friction velocity of the unmanipulated boundary layer (the reference case), denoted as , and the kinematic viscosity (). Note that is varying with . In the calculation of the resolution, a representative is taken from a position , where is the length of the computational box, which corresponds to the location where the oscillating region ends. The grid resolution is consistent with other DNS, which has been proven to be sufficient [29]. In addition, with the same resolution and Reynolds number another case was simulated for which experimental data exist [28] and almost perfect agreement was obtained with respect to both magnitude and spatial development; see [26].
An illustration for the current problem setup is shown in Figure 2. After an initial trip of the flow, the flow undergoes transition to fully turbulent flow. The tripping has been described in detail in [29] and the quality of this methodology has been thoroughly investigated by Schlatter and ĆrlĆ¼ [30]. A section of the resulting turbulent flow is then subjected to wall oscillations, as described below. At the end of the domain, there is a fringe region to satisfy periodic boundary condition as discussed in the previous section. Minimal upstream influence is ensured by the design of the fringe function [31].
The implementation of the oscillating spanwise wall velocity is identical to the one used in [29]. The wall oscillation is applied in the spanwise direction at a particular region in streamwise direction. Therefore, a profile function is utilized to define the domain where the oscillation occurs. Furthermore, the gradual increase of the profile function prevents the spurious oscillations (Gibbās phenomena) which might otherwise be introduced due to a discontinuous jump in the velocity around the starting point of wall forcing.
The form of the wall velocity boundary condition is given bywhere is the same profile function as that used for fringe region; see (6); and is the maximum wall velocity. The parameter is the angular frequency of the wall oscillation, which is related to the period through . The parameter determines when the oscillation starts.
The sampling time for the reference case was in time units () and started only after a stationary flow (in the statistical sense) was reached after time units. In the cases with wall forcing, the reference case was used as the initialized flow field and additional 4300 time units were simulated before statistics were sampled during another 33000 time units, corresponding to 478 periods of oscillation. The parameter in (6) is 38000 ().
In the simulations presented here, the angular frequency () of the wall oscillation is set to and maximum wall velocity () is set as , which in wall units corresponds to and , calculated as and , based on at the start of oscillations (), which is equal to . Note that the value of and changes with streamwise position because local varies, whereas and are kept fixed. At the end of the oscillating part () the friction velocity of the reference case is and gives a value of and at that position. The parameters pertaining to the oscillations are summarized in Table 2.
3. Results
Before presenting the results in Sections 3.1 to 3.3, various definitions needed in the following are presented. The friction velocity is defined aswhere is the kinematic viscosity.
The resulting drag reduction (DR) is calculated fromwhere is the skin friction of the reference case.
An alternative expression of the DR will in this paper be used according to
Hence, we may write
3.1. Drag Reduction in Channel and Boundary Layer Flows
In this section focus will be on the maximum DR in the boundary layer and its dependency on the Reynolds number where the oscillation starts. In addition, data from channel flow DR in the literature will be utilized as there are more data to compare with for this flow geometry. Note that one unique DR value exists for a channel flow, while the DR is varying downstream for a boundary layer. Hence, in order to compare the two cases, we specify the maximum DR in boundary layer case for comparison with channel flow data.
Three different sets of channel flow data will be presented here. They are chosen because up to four different Reynolds numbers have been investigated for identical oscillation parameters.
For comparison with the boundary layer it is necessary to convert the Reynolds number used in channel flow configuration, , to , which is done by using the relationas given by Schlatter and ĆrlĆ¼ [32].
The data is presented in Figure 3 and is colour-coded according to the oscillation parameters used in the DNS/experiments. Note that black, red, and magenta constitute separate sets of parameters, while green, blue, and cyan are data from investigations using identical values of and . In addition, data from DNS/experiments is indicated by symbols, while the solid and broken lines represent the relations proposed.
The boundary layer data shown in Figure 3 as circles is taken from [26] and is compared with channel flow (cross) from the DNS by [27] as well as experimental boundary layer data (plus) from [28]. The DNS data from [26] show good agreement with the other channel flow simulations and boundary layer experiments at identical oscillation parameters of and . Note that the simulations and experiments of the boundary layer were too short in those investigations to provide any significant information about the variation downstream with increasing . Thus, we treat the maximum DR obtained from the boundary layer in these cases as equivalent to the uniquely obtained value of DR for the channel flow. (In contrast, the present boundary layer simulation, shown with triangles in Figure 3, is much longer and the downstream variation is discussed in the next two sections.)
Connect all the data points (in black colour in Figure 3) with a curve obtained aswhere and . This relation suggests that as becomes larger, DR approaches a constant value of 23.44% for this particular set of oscillation parameters.
Another two data sets have been produced at higher Reynolds numbers for the channel flow geometry with identical oscillation parameters at and . Touber and Leschziner [7] (diamond) reached , while Hurst et al. [9] increased to (squares). It is interesting to note the very different resulting values of the DR in these two cases. Both of these research groups have claimed that a power law relation, as proposed by Choi et al. [12] according to (1), fits their data with and , indicated with the broken lines in Figure 3. These profiles are derived by using the relation between and above, (see (13)), to calculate the values of in the expression,which are given in Table 3. On the other hand, the herein proposed relation (14) results in slightly more accurate fit of the data. For the parameters derived, please see Table 3. In addition, one more channel flow simulation, taken from [27] with and , is shown (cross) and differs as well from the two other sets of data (even if only a single Reynolds number is available from this simulation).
Choi et al. [12] have provided results for DR in a channel for different with and , and the data is presented (stars) in Figure 3. The curve fitting in this case yields and in expression (14). The fact that the coefficients in (14) depend on the oscillation parameters is consistent with the investigation of the power law relating DR with Reynolds number for channel flow. Furthermore, the analysis by Belan and Quadrio [14] indicated clearly that a constant DR is obtained in the limit of very high Reynolds numbers, which is also consistent with findings of Iwamoto et al. [33].
After this discussion on data found in the literature we turn to the present data, indicated as triangles in Figure 3. This set of data is different from the others in that the downstream variation is available, not only the maximum DR. There seems to be an agreement with the power law in this case, which will be discussed in more detail in the next section.
Finally, we may observe that although the rapid decay of DR at low Reynolds numbers is captured by the power law, the trend at higher Reynolds number where the values seem to reach a constant value has no correspondence to the power law which instead indicates a continuous decaying DR at high Reynolds numbers. In addition, the power law suggests that the DR approaches zero in the limit of infinite Reynolds number, which one might find peculiar. However, this is consistent with other empirical formulas such as the relation between skin friction and Reynolds number,proposed by [34] which indicates vanishing skin friction as the Reynolds number grows. In short, because there is no drag to reduce when .
3.2. Drag Reduction of the Boundary Layer and Dependency on the Reynolds Number
In order to explain the decreasing DR with Reynolds number, two different linear approximations have been used in the past, which have been focused on the fact that the upward shift of the velocity profile () is constant, regardless of Reynolds number. The shift is equal to the change in intercept of the logarithmic law. The connection between the DR and upward shift of the velocity profile is of course nothing new since it has been observed for a long time. In fact, simplified relations between the shift and the DR, based on various degree of approximations, have been proposed. For example, [35] derived an expression,by differentiating the logarithmic friction law. This expression can also be found in [36] and takes into account the growth of the boundary layer.
A simplified version was used by [37] which yields
The two expressions above were derived in the context of DR by means of riblets, in which case is constant and linearly related to the protrusion height. Furthermore, in the case of (18), by neglecting the growth of the boundary layer, we can express the velocity profile shift assince the wake function is assumed to be unaltered by the DR. Equation (19) can be interpreted as an equating of the shift of the velocity profile in the logarithmic region with the shift of the wake region.
In order to test two expressions (17) and (18) above, we simply insert the spatially varying according to (19) from the DNS into the expressions and plot them in Figure 4. From Figure 4 we note that expressions (17) and (18) yield very similar results, probably due to the limited variation of the Reynolds number, which translates to a small variation of the boundary layer thickness.
In the following a new expression for relating DR and is derived by observing that the right-hand side of (19) constitutes a linearization of the DR.
The starting point is to rewrite (19) as
Multiplying (20) with yields
We wish to replace with on the right-hand side, so (20) is rewritten asfrom which we obtainwich, inserted in the right-hand side of (21), gives
Using Taylor series expansion of the left-hand side of (24) and neglecting higher order terms, we obtain
Utilizing the approximative relation (25), the following expression for the DR is obtained:
The DR obtained from (26) is illustrated by the blue broken curve in Figure 4 and is more accurate than the other two expressions ((17) and (18)) found in the literature. Note that by neglecting the first term on the right-hand side in (22), the original expression (18) is reproduced. Hence, terms of the order are neglected in the original approach, while in the current derivation only higher order terms are neglected by the Taylor series truncation in (25).
The equivalent expression, taking the growth of the boundary layer into account, may also be obtained using the same technique applied to the logarithmic friction law and iswhich contains an identical correction as compared to expression (17) by a factor in the denominator. However, similar to the original expressions, the inclusion of the growth of the boundary layer via the term has very little influence on the result (not shown in the figure).
3.3. The Velocity Shift for the Boundary Layer Flow
In the present case, where the DR is obtained by the oscillating wall technique, needs to be investigated. The downstream variation of the velocity shift according to (19) is shown in Figure 5 as the black curve, and is clearly not a constant. However, when calculating the change in the intercept of the logarithmic layer, the values seem to reach a constant value, as indicated by the blue curve. The intercept is here calculated by considering the indicator function,which for low Reynolds number flow never reaches a constant value, although a distinct local minimum is always found. By finding the local minimum and then at that position calculating the difference between the velocity from the reference and the drag reduced profiles, the shift of the logarithmic profile is derived. This shift is shown as the blue line in Figure 5 and is approaching a constant value. The red line in the figure is the constant value taken from the expression,derived from the relation in [11] where the constant was found to be . With the value of in the present case, expression (29) yields . from expression (19) is slowly approaching the value given by the velocity profiles, which indicates that the wake function is slowly reaching a new equilibrium. This scenario is consistent with the findings in [11].
Finally, to investigate if the constant value of yields the desired variation of DR with Reynolds number by utilizing (26), the DR is plotted as the blue solid line in Figure 4. It does not seem unlikely that all three curves in Figure 5 converge to the same value at high Reynolds number (futher downstream), which would in turn yield a correct prediction of the DR and its variation with through (26).
Note that, by using (16) in (26) or (27), we obtainfor high Reynolds numbers.
This variation is shown by the cyan solid line in Figure 4 to indicate that this indeed gives identical Reynolds number variation as (26).
4. Conclusions
Direct numerical simulations have been performed to study the effect of wall oscillation on the spatially developing turbulent boundary layer. The resulting drag reduction and its variation in the streamwise direction are investigated by comparison with earlier short boundary layer simulations and experiments as well as channel flow simulations. In addition, new relations regarding the prediction of the DR at higher Reynolds numbers have been formulated.
Two separate relations between DR and Reynolds number are considered here. One is the DR which refers to the uniquely determined DR in channel flow or the maximum DR obtained in the boundary layer flow. The other is the downstream development of DR in a boundary layer flow.
The dependency of peak drag reduction values on the Reynolds number at which the oscillation commences is similar to that of channel flows, although a conclusive relation between DR and Reynolds number seems elusive as the relation most probably varies for different oscillation parameters. On the other hand, we have in the present work diversified the analysis and have shown that an alternative relation of type (14) can describe the current set of available data. The relation given here suggests that drag reduction attains a constant value at higher Reynolds numbers. The conclusions that can be drawn from this part of the study are the following:(i)The influence the Reynolds number has on the drag reduction varies widely in the data taken from internal channel flows in literature, even with identical oscillation parameters and flow geometry.(ii)Both the traditional relation (15) and the herein proposed expression (14) can be made to fit the data.
A relation between DR and the shift of the logarithmic velocity profile is derived which differs from the expressions found in literature. The conclusions drawn are the following:(i)The herein derived equation (26) yields DR values in more agreement with DNS data compared to (18).(ii)The growth of the boundary layer does not need to be taken into account when estimating the drag reduction of a turbulent boundary layer.
A new estimation of the DR variation for a spatially developing boundary layer was derived by combining (26) or (27) at high Reynolds numbers with the empirical relation (16), giving . Regarding the presented DNS data of the downstream behavior of the present long simulation of the drag reduced boundary layer, the following conclusions can be drawn:(i)The DR is influenced by spatial transient even for this long computational simulation and the theoretically obtained (see above) relation is only slowly approached downstream.(ii)The shift of the velocity profile according to (19) is not constantly downstream, probably due to the fact that the wake function is still adjusting to the new turbulent stage near the wall.(iii)The shift of the velocity profile calculated from the velocity difference at the position of the local minimum of the indicator function (28) is approaching a constant value identical to expression (29) from [11].
Further validation of our numerical results has been planned to use state-of-the-art three-dimensional PIV measurements [38ā41].
Conflict of Interests
The authors declare that there is no conflict of interests regarding the publication of this paper.
Acknowledgment
This research work is supported in part by Singapore MOE Tier-2 grant MOE2012-T2-1-030.