# Drag Reduction of a Turbulent Boundary Layer over an Oscillating Wall and Its Variation with Reynolds Number

Authors: Martin Skote; Maneesh Mishra; Yanhua Wu
Journal: International Journal of Aerospace Engineering
 (2015)

---

## 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.


---

## Body

## 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 form(1)
D
R
~
R
e
τ
α
,where R
e
τ
=
h
u
τ
/
ν is the Reynolds number based on channel half width (h), friction velocity (u
τ), and kinematic viscosity (ν). This relation was first suggested by Choi et al. [12]. Channel flow simulations with optimal oscillation parameters have subsequently suggested that α
=
-
0.2. 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 type(2)
D
R
~
γ
+
R
e
τ
αmatches the data better. Hence, the DR approaches a constant value as R
e
τ
→
∞.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 Section2, 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 Section2.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,(3)
F
i
=
λ
x
u
~
i
-
u
i
, to the Navier-Stokes equations:(4)
∂
u
i
∂
t
+
u
j
∂
u
i
∂
x
j
=
-
1
ρ
∂
p
∂
x
i
+
ν
∂
2
u
i
∂
x
j
2
+
F
i
.The forceF
i applies to the fringe region only where λ
(
x
) is the strength of the forcing and u
~
i is the laminar inflow velocity profile to which the solution u
i is forced. The fringe function is required to have minimum upstream influence and is designed as(5)
λ
x
=
λ
max
⁡
f
x
, with (6)
f
x
=
S
x
-
x
start
Δ
x
rise
-
S
x
-
x
end
Δ
x
fall
+
1
.Hereλ
max
⁡ is the maximum strength of the fringe, x
start and x
end denote the spatial extent of the region where the fringe is nonzero, and Δ
x
rise and Δ
x
fall are the rise and fall distance of the fringe function, respectively. Figure 1 shows how the fringe function varies with x with Δ
x
rise
=
3 and Δ
x
fall
=
3. S
(
η
) is a continuous step function that varies from zero for η
≤
0 to unity for η
≥
1 and is given by(7)
S
η
=
0
,
η
≤
0
,
1
1
+
e
1
/
η
-
1
+
1
/
η
,
0
<
η
<
1
,
1
,
η
≥
1
.The advantage of this expression is that S
(
η
) has continuous derivatives of all orders.Figure 1
Functionf
(
x
) with Δ
x
r
i
s
e
=
Δ
x
f
a
l
l
=
3.The functionf
(
x
) 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 (U
∞) and the displacement thickness (δ
∗) at the starting position of the simulation (x
=
0), where the flow is laminar. The Reynolds number is set by specifying R
e
δ
∗
=
U
∞
δ
∗
/
ν at the laminar inlet (x
=
0).Wall oscillation is imposed over a long segment of the wall to examine the influence ofR
e
Θ (R
e
Θ
=
U
Θ
/
ν, Θ 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 (δ
x
=
0
∗). 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 u
τ
0, and the kinematic viscosity (ν). Note that u
τ
0 is varying with x. In the calculation of the resolution, a representative u
τ
0 is taken from a position x
=
2
/
3
L
x, where L
x 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].Table 1
Numerical parameters for the simulations presented.L and N refer to the size (measured in displacement thickness at x
=
0) and the number of modes, respectively.
L
x
L
y
L
z
N
x
N
y
N
z
Δ
x
+
Δ
y
m
i
n
+
Δ
z
+
2400
75
34
3200
301
144
15.8
0.04
4.96An illustration for the current problem setup is shown in Figure2. 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].Figure 2
Illustration for the problem setup.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 f
(
x
) 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 by(8)
W
y
=
0
=
f
x
W
m
sin
⁡
ω
t
-
t
start
,where f
(
x
) is the same profile function as that used for fringe region; see (6); and W
m is the maximum wall velocity. The parameter ω is the angular frequency of the wall oscillation, which is related to the period through ω
=
2
π
/
T. The parameter t
start determines when the oscillation starts.The sampling time for the reference case was24800 in time units (δ
x
=
0
∗
/
U
∞) and started only after a stationary flow (in the statistical sense) was reached after 13200 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 t
start in (6) is 38000 (=
13200
+
24800).In the simulations presented here, the angular frequency (ω) of the wall oscillation is set to 0.09099 and maximum wall velocity (W
m) is set as 0.5028, which in wall units corresponds to T
+
=
76.4 and W
m
+
=
10.7, calculated as W
m
+
=
W
m
/
u
τ
0 and T
+
=
T
(
u
τ
0
)
2
/
ν, based on u
τ
0 at the start of oscillations (x
=
500), which is equal to 0.047. Note that the value of W
m
+ and T
+ changes with streamwise position because local u
τ
0 varies, whereas W
m and T are kept fixed. At the end of the oscillating part (x
=
1600) the friction velocity of the reference case is u
τ
0
=
0.042 and gives a value of W
m
+
=
12.0 and T
+
=
60.8 at that position. The parameters pertaining to the oscillations are summarized in Table 2.Table 2
Parameters regarding the wall oscillations for the simulation presented.
W
m
+
T
+
x
s
t
a
r
t
x
e
n
d
R
e
Θ
s
t
a
r
t
R
e
Θ
e
n
d
10.7
76.4
500
1600
963
1989


## 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,(3)
F
i
=
λ
x
u
~
i
-
u
i
, to the Navier-Stokes equations:(4)
∂
u
i
∂
t
+
u
j
∂
u
i
∂
x
j
=
-
1
ρ
∂
p
∂
x
i
+
ν
∂
2
u
i
∂
x
j
2
+
F
i
.The forceF
i applies to the fringe region only where λ
(
x
) is the strength of the forcing and u
~
i is the laminar inflow velocity profile to which the solution u
i is forced. The fringe function is required to have minimum upstream influence and is designed as(5)
λ
x
=
λ
max
⁡
f
x
, with (6)
f
x
=
S
x
-
x
start
Δ
x
rise
-
S
x
-
x
end
Δ
x
fall
+
1
.Hereλ
max
⁡ is the maximum strength of the fringe, x
start and x
end denote the spatial extent of the region where the fringe is nonzero, and Δ
x
rise and Δ
x
fall are the rise and fall distance of the fringe function, respectively. Figure 1 shows how the fringe function varies with x with Δ
x
rise
=
3 and Δ
x
fall
=
3. S
(
η
) is a continuous step function that varies from zero for η
≤
0 to unity for η
≥
1 and is given by(7)
S
η
=
0
,
η
≤
0
,
1
1
+
e
1
/
η
-
1
+
1
/
η
,
0
<
η
<
1
,
1
,
η
≥
1
.The advantage of this expression is that S
(
η
) has continuous derivatives of all orders.Figure 1
Functionf
(
x
) with Δ
x
r
i
s
e
=
Δ
x
f
a
l
l
=
3.The functionf
(
x
) 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 (U
∞) and the displacement thickness (δ
∗) at the starting position of the simulation (x
=
0), where the flow is laminar. The Reynolds number is set by specifying R
e
δ
∗
=
U
∞
δ
∗
/
ν at the laminar inlet (x
=
0).Wall oscillation is imposed over a long segment of the wall to examine the influence ofR
e
Θ (R
e
Θ
=
U
Θ
/
ν, Θ 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 (δ
x
=
0
∗). 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 u
τ
0, and the kinematic viscosity (ν). Note that u
τ
0 is varying with x. In the calculation of the resolution, a representative u
τ
0 is taken from a position x
=
2
/
3
L
x, where L
x 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].Table 1
Numerical parameters for the simulations presented.L and N refer to the size (measured in displacement thickness at x
=
0) and the number of modes, respectively.
L
x
L
y
L
z
N
x
N
y
N
z
Δ
x
+
Δ
y
m
i
n
+
Δ
z
+
2400
75
34
3200
301
144
15.8
0.04
4.96An illustration for the current problem setup is shown in Figure2. 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].Figure 2
Illustration for the problem setup.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 f
(
x
) 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 by(8)
W
y
=
0
=
f
x
W
m
sin
⁡
ω
t
-
t
start
,where f
(
x
) is the same profile function as that used for fringe region; see (6); and W
m is the maximum wall velocity. The parameter ω is the angular frequency of the wall oscillation, which is related to the period through ω
=
2
π
/
T. The parameter t
start determines when the oscillation starts.The sampling time for the reference case was24800 in time units (δ
x
=
0
∗
/
U
∞) and started only after a stationary flow (in the statistical sense) was reached after 13200 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 t
start in (6) is 38000 (=
13200
+
24800).In the simulations presented here, the angular frequency (ω) of the wall oscillation is set to 0.09099 and maximum wall velocity (W
m) is set as 0.5028, which in wall units corresponds to T
+
=
76.4 and W
m
+
=
10.7, calculated as W
m
+
=
W
m
/
u
τ
0 and T
+
=
T
(
u
τ
0
)
2
/
ν, based on u
τ
0 at the start of oscillations (x
=
500), which is equal to 0.047. Note that the value of W
m
+ and T
+ changes with streamwise position because local u
τ
0 varies, whereas W
m and T are kept fixed. At the end of the oscillating part (x
=
1600) the friction velocity of the reference case is u
τ
0
=
0.042 and gives a value of W
m
+
=
12.0 and T
+
=
60.8 at that position. The parameters pertaining to the oscillations are summarized in Table 2.Table 2
Parameters regarding the wall oscillations for the simulation presented.
W
m
+
T
+
x
s
t
a
r
t
x
e
n
d
R
e
Θ
s
t
a
r
t
R
e
Θ
e
n
d
10.7
76.4
500
1600
963
1989


## 3. Results

Before presenting the results in Sections3.1 to 3.3, various definitions needed in the following are presented. The friction velocity is defined as(9)
u
τ
≡
ν
∂
u
∂
y
y
=
0
,where ν is the kinematic viscosity.The resulting drag reduction (DR) is calculated from(10)
D
R
%
=
100
C
f
0
-
C
f
C
f
0
,where C
f
0
=
2
(
u
τ
0
/
u
∞
)
2 is the skin friction of the reference case.An alternative expression of the DR will in this paper be used according to(11)
D
=
C
f
0
-
C
f
C
f
0
.Hence, we may write(12)
D
=
1
-
u
τ
u
τ
0
2
.


### 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,R
e
τ, to R
e
Θ, which is done by using the relation(13)
R
e
τ
=
1.13
R
e
Θ
0.843
,as given by Schlatter and Örlü [32].The data is presented in Figure3 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 W
m
+
=
12 and T
+
=
100. In addition, data from DNS/experiments is indicated by symbols, while the solid and broken lines represent the relations proposed.Figure 3
DR as a function ofR
e
Θ. (–) Profiles according to (14). (– –) Profiles according to (15). Black colour indicates data for oscillations with W
m
+
=
11.3 and T
+
=
67: (◯) DNS from [26], (×) for channel flow from [27], and (+) for boundary layer experiments from [28]. Red colour indicates data for oscillations with W
m
+
=
10 and T
+
=
50: (∗) DNS for channel flow from [12]. Green, blue, and cyan colours indicate data for oscillations with W
m
+
=
12 and T
+
=
100: (⋄) DNS for channel flow from [7], (□) DNS for channel flow from [9], and (×) for channel flow from [27]. Magenta colour indicates data (▽) from the present simulation of a spatially developing boundary layer over oscillations with W
m
+
=
10.7 and T
+
=
76.4.The boundary layer data shown in Figure3 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 W
m
+
=
11.3 and T
+
=
67. 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 R
e
Θ. 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 Figure3) with a curve obtained as(14)
D
R
=
A
R
e
Θ
+
B
,where A
=
2884 and B
=
23.44. This relation suggests that as R
e
Θ 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 atW
m
+
=
12 and T
+
=
100. Touber and Leschziner [7] (diamond) reached R
e
τ
=
1000, while Hurst et al. [9] increased to R
e
τ
=
1600 (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 α
=
-
0.2 and α
=
-
0.26, indicated with the broken lines in Figure 3. These profiles are derived by using the relation between R
e
Θ and R
e
τ above, (see (13)), to calculate the values of β in the expression,(15)
D
R
=
C
R
e
Θ
β
,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 W
m
+
=
12 and T
+
=
100, 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).Table 3
Parameters for the profiles presented in Figure3.
A
/
R
e
Θ
+
B
C
R
e
Θ
β
Line type
Solid
Broken
Colour
W
m
+
T
+
A
B
C
β
Black
11.3
67.0
2884
23.44
Red
10.0
50.0
2679
20.12
Green
12.0
100.0
7143
21.97
140
−0.22
Blue
12.0
100.0
5068
27.66
111
−0.17Choi et al. [12] have provided results for DR in a channel for different R
e
τ with W
m
+
=
10 and T
+
=
50, and the data is presented (stars) in Figure 3. The curve fitting in this case yields A
=
2679 and B
=
20.12 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 Figure3. 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,(16)
C
f
=
0.024
R
e
Θ
-
0.25
,proposed by [34] which indicates vanishing skin friction as the Reynolds number grows. In short, D
R
→
0 because there is no drag to reduce when R
e
Θ
→
∞.


### 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 (Δ
U
+) 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,(17)
D
=
Δ
U
+
2
C
f
0
-
1
/
2
+
1
/
2
κ
,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(18)
D
=
2
C
f
0
1
/
2
Δ
U
+
.The two expressions above were derived in the context of DR by means of riblets, in which caseΔ
U
+ 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 as(19)
Δ
U
+
=
U
∞
u
τ
-
U
∞
u
τ
0
,since 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 Δ
U
+ 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.Figure 4
The variation of DR in the downstream direction. (Black, solid) DR from DNS. (Red, broken) DR from (18). (Green, broken) DR from (17). (Blue, broken) DR from (26). (Blue, solid) DR from (26) with Δ
U
+ from expression (29) (Δ
U
+
=
2.25). (Cyan, solid) DR from expression (30).In the following a new expression for relating DR andΔ
U
+ 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(20)
2
C
f
-
2
C
f
0
=
Δ
U
+
.Multiplying (20) with C
f
/
2 yields(21)
1
-
C
f
C
f
0
=
Δ
U
+
2
C
f
.We wish to replaceC
f with C
f
0 on the right-hand side, so (20) is rewritten as(22)
2
C
f
=
Δ
U
+
+
2
C
f
0
,from which we obtain(23)
C
f
2
=
1
Δ
U
+
+
2
/
C
f
0
,wich, inserted in the right-hand side of (21), gives(24)
1
-
C
f
C
f
0
=
Δ
U
+
Δ
U
+
+
2
/
C
f
0
.Using Taylor series expansion of the left-hand side of (24) and neglecting higher order terms, we obtain(25)
1
-
C
f
C
f
0
=
1
2
1
-
C
f
C
f
0
.Utilizing the approximative relation (25), the following expression for the DR is obtained:(26)
D
=
Δ
U
+
1
/
2
Δ
U
+
+
2
C
f
0
-
1
/
2
.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 Δ
U
+ 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 is(27)
D
=
Δ
U
+
1
/
2
Δ
U
+
+
2
C
f
0
-
1
/
2
+
1
/
2
κ
,which contains an identical correction as compared to expression (17) by a factor Δ
U
+
/
2 in the denominator. However, similar to the original expressions, the inclusion of the growth of the boundary layer via the term 1
/
(
2
κ
) 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,Δ
U
+ needs to be investigated. The downstream variation of the velocity shift according to (19) is shown in Figure 5 as the black curve, and Δ
U
+ 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,(28)
ς
=
y
+
d
u
+
d
y
+
,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,(29)
Δ
U
+
=
k
D
D
,derived from the relation in [11] where the constant was found to be k
D
=
9.0. With the value of D
=
0.25 in the present case, expression (29) yields Δ
U
+
=
2.25. Δ
U
+ 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].Figure 5
The variation of the velocity shiftΔ
U
+ in the downstream direction (x). (Black) Δ
U
+ from (19). (Blue) Δ
U
+ taken from the velocity profiles at the position where the local minimum of (28) is found. (Red) Δ
U
+ from expression (29).Finally, to investigate if the constant value ofΔ
U
+ 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 R
e
Θ through (26).Note that, by using (16) in (26) or (27), we obtain(30)
D
R
~
R
e
Θ
-
0.125for high Reynolds numbers.This variation is shown by the cyan solid line in Figure4 to indicate that this indeed gives identical Reynolds number variation as (26).


## 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,R
e
τ, to R
e
Θ, which is done by using the relation(13)
R
e
τ
=
1.13
R
e
Θ
0.843
,as given by Schlatter and Örlü [32].The data is presented in Figure3 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 W
m
+
=
12 and T
+
=
100. In addition, data from DNS/experiments is indicated by symbols, while the solid and broken lines represent the relations proposed.Figure 3
DR as a function ofR
e
Θ. (–) Profiles according to (14). (– –) Profiles according to (15). Black colour indicates data for oscillations with W
m
+
=
11.3 and T
+
=
67: (◯) DNS from [26], (×) for channel flow from [27], and (+) for boundary layer experiments from [28]. Red colour indicates data for oscillations with W
m
+
=
10 and T
+
=
50: (∗) DNS for channel flow from [12]. Green, blue, and cyan colours indicate data for oscillations with W
m
+
=
12 and T
+
=
100: (⋄) DNS for channel flow from [7], (□) DNS for channel flow from [9], and (×) for channel flow from [27]. Magenta colour indicates data (▽) from the present simulation of a spatially developing boundary layer over oscillations with W
m
+
=
10.7 and T
+
=
76.4.The boundary layer data shown in Figure3 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 W
m
+
=
11.3 and T
+
=
67. 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 R
e
Θ. 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 Figure3) with a curve obtained as(14)
D
R
=
A
R
e
Θ
+
B
,where A
=
2884 and B
=
23.44. This relation suggests that as R
e
Θ 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 atW
m
+
=
12 and T
+
=
100. Touber and Leschziner [7] (diamond) reached R
e
τ
=
1000, while Hurst et al. [9] increased to R
e
τ
=
1600 (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 α
=
-
0.2 and α
=
-
0.26, indicated with the broken lines in Figure 3. These profiles are derived by using the relation between R
e
Θ and R
e
τ above, (see (13)), to calculate the values of β in the expression,(15)
D
R
=
C
R
e
Θ
β
,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 W
m
+
=
12 and T
+
=
100, 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).Table 3
Parameters for the profiles presented in Figure3.
A
/
R
e
Θ
+
B
C
R
e
Θ
β
Line type
Solid
Broken
Colour
W
m
+
T
+
A
B
C
β
Black
11.3
67.0
2884
23.44
Red
10.0
50.0
2679
20.12
Green
12.0
100.0
7143
21.97
140
−0.22
Blue
12.0
100.0
5068
27.66
111
−0.17Choi et al. [12] have provided results for DR in a channel for different R
e
τ with W
m
+
=
10 and T
+
=
50, and the data is presented (stars) in Figure 3. The curve fitting in this case yields A
=
2679 and B
=
20.12 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 Figure3. 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,(16)
C
f
=
0.024
R
e
Θ
-
0.25
,proposed by [34] which indicates vanishing skin friction as the Reynolds number grows. In short, D
R
→
0 because there is no drag to reduce when R
e
Θ
→
∞.


## 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 (Δ
U
+) 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,(17)
D
=
Δ
U
+
2
C
f
0
-
1
/
2
+
1
/
2
κ
,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(18)
D
=
2
C
f
0
1
/
2
Δ
U
+
.The two expressions above were derived in the context of DR by means of riblets, in which caseΔ
U
+ 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 as(19)
Δ
U
+
=
U
∞
u
τ
-
U
∞
u
τ
0
,since 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 Δ
U
+ 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.Figure 4
The variation of DR in the downstream direction. (Black, solid) DR from DNS. (Red, broken) DR from (18). (Green, broken) DR from (17). (Blue, broken) DR from (26). (Blue, solid) DR from (26) with Δ
U
+ from expression (29) (Δ
U
+
=
2.25). (Cyan, solid) DR from expression (30).In the following a new expression for relating DR andΔ
U
+ 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(20)
2
C
f
-
2
C
f
0
=
Δ
U
+
.Multiplying (20) with C
f
/
2 yields(21)
1
-
C
f
C
f
0
=
Δ
U
+
2
C
f
.We wish to replaceC
f with C
f
0 on the right-hand side, so (20) is rewritten as(22)
2
C
f
=
Δ
U
+
+
2
C
f
0
,from which we obtain(23)
C
f
2
=
1
Δ
U
+
+
2
/
C
f
0
,wich, inserted in the right-hand side of (21), gives(24)
1
-
C
f
C
f
0
=
Δ
U
+
Δ
U
+
+
2
/
C
f
0
.Using Taylor series expansion of the left-hand side of (24) and neglecting higher order terms, we obtain(25)
1
-
C
f
C
f
0
=
1
2
1
-
C
f
C
f
0
.Utilizing the approximative relation (25), the following expression for the DR is obtained:(26)
D
=
Δ
U
+
1
/
2
Δ
U
+
+
2
C
f
0
-
1
/
2
.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 Δ
U
+ 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 is(27)
D
=
Δ
U
+
1
/
2
Δ
U
+
+
2
C
f
0
-
1
/
2
+
1
/
2
κ
,which contains an identical correction as compared to expression (17) by a factor Δ
U
+
/
2 in the denominator. However, similar to the original expressions, the inclusion of the growth of the boundary layer via the term 1
/
(
2
κ
) 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,Δ
U
+ needs to be investigated. The downstream variation of the velocity shift according to (19) is shown in Figure 5 as the black curve, and Δ
U
+ 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,(28)
ς
=
y
+
d
u
+
d
y
+
,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,(29)
Δ
U
+
=
k
D
D
,derived from the relation in [11] where the constant was found to be k
D
=
9.0. With the value of D
=
0.25 in the present case, expression (29) yields Δ
U
+
=
2.25. Δ
U
+ 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].Figure 5
The variation of the velocity shiftΔ
U
+ in the downstream direction (x). (Black) Δ
U
+ from (19). (Blue) Δ
U
+ taken from the velocity profiles at the position where the local minimum of (28) is found. (Red) Δ
U
+ from expression (29).Finally, to investigate if the constant value ofΔ
U
+ 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 R
e
Θ through (26).Note that, by using (16) in (26) or (27), we obtain(30)
D
R
~
R
e
Θ
-
0.125for high Reynolds numbers.This variation is shown by the cyan solid line in Figure4 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 D
R
~
R
e
Θ
-
0.125. 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].



---

*Source: 891037-2015-12-14.xml*