The minimal-span channel for rough-wall turbulent flows
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi, R. García-Mayoral
aa r X i v : . [ phy s i c s . f l u - dyn ] M a r Under consideration for publication in J. Fluid Mech. The minimal-span channel for rough-wallturbulent flows
M. MacDonald † , D. Chung , N. Hutchins , L. Chan , A. Ooi andR. Garc´ıa-Mayoral Department of Mechanical Engineering, University of Melbourne, Victoria 3010, Australia Department of Mechanical Engineering, Universiti Tenaga Nasional, Kajang 43000, Malaysia Department of Engineering, University of Cambridge, Cambridge CB2 1PZ, UK(Received ?; revised ?; accepted ?)
Roughness predominantly alters the near-wall region of turbulent flow while the outerlayer remains similar with respect to the wall shear stress. This makes it a prime candi-date for the minimal-span channel, which only captures the near-wall flow by restrictingthe spanwise channel width to be of the order of a few hundred viscous units. Recently,Chung et al. ( J. Fluid Mech. , vol. 773, 2015, pp. 418–431) showed that a minimal-spanchannel can accurately characterise the hydraulic behaviour of roughness. Following this,we aim to investigate the fundamental dynamics of the minimal-span channel frameworkwith an eye towards further improving performance. The streamwise domain length ofthe channel is investigated with the minimum length found to be three times the span-wise width or 1000 viscous units, whichever is longer. The outer layer of the minimalchannel is inherently unphysical and as such alterations to it can be performed so longas the near-wall flow, which is the same as in a full-span channel, remains unchanged.Firstly, a half-height (open) channel with slip wall is shown to reproduce the near-wallbehaviour seen in a standard channel, but with half the number of grid points. Next,a forcing model is introduced into the outer layer of a half-height channel. This re-duces the high streamwise velocity associated with the minimal channel and allows fora larger computational time step. Finally, an investigation is conducted to see if varyingthe roughness Reynolds number with time is a feasible method for obtaining the fullhydraulic behaviour of a rough surface. Currently, multiple steady simulations at fixedroughness Reynolds numbers are needed to obtain this behaviour. The results indicatethat the non-dimensional pressure gradient parameter must be kept below 0 . .
07 toensure that pressure gradient effects do not lead to an inaccurate roughness function. Anempirical costing argument is developed to determine the cost in terms of CPU hours ofminimal-span channel simulations a priori . This argument involves counting the numberof eddy lifespans in the channel, which is then related to the statistical uncertainty of thestreamwise velocity. For a given statistical uncertainty in the roughness function, this canthen be used to determine the simulation run time. Following this, a finite-volume codewith a body-fitted grid is used to determine the roughness function for square-basedpyramids using the above insights. Comparisons to experimental studies for the sameroughness geometry are made and good agreement is observed.
Key words: † Email address for correspondence: [email protected]
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral
1. Introduction
Conventional Direct Numerical Simulations (DNS) of wall-bounded turbulent flowsrepresent a challenging computational problem, as both small and large scales need tobe represented. The former requires a fine grid to resolve the small viscous scales, whilethe latter requires a large domain to capture the large outer-layer motions. However,pioneering work by Jim´enez & Moin (1991) and Hamilton et al. (1995) into minimal-flow units showed that a small computational domain can be used to exclusively capturethe turbulent near-wall cycle, independent of the large outer scales. This is achievedby simply restricting the domain of the channel to a small size, where the spanwiseand streamwise lengths are prescribed in terms of viscous units. Jim´enez & Moin (1991)showed that turbulence could be maintained in the form of the near-wall cycle of thebuffer layer when the spanwise domain width was only on the order of 100 ν/U τ ; here ν is the kinematic viscosity and U τ = p τ w /ρ is the friction velocity, defined using the wallshear stress τ w and the fluid density ρ . This was further supported by future studiesinto minimal-flow units (Jim´enez & Pinelli 1999; Flores & Jim´enez 2010; Hwang 2013;Chung et al. z c . Above this point, the streamwisevelocity increases compared to a full-span channel. This unphysical increase occurs asthe narrow spanwise domain width of the minimal-span channel acts as a filter whichlimits the largest spanwise scale of energy-containing eddies. The flow does, however,retain turbulent scales smaller than the spanwise domain width, so it is not laminar(Jim´enez & Pinelli 1999). The critical value where the minimal-span channel departsfrom the full-span channel has been shown to scale with the spanwise domain width, z c ≈ . L y (Chung et al. z c as ‘healthy’ turbulence, as it is in a full-span channel.In the context of roughness, the central question is how the geometry of a rough surfaceis related to its hydraulic behaviour; namely, what value the (Hama) roughness function,∆ U + , takes (Hama 1954). This quantity reflects the flow retardation, or velocity shift,that the roughness imposes on the flow when compared to a smooth wall, for matchedfriction Reynolds numbers, Re τ = U τ h/ν (here h is the half-channel height, boundarylayer thickness, or pipe radius). For a given surface, ∆ U + can be obtained from varioussemi-empirical models and approximations, or directly from full-scale experiments andnumerical simulations. The former are of varying accuracy and depend on the model se-lected and the rough surface in question, while the latter can be prohibitively expensive.In particular, full-scale numerical simulations suffer from the drawbacks mentioned above,in which both small- and large-scale motions need to be captured. However, roughness isthought to primarily alter only the near-wall flow, in a region called the roughness sub-layer which typically extends 3–5 k from the wall, where k is some characteristic heightof the roughness (Raupach et al. U τ . This is the basis ofTownsend’s outer-layer similarity hypothesis (Townsend 1976) which has received signif-icant attention and has been supported by several rough-wall studies (Flack et al. et al. et al. he minimal-span channel for rough-wall turbulent flows k s , is a single dynamic parameter that is used todescribe the hydraulic behaviour of roughness. It is defined so that the roughness func-tion collapses for all data in the fully rough regime, when the friction factor no longerdepends on the bulk-velocity Reynolds number (Jim´enez 2004, who called this k s ∞ ).However, each rough surface will have a unique behaviour in the transitionally roughregime. It is an expensive process to determine the full hydraulic behaviour, as a rangeof simulations need to be conducted for the different roughness Reynolds numbers, eachrequiring a different body-fitted grid and its own initialisation period. It would there-fore be desirable to conduct a simulation in which only a single computational grid isused, and the bulk velocity is changed over time to sweep through a range of roughnessReynolds numbers. This is a similar approach to how experimental studies are performedin that a single rough surface with fixed k/h is tested at multiple flow speeds, so thatthe roughness Reynolds number varies. In order to obtain statistics at a desired frictionReynolds number in a temporally evolving flow, statistics are averaged over a small win-dow in which the instantaneous friction Reynolds number is close to the desired one.Within this window, the flow is assumed to be quasi steady. This would therefore gen-erate a near-continuous curve of ∆ U + versus k + , rather than the conventional (steady)approach which only generates a few data points. An important consideration which willbe investigated here is whether acceleration effects from the changing mass flux becomesignificant and distort the estimation of ∆ U + , that is, how quickly can the bulk velocitybe varied such that the quasi-steady assumption remains approximately valid.Recently, we have applied the idea of minimal-span channels to the roughness problemin a proof-of-concept study (Chung et al. § § § § § a priori . The insights gained in this paper arethen used in §
2. Numerical Method
The majority of Direct Numerical Simulations in this paper are conducted using thesame finite difference code as in Chung et al. (2015), which uses a fully conservativefourth-order staggered-grid scheme. Time integration is performed using the third-order
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral low-storage Runge–Kutta scheme of Spalart et al. (1991), and the fractional-step methodof Kim & Moin (1985) is used. The Navier–Stokes equations are ∇ · u = 0 , ∂ u ∂t + ∇ · ( uu ) = − ρ ∇ p + ∇ · ( ν ∇ u ) + F + G + K , (2.1)where u = ( u, v, w ) is the velocity in the streamwise ( x ), spanwise ( y ) and wall-normal ( z )directions, t is time, and p is pressure. G = G x ( t ) i is the spatially invariant, time-varyingstreamwise forcing term which drives the flow at constant mass flux. The flow is solvedin a reference frame in which the mass flux is zero at all times, although all equations inthis paper are given in the stationary-wall reference frame. Solving in a zero-mass-fluxreference frame permits a larger computational time step (Lundbladh et al. et al. F = − αF ( z, k ) u | u | i , is basedon the work of Busse & Sandham (2012), which applies a forcing in the streamwisedirection that opposes the flow. The roughness factor α = 1 / (40 k ) is kept constant forall simulations, with the roughness height k = h/
40. The function F ( z, k ) is a simplestep function which applies the roughness forcing adjacent to the top and bottom no-slipchannel walls at z = 0 and z = 2 h , F ( z, k ) = ( , if z < k or 2 h − k < z. , otherwise . (2.2)This roughness model removes any spanwise or streamwise roughness length scales fromthe problem and allows for simpler computations as the same smooth-wall grid can beused. This is a similar model to Borrell (2015), although there the author used an F ∝ − u i scaling. The final term in (2.1), K , is a forcing function which damps the velocity in theouter layer of the minimal channel. It has the form K i = − γ Γ( z, z d ) (cid:16) u i ( x, y, z, t ) − h u i i ( z d , t ) (cid:17)(cid:12)(cid:12)(cid:12) u i ( x, y, z, t ) − h u i i ( z d , t ) (cid:12)(cid:12)(cid:12) , (2.3)(no summation over i ) where angled brackets denote the spatial average of the instan-taneous velocity over a wall-parallel plane. The factor γ has units of inverse time andin this study is set according to γh/U τ ≈
1. The parameter Γ( z, z d ) is similar to thefunction F ( z, k ) of the roughness forcing model in that it indicates where the dampingis applied, namely in the outer layer of the channel,Γ( z, z d ) = ( , if z > z d and z < h − z d . , otherwise . (2.4)The value of z d should be greater than the critical wall-normal location, z c , of theminimal-span channel, so that the forcing does not contaminate the healthy turbulenceof the near-wall flow. A step function for Γ is used to minimise the parameter space,as well as to ensure the location of z d is unambiguous. Several different values of z d aretested and will be presented in the following section. In (2.3), the term h u i i ( z d , t ) is thewall-parallel spatially averaged velocity at z d . This is present to reduce the magnitude ofthe streamwise velocity; the minimal channel has a very high streamwise velocity whichpresents a time step restriction from the CFL number. The forcing is in this form so thatthe velocity is forced to remain at the same level as at z = z d . This outer-layer dampingis similar to the masks employed in Jim´enez & Pinelli (1999), which damp fluctuationsin the outer-layer region of minimal channels. The current forcing model simply sets thestreamwise velocity to be approximately that at z = z d , however it will be shown to workreasonably well for the purpose of obtaining ∆ U + . Note that the average velocity in the he minimal-span channel for rough-wall turbulent flows L x L y L z = h z = 0, no-slip wall z = kz = z d z = h , slip wallFlowOuter-layer damping ❍❍❍❍ Roughness forcing ❍❍❍
Figure 1: (Colour online) Half-height minimal-span channel, showing the roughness forc-ing zone (2.2) and the outer-layer damping zone (2.3). No-slip wall at z = 0 and free-slipwall at z = h .spanwise and wall-normal directions is zero, so that the forcing in these directions cansimply be K i = − γ Γ( z, z d ) u i | u i | , for i = 2 ,
3. Figure 1 shows the two forcing regions thatare employed in the current study.The streamwise and spanwise grid is evenly spaced, while the wall-normal grid isstretched with a cosine mapping. Periodic boundary conditions are applied in the span-wise and streamwise directions. In the case of standard-height channels, no-slip walls arelocated at z = 0 and z = 2 h . For clarity, the word ‘full’ will be used to refer to thespan (full span), while the case with two no-slip walls will be referred to as ‘standard’height. For half-height (open) channels, a no-slip wall is still positioned at z = 0, howevernow the top domain surface is a free-slip wall with ∂u/∂z = ∂v/∂z = 0, positioned at z = h . This slip wall maintains the impermeability constraint, w = 0. The spanwisewidth, L y , of the minimal-span channel satisfies the guidelines of Chung et al. (2015),namely that L y & max(100 ν/U τ , k/ . , λ r,y ), where λ r,y is the spanwise length scale ofroughness elements. Because the homogeneous roughness forcing model has no spanwiselength scale then the final constraint can be ignored. Simulations are conducted at afriction Reynolds number of Re τ ≈ Re τ ≈ Temporal sweep
As discussed in the introduction, instead of conducting multiple steady simulations ofdifferent roughness Reynolds numbers in order to determine the equivalent sandgrainroughness, a single unsteady simulation is to be performed in which the bulk velocity isvaried. In particular, we will investigate different rates of change of the bulk velocity andthe effects that this acceleration has on the flow. The sweep will start with the highestfriction number to be tested, Re τ,start and will then be decelerated via an adverse pressuregradient. This ensures an adequate grid resolution at the start of simulation. At the finalfriction Reynolds number, Re τ,end < Re τ,start , the grid resolution would be finer thannecessary, which for the current simulations would have 4–5 times more cells than if aconventional steady simulation were conducted at this final Reynolds number.Presently, we vary the bulk velocity U b ( t ) with time,d U b d t = A = U b,end − U b,start ∆ T = const., (2.5)where U b,end is the desired end point of the simulation (the bulk velocity corresponding M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral tU τ /h R e τ tU τ /h ∆ p . . . . . . ( a ) A A A ( b ) ∆ p,Re τ =180 ≈ . ∆ p,Re τ =180 ≈ . ∆ p,Re τ =180 ≈ . Figure 2: ( a ) Instantaneous Reynolds number and ( b ) pressure gradient parameter (2.6)of the smooth-wall temporal sweep, against viscous time. The cases with faster sweepshave multiple runs (shown by grey lines) and are ensemble averaged. A is the linear rateof change of U b (2.5).to Re τ = 180), and ∆ T is a time scale which should be sufficiently long enough thatacceleration effects are not significant to the flow. The initial bulk velocity, U b,start , isset corresponding to a friction Reynolds number of Re τ,start = 590 and different valuesof A are tested. An initial value of A is selected such that over the course of the sweep∆ T U τ /h = 30 ⇒ U τ /hA ≈
30, with subsequent runs using 2 A and 4 A . In these casesof a higher rate of change of U b , multiple sweeps are run from different initial conditionswith the results ensemble averaged, to ensure that statistics are obtained over the sameamount of simulation time. For example, the case where the gradient is quadrupled, foursweeps are conducted with four different initial conditions. This is visualised in figure2( a ), which shows the change in the friction Reynolds number as a function of time.An important consideration with this technique is that the friction velocity and hencethe normalised spanwise domain width L + y will vary with time. The domain width mustremain larger than 100 viscous units at all times to ensure the turbulent flow is sustained,which may necessitate multiple stages in the sweep if L + y becomes too small. Each stagein the sweep should be set up such that L + y & ∆ p = νρU τ d p d x . (2.6)Various experimental studies of decelerating boundary layers show that mild deceler-ations of ∆ p < .
01 have small ( < U τ from methods basedon the assumption of zero pressure gradients (Patel 1965; Jones et al. et al. (2014), meanwhile, compared a step acceleration with a slow ramp up inwhich the bulk velocity was linearly varied in a channel such that − ∆ p ≈ .
73 (favourablepressure gradient). Acceleration effects were still seen in this slow ramp up. Note thatfor steady channel flow, the pressure gradient parameter is ∆ p = − /Re τ <
0, whichfor the current simulations at Re τ = 590 takes a value of − . b ), the linear variation in U b results in the pressure gradient parameter having a maximum value at the end ofthe simulation, when Re τ ≈ he minimal-span channel for rough-wall turbulent flows ID Re τ L x h L + x L + y L z h N x N y N z ∆ x + ∆ y + ∆ z + w ∆ z + h ∆ t + S ∆ t + R Full-span channelFS 590 2 π . π
190 118 2 24 24 256 7.7 4.9 0.04 7.2 0.28 0.32MS2 590 0 . π
300 118 2 48 24 256 6.2 4.9 0.04 7.2 0.24 0.27MS3 590 0 . π
460 118 2 48 24 256 9.7 4.9 0.04 7.2 0.36 0.42MS4 590 0 . π
930 118 2 96 24 256 9.7 4.9 0.04 7.2 0.36 0.42MS5 590 π . π
190 354 2 24 72 256 7.7 4.9 0.04 7.2 0.37 0.45MSW2 590 0 . π
300 354 2 48 72 256 6.2 4.9 0.04 7.2 0.31 0.40MSW3 590 0 . π
460 354 2 48 72 256 9.7 4.9 0.04 7.2 0.46 0.52MSW4 590 0 . π
930 354 2 96 72 256 9.7 4.9 0.04 7.2 0.43 0.46MSW5 590 π Table 1: Description of the simulations performed with the finite difference code to inves-tigate streamwise domain length. All simulations done as both smooth wall and modelledrough wall. Entries: Re τ , nominal friction Reynolds numbers; L , domain length; N , num-ber of grid points; ∆, grid-spacing in viscous units, subscript w and h refers to the wall-normal grid spacing at the wall and channel centre; L z /h = 2 denotes standard-height(two no-slip walls), ∆ t + S and ∆ t + R is the smooth- and rough-wall time step. Roughnessheight k = h/
40 so k + ≈
15 at Re τ = 590.their maximum pressure gradient parameter, which are ∆ p,Re τ =180 ≈ (0 . , . , . ∆ p < .
01) and noticeable ( − ∆ p > .
73) pressure gradient effects, and are substantiallylarger than the steady value ( − ∆ p = 0 .
3. Results
Varying streamwise domain length (table 1)
Wall-bounded turbulent flows are often characterised by very long structures on the orderof tens of channel half heights (Kim & Adrian 1999; Lozano-Dur´an & Jim´enez 2014).This becomes very expensive to simulate for conventional full-span channel simulations,so shorter domain lengths of L x /h = 2 π are often used (Chin et al. et al. h .The minimal-span channel only captures the inner-layer flow which is not dependent on h , suggesting a shorter domain length can be used. The near-wall cycle of the buffer layerproduces streaky structures with streamwise lengths of 1000 viscous units (Kline et al. k x = 0) mode and hence the domainlength doesn’t have to be this long. These streaks are accompanied by quasi-streamwisevortices with streamwise lengths of 200–300 viscous units (Jeong et al. M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral z + U + z + U + ( a ) ( b )Smooth, L + y ≈ L + x ≈ , ❍❍ L + x > ❅❅ Rough, L + y ≈ z + U + S − U + R L + x ≈ , ❈❈ z + U + z + U + ( c ) ( d )Smooth, L + y ≈ L + x ≈ ❅❅ Rough, L + y ≈ L + x ≈ ❅❅ z + U + S − U + R Figure 3: Streamwise velocity profile for varying streamwise length for ( a , c ) smooth- and( b , d ) rough-wall standard-height channels. Spanwise domain width is ( a , c ) L + y ≈ b , d ) L + y ≈ k = h/
40; vertical dotted line indicates critical height z c = 0 . L y . Insets of ( b , d ) shows difference in smooth and rough wall flows.rowest minimal-span channels with L + y ≈
100 therefore require L + x ≈ et al. et al. λ y = L y and λ x = 2–3 L y .Two different spanwise domain widths are investigated for varying the streamwiselength (table 1). Firstly, the smallest width of L + y ≈
120 would have the largest eddieshaving streamwise lengths of λ + x = 240–350. Various streamwise domain lengths of L + x =(190 , , , , L + y ≈
350 would be able to capture larger eddies, the largestof which would have a streamwise length of approximately λ + x = 710–1100. The sameset of streamwise domain lengths are simulated for this wider domain. A standard-heightchannel is used, with no outer-layer damping, K i = 0.Figure 3 shows the mean velocity profile for the two spanwise widths, for all thestreamwise lengths tested. The minimal channels (grey lines) tend to agree with the full- he minimal-span channel for rough-wall turbulent flows z c = 0 . L y (verticaldotted line), at which point the velocity profile of the minimal channels increases. Forthe smaller domain width (figure 3 a , b ), the shortest streamwise lengths of L + x ≈
190 and300 result in a reduction in the smooth-wall velocity above z + &
20. There is also a slightincrease in the centreline velocity. As a result, the difference in smooth- and rough-wallvelocity (inset of figure 3 b ) decreases relative to the full-span and longer minimal-spancases. This effect is not due to the differences in grid resolution of the two shortestchannels (table 1), as the wider channel has the same differences in grid resolution butwill be seen to not have this effect. Little discernible difference can be seen between thelonger streamwise lengths of L + x > k x = 0 mode.The case of L + x ≈
300 = 2 . L + y produces an incorrect profile and this suggests that aminimum streamwise length of L + x = 3 L + y is required, especially for narrow channels inwhich z c is closer to the buffer layer than log layer of the flow.The smallest streamwise length of L + x ≈
190 was unable to sustain turbulence for aprolonged period. It would decay to a laminar state after T + = T U τ /ν ≈ . × inthe smooth-wall channel, and after T + ≈ . × in the rough-wall channel. The datashown in the previous figures for this streamwise length is averaged only when the flowis turbulent. This behaviour is similar to what was observed in Jim´enez & Moin (1991),who showed that turbulence could not be maintained in channels with streamwise lengthsof around 200 viscous units. As discussed above, this is because it is shorter than thequasi-streamwise vortices which have streamwise lengths of λ + x ≈ L + y ≈
354 (figure 3 c , d ) results in a centrelinevelocity that is closer to the full-span channel, as these wider channels capture a widerrange of turbulent motions. However, the case with the shortest streamwise length of L + x ≈ ≈ . L + y has a larger centreline velocity, with the critical wall-normal heightappearing lower than z c = 0 . L y . This somewhat resembles the channels with a narrowerspanwise width of L + y ≈
120 (figure 3 a , b ), suggesting that the very short domain lengthrestricts the size of the largest eddies. Their spanwise width would now be smaller thanthe width of the domain, i.e. the streamwise length is now the limiting length scale.Interestingly, the streamwise length of L + x ≈
300 = 0 . L + y appears to agree quite wellwith the cases with longer lengths. This is in contrast to the narrow domain (figure3 a ), which shows a clear difference for L + x ≈ ≈ . L + y . Even the case with L + x ≈ ≈ . L + y is producing a velocity profile and roughness function comparable to thatof channels with longer streamwise lengths, despite not having the requisite scaling of L + x & L + y discussed above. This is similar to the results of Toh & Itano (2005) wholooked at wide spanwise channels with very short streamwise lengths. The velocity profilesthey obtained look similar to a full-span, full-length channel, with no apparent increasein streamwise velocity in the outer layer that is characteristic of minimal-span channels.A possible explanation for the results seen here is that the wall-normal critical location z + c ≈
140 = 0 . h + is outside the log layer, which is generally believed to end at z ≈ . h (Marusic et al. L y is no longerappropriate in the outer layer. In this case, the largest eddy at the edge of the logarithmiclayer would have a streamwise length of approximately λ + x ≈ L + x ≈
460 produces a reasonable velocity profile.Premultiplied one-dimensional streamwise energy spectra of streamwise velocity, k x E + uu ,are shown in figure 4, where u ′ rms = R ∞ k x E uu d log( k x ). Only the smooth-wall flow isshown as the rough-wall data exhibits the same qualitative trends as the smooth walland so is omitted. This occurs because the rough-wall flow is in the transitionally rough0 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral λ + x z + λ + x z + Increasing width L y ✲ ( a ) L + y ≈
120 ( b ) L + y ≈ λ + x = 3 L + y λ + x = 3 L + y λ + x z + λ + x z + I n c r e a s i n g l e n g t h L x ( c ) L + y ≈
120 ( d ) L + y ≈ λ + x = 3 L + y λ + x = 3 L + y λ + x z + λ + x z + ❄ ( e ) L + y ≈
120 ( f ) L + y ≈ λ + x = 3 L + y λ + x = 3 L + y λ + x z + λ + x z + ( g ) L + y ≈
120 ( h ) L + y ≈ λ + x = 3 L + y λ + x = 3 L + y Figure 4: Smooth-wall premultiplied one-dimensional streamwise energy spectra ofstreamwise velocity, k x E + uu , for varying streamwise lengths of ( a , b ) L + x ≈ c , d ) L + x ≈ e , f ) L + x ≈ g , h ) L + x ≈ a , c , e , g ) L + y ≈ b , d , f , h ) L + y ≈ k x E + uu = 0 .
5, 1.0, 1.5, 2.0. Inset box gives channelscale, arrow shows direction of streamwise flow. Vertical dashed line shows λ + x = 3 L + y ,horizontal dashed line shows z + c = 0 . L + y . Square dashed box of wider channel ( b , d , f , h )shows the additional captured region. he minimal-span channel for rough-wall turbulent flows λ + y z + λ + y z + Increasing width L y ✲ ( a ) L + x ≈
190 ( b ) L + x ≈ λ + y = L + x / λ + y = L + x / λ + y z + λ + y z + I n c r e a s i n g l e n g t h L x ( c ) L + x ≈
460 ( d ) L + x ≈ λ + y = L + x / λ + y = L + x / λ + y z + λ + y z + ❄ ( e ) L + x ≈
930 ( f ) L + x ≈ λ + y = L + x / λ + y = L + x / λ + y z + λ + y z + ( g ) L + x ≈ h ) L + x ≈ λ + y = L + x / λ + y = L + x / Figure 5: Smooth-wall premultiplied one-dimensional spanwise energy spectra of stream-wise velocity, k y E + uu , for varying streamwise lengths of ( a , b ) L + x ≈ c , d ) L + x ≈ e , f ) L + x ≈ g , h ) L + x ≈ a , c , e , g ) L + y ≈ b , d , f , h ) L + y ≈ k y E + uu = 1 .
0, 2.0, 3.0, 4.0. Inset box gives channel scale, arrowshows direction of streamwise flow. Vertical dashed line shows λ + y = L + x /
3, horizontaldashed line shows z + c = 0 . L + y .2 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral regime, which leads to a weakening of the near-wall cycle (MacDonald et al. λ x ≈ λ y would likely still hold as this is generated by a self-sustaining process in the logarithmiclayer (Flores & Jim´enez 2010; Hwang 2015). The vertical dashed lines correspond to thelongest streamwise scale based on the log-layer scaling of 3 L + y discussed above, whilethe termination of the line contours of the minimal channels show the channel stream-wise length. The dashed square in the wider channels (figure 4 b , d , f , h ) shows the extracaptured region due to having a larger L + y , compared to the narrower minimal channels(figure 4 a , c , e , g ).It is clear that for the narrowest channel with L + y ≈ L + x ≈
190 = 1 . L + y (figure 4 a ) is too short. This is seen to result in an increasein turbulent energy above z + c (horizontal dashed line), relative to the full-span channel(shaded contour). A similar effect is observed for L + x ≈
300 = 2 . L + y (not shown). Thisincrease in energy at smaller streamwise length scales above z + c is in agreement with pre-vious minimal-channel simulations (Hwang 2013). When the streamwise domain lengthexceeds the 3 L y scaling with L + x ≈
460 = 3 . L + y (figure 4 c ) then reasonable agreementwith the full-span channel is observed, despite the narrow range of scales captured bythe minimal channel. Further increases to the streamwise length for this domain width(figure 4 e , g ) do not improve the collapse with the full-span channel, especially near z + c ,as no additional turbulent motions are captured according to the 3 L + y scaling discussedabove. The increased length is however able to capture the near-wall cycle with peak at z + ≈
15 and λ + x ≈ g ).A similar picture emerges for the wider minimal channel with L + y ≈ L + x ≈ ≈ . L + y (figure 4 b ), is unable to capture much of theturbulent motions, similar to that observed in the narrower minimal channel of figure4( a ). A streamwise domain length of L + x ≈ ≈ . L + y (figure 4 d ) has better agreementwith the full-span channel, although there is still some discrepancy. The 3 L y scaling isalmost reached in figure 4( f ) with L x + ≈ ≈ . L y , and excellent agreement isobserved with the full-span channel, with only a slight increase in energy above z + c .Further increasing L + x above 3 L y (figure 4 h ) does not provide any improvement otherthan capturing more of the near-wall cycle, further supporting this scaling argument.The premultiplied spanwise energy spectra of streamwise velocity are shown in figure 5.Here, the vertical dashed line now shows the widest spanwise length scale that can existbased on the streamwise domain length using the λ y = 3 L x scaling. For the narrowerchannel with L + y ≈ L + x ≈
460 = 3 . L + y (figure 5 c ), andfurther increases to L + x (figure 5 e , g ) do not result in any change to the turbulent energy.Similarly for the wider channel with L + y ≈ f ) with L + x ≈
930 = 2 . L + y and increasing L + x to 1850 viscous units (figure 5 h )does not result in any substantial change to the spectra. The cases with L x . L y (figure 5 a , b , d ) result in enhanced turbulence energy, particularly in the near-wall peakat z + ≈
15 and λ + y ≈ L x = 3 L y .This is due to quasi-streamwise vortices that exist in the logarithmic layer (where neitherviscosity or the channel half height are relevant length scales), suggesting this scaling isindependent of Re τ . Smaller streamwise lengths than this scaling lead to poor agreementbetween the minimal and full-span channel, as seen in figures 4 and 5. However, for verynarrow channels this can result in a streamwise domain length less than 1000 viscous he minimal-span channel for rough-wall turbulent flows ID Re τ L x h L + x L + y L z h N x N y N z ∆ x + ∆ y + ∆ z + w ∆ z + h ∆ t + S ∆ t + R Full-span channelFS 590 2 π π π π Table 2: Full-span and minimal-span channel simulations, for standard-height (two no-slip walls, L z /h = 2) and half-height (one no-slip and one slip wall, L z /h = 1) channels.Entries are same as table 1. z + U + z + U + z + U + S − U + R ( a ) ( b )Full-span z + U + S − U + R Minimal-span
Figure 6: Mean streamwise velocity against wall-normal position for ( a ) full-span channeland ( b ) minimal-span channel. Line styles are: black, standard-height channel; grey, half-height channel; solid, smooth wall; dashed, rough-wall model. Vertical dashed line showsroughness crest, k = h/
40, vertical dotted line shows critical height, z c = 0 . L y . Insetshows difference in smooth- and rough-wall velocity, U + S − U + R .units. While it appears the buffer-layer streaks do not need to be captured in the domain,having such a small streamwise domain exacerbates the bursty nature of the minimal-spanchannel (Jim´enez 2015). Only one or maybe two quasi-streamwise vortices are presentin the domain, which makes obtaining converged statistics difficult as the simulationneeds to be run for a significantly long time, an issue discussed in detail in §
4. Webelieve that a minimum streamwise length of approximately 1000 viscous units seems areasonable length in these cases, so that several of the smallest quasi-streamwise vorticesare present. As such, we recommend that L x & max( 3 L y , ν/U τ , λ r,x ), where thelast requirement pertains to the streamwise length scale, λ r,x , of the roughness which isabsent in the current roughness forcing model.3.2. Half-height channel (table 2)
The previous section looked at the effect of the streamwise domain. Now we will considerthe effect of the wall-normal domain, particularly in terms of the outer-layer flow whichis unphysical in minimal channels. First, we will consider a half-height (open) channelwhich consists of a slip wall positioned at the channel centreline. Intuitively, changing4
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral z + u ′ + i , r m s z + u ′ + i , r m s ( a ) ( b )Full-span u ′ + x,rms u ′ + y,rms u ′ + z,rms Minimal-span u ′ + x,rms u ′ + y,rms u ′ + z,rms Figure 7: Streamwise, spanwise, and wall-normal root-mean-squared velocity fluctuationsagainst wall-normal position for ( a ) full-span channel and ( b ) minimal-span channel. Linestyles are same as figure 6.the boundary condition at z = h is unlikely to effect the flow at z = z c , given thedistance between these two heights ( z c /h ≈ . L + x ≈ K i = 0.Figure 6 shows the effect of the half-height channel in terms of the mean velocity profile,when compared to the standard-height channel. This is done for both full-span (figure6 a ) and minimal-span (figure 6 b ) channels, for smooth and rough walls. For clarity, thefull-span and minimal-span channels are shown in different figures, however the near-wall behaviour is identical up until z c = 0 . L y ⇒ z + c ≈
50, as observed previously (figure3). Both sets of figures show that the use of the half-height channel with slip wall hasa negligible effect on the flow. The main difference is seen in the wake region, wherethe half-height channel restricts the outer-layer structures, resulting in a slight decreasein the mean velocity. Crucially, this difference is the same for both smooth and roughwalls, meaning that the difference between them, U + S − U + R , is the same. Moreover, thechange in the half-height channel occurs above the critical wall-normal location, wherethe minimal-span flow is already altered compared to the full-span flow.The streamwise, spanwise, and wall-normal root-mean-squared velocity fluctuationsare shown in figure 7 for both full-span and minimal-span channels, comparing standard-height and half-height channels. For the full-span channel (figure 7 a ), these velocityfluctuations show very good agreement between the standard-height (black lines) andhalf-height (grey lines) channels in the near-wall region, in agreement with previousfull-span half-height channel studies (Handler et al. z + &
40 relative to the standard-heightchannels, however the difference is relatively small. Moreover, outer-layer similarity isstill maintained between the smooth- and rough-wall channels, suggesting the effect onthe difference between the two flows will be minor.For the minimal-span channel (figure 7 b ), a slightly different picture emerges. Firstly,it should be noted that a standard-height minimal-span channel (black lines of figure 7 b ) he minimal-span channel for rough-wall turbulent flows λ + x z + z + p ′ + r m s ( a ) ( b ) Figure 8: ( a ) Smooth-wall premultiplied one-dimensional streamwise energy spectrum ofpressure, k x E + pp . Shaded contour is full-span channel, line contour is minimal-span half-height channel. Horizontal dashed line shows z + c = 0 . L + y , vertical dashed line shows 3 L + y ( § b )Pressure root-mean-squared fluctuations against wall-normal position. Line styles: black,full-span channel; grey solid, minimal-span half-height channel; grey dashed, minimal-span channel with altered velocity field u a (see text). Vertical dashed line shows z + c .have enhanced streamwise and wall-normal velocity fluctuations compared to its full-spancounterpart (black lines of figure 7 a ) in the outer layer. This has been noted in otherminimal-span channel studies (e.g. Hwang 2013; MacDonald et al. b ) are damped from z + &
20 compared to the standard-height minimal-span channel. This is the opposite ofwhat occurred in a full-span channel. Interestingly, Hwang (2013) found that when thespanwise uniform eddies (that is, the k y = 0 mode) were filtered out of a minimal-spanchannel, a similar behaviour was observed in that the wall-normal and streamwise velocityfluctuations were damped compared to the unfiltered channel. It is unclear whether thehalf-height channel is performing a similar operation to this k y = 0 filtering, as thiswould imply the impermeability constraint is limiting the infinite spanwise motions insome way. However, it is clear that the near-wall flow is preserved, despite the impositionof an outer-layer boundary condition. This shows that a half-height channel flow can besimulated without significant near-wall detriment.The pressure fluctuations, meanwhile, show a substantial increase in the minimal-span channel compared to the full-span channel. The premultiplied streamwise spectraof pressure is shown in figure 8( a ), and it is clear that the fluctuations are larger at longerwavelengths ( λ + x ≈ z + c . The fluctuations at wavelengthsshorter than 3 L + y and below z + c are, however, in reasonable agreement with the full-span channel. Despite this, because the root-mean-squared pressure fluctuations are theintegral of this spectrum ( p ′ rms = R ∞ k x E pp d log( k x )), then the profile is nearly an orderof magnitude larger in the minimal-span channel (grey line, figure 8 b ). Note that themean pressure is still correct, as this must obey the averaged wall-normal momentumequation, w + p = C , where C is a constant.A possible solution to this behaviour is to consider decomposing the pressure into itsrapid and slow parts (Kim 1989). Here, the rapid pressure emerges due to the meanstreamwise velocity gradient while the slow pressure is due to the velocity fluctuationgradients. The mean velocity gradient dU/dz in the minimal-span channel is enhanced6 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral ID Re τ L x h L + x L + y L z h N x N y N z ∆ x + ∆ y + ∆ z + w ∆ z + h z + d ∆ t + S ∆ t + R Minimal-span, half-height channel with outer-layer dampingMSHD1 590 2 π π π π π π π π π Re τ )MSHRD1 2000 0 . π . π . π Table 3: Minimal-span simulations with outer-layer damping. Refer to table 1 for def-initions. z + d indicates location where damping starts (figure 1), cases with no dampingindicated by a hyphen.above z c , which means the rapid pressure will also be larger. If we instead define analtered streamwise velocity field whereby u a ( x, y, z ) = u ( x, y, z ) − ( U ( z ) − U ( z c )) for z > z c where U ( z c ) is the mean streamwise velocity at z = z c , this will then set the mean dU a /dz = 0. If the other two velocity components are unchanged and u a = u for z z c ,then the rapid pressure would be zero above z c . The pressure field arising from thisaltered velocity field, ∇ p a = − u ai,j u aj,i can then be computed and is seen to be in muchbetter agreement with the full-span channel (figure 8 b ). There is a slight difference in thenear-wall peak near z + ≈
30, however this may be due to the discontinuity in dU a /dz at z = z c . This correction to pressure is an additional step that needs to be performedwhenever pressure statistics are to be outputted. The unaltered (discrete) pressure is stillnecessary at each time step to ensure the flow remains divergence free.3.3. Outer-layer damping (table 3)
The previous section showed that a half-height (open) channel did not significantly alterthe healthy near-wall flow. A more aggressive alteration is now attempted in which theouter-layer flow is damped through use of the K i forcing term (2.3) in a half-height chan-nel. This will reduce the large centreline velocity (figure 6 b ) of the minimal-span channel,which places a restriction on the time step due to the CFL number. The streamwise do-main length is fixed at L + x ≈ z d , is varied in channels of two different spanwisewidths; L + y ≈
120 and L + y ≈ z + c ≈
47 and z + c ≈
94, respectively.Figure 9 shows the mean streamwise velocity profile for smooth- and rough-wall minimal-span channels, for both of the spanwise widths at Re τ = 590. The smallest spanwisewidth, L + y ≈ ⇒ z + c ≈
47 (figure 9 a , b ) has four different positions for z d examined,with values of z + d = (100 , , , z > z d to the value of U ( z = z d ). It is clear in the smooth-wall he minimal-span channel for rough-wall turbulent flows z + U + z + U + ( a ) ( b )Smooth, L + y ≈ z d ❆❆❯ ✁✁✁✕ Rough, L + y ≈ ✓✓✓✓✼ z + U + S − U + R ❆❆❯ z + U + z + U + ( c ) ( d )Smooth, L + y ≈ (cid:0)(cid:0)✒ Rough, L + y ≈ (cid:0)(cid:0)✒ z + U + S − U + R Figure 9: Streamwise velocity profile for outer-layer damping for ( a , c ) smooth- and( b , d ) rough-wall half-height channels. Spanwise domain width is ( a , b ) L + y ≈
120 and( c , d ) L + y ≈
240 (table 3). Line styles are: black, full-span channel; grey, minimal-span channel with outer-layer damping (2.3) for varying damping heights of ( a , b ) z + d = (100 , , , c , d ) z + d = (200 , , +++ ; light-grey, minimal-span with no damping. Vertical dashed line indicates the roughness height k = h/
40. Insetof ( b , d ) shows the difference in smooth- and rough-wall velocity. Arrows indicates trendof increasing z d .flow (figure 9 a ) that the two smallest values of z + d = 100 and 150 are too close to thewall, contaminating the healthy near-wall flow and resulting in an increase in the stream-wise velocity below z d . In particular, z + d = 100 results in an obvious overestimation of∆ U + , as seen in the inset of figure 9( b ). A value of z + d = 150 has a similar, although notas strong effect. This is clearly not desirable and suggests that values of z + d >
200 arenecessary to avoid contamination of the healthy near-wall flow.To investigate how the location of z d is related to z c will require cases with a largerspanwise domain width, as z c ≈ . L y . This is done in figure 9( c , d ), where the spanwisewidth is now L + y ≈ ⇒ z + c ≈
94. Here, three different positions for z d are analysed,with values of z + d = (200 , , z d positions. This suggests that, even though z c has doubled, westill require the damping to begin at least 200 viscous-units away from the wall. Havingthe outer-layer damping begin closer to the wall (i.e. reducing z d ) leads to its effectspermeating the entire near-wall region, increasing the viscous stress in the log-layer andeven buffer layer. This results in an overestimation of the mean streamwise velocity. Thestrong effect of the damping when z d is small could also be due to the strength of the8 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral z + U + z + U + ( a ) ( b )Smooth Inc. z d (cid:0)(cid:0)✒ Rough (cid:0)(cid:0)✒ z + U + S − U + R Figure 10: Streamwise velocity profile for ( a ) smooth- and ( b ) rough-wall half-heightchannels at Re τ = 2000 (table 3). Line styles are same as figure 9, however full-spansmooth-wall channel data taken from Hoyas & Jim´enez (2006) and damping heights are z + d = (200 , ××× .damping, γh/U τ ≈ Re τ = 2000. Since the roughnessheight is fixed as a ratio of the channel half-height, k = h/
40, then the roughness Reynoldsnumber increases to k + = 50. The spanwise width also needs to be increased to submergethe roughness sublayer in healthy turbulence (Chung et al. L + y = 300 so that the critical wall-normal location is z + c ≈
120 = 2 . k . Figure 10 showsthe mean velocity profile for these higher Re τ simulations in which z + d = (200 , z + ≈
120 = 0 . L + y (figure 10 a ). Here, both positions of z d produce velocity profiles which compare well withthe minimal-span channel with no outer-layer damping, up to the forcing location z d .It is also of interest to see what effect the outer-layer damping has on second-orderturbulence statistics. This is shown in figure 11 for the streamwise and wall-normalvelocity fluctuations (figure 11 a ) and the Reynolds shear stress (figure 11 b ), for the caseof outer forcing with z + d = 200 and L + y = 120. This is compared to the minimal-spanchannel with the same spanwise width, but with no outer-layer damping. The turbulenceintensity up until z + ≈
30 agrees reasonably well between the case with forcing andwithout, however because the forcing term scales as ( u i − h u i i ) , fluctuations are almostzero above z > z d . The fact that they are not exactly zero is likely because only thespatially averaged velocity h u i i is used in the forcing term, and so will vary slightly withtime given the small spatial averaging domain of the minimal channel.The Reynolds shear stress (figure 11 b ) is almost zero above z + > γ . This is different to a half-height channel at a lower friction Reynoldsnumber, in which the impermeability constraint confines the structures at the slip wall.An advantage of higher Re τ with damping compared to lower Re τ without damping isthat the higher friction Reynolds number of the flow with damping will reduce the effectof low friction Reynolds number flows. Chan et al. (2015) showed that the roughness he minimal-span channel for rough-wall turbulent flows z + u ′ + i , r m s z + − h u ′ w ′ i + . . . . ( a ) ( b ) u ′ + x,rms u ′ + z,rms Figure 11: Effect of the outer-layer damping on ( a ) the root mean squared velocityfluctuations in the streamwise and wall-normal directions, and ( b ) the Reynolds stress.Line styles: black, minimal-span, half-height channel (no outer-layer damping); grey,minimal-span, half-height channel with outer-layer damping starting at z + d = 200 shownby ××× ; solid, smooth wall; dashed, rough wall. Vertical dashed line shows roughness crest,vertical dotted line shows critical height z c = 0 . L y .function ∆ U + was overestimated for Re τ ≈
180 compared to Re τ ≈
360 and Re τ ≈ § ∆ p = − /Re τ . It would thereforebe preferable to conduct higher friction Reynolds number simulations with damping,rather than a lower friction Reynolds number simulations. While not investigated here,it also seems reasonable that the grid spacings above z d could be relaxed without alteringthe near-wall flow.This damping is similar to the masks, or sponge regions, employed in Jim´enez & Pinelli(1999). The authors explicitly filtered out disturbances in the outer layer, resulting inlaminar-like flow in the outer layer but still maintaining a healthy near-wall cycle. Thepresent damping is similar, although some fluctuations are still present in the outer layer(figure 11 a ) and the velocity is fixed rather than forming a laminar velocity profile. Abenefit of this damping is that it reduces the computational time step, which is limitedby the CFL number,
CFL = max i ( | u i | ∆ t/ ∆ x i ). Here | u i | is the maximum instan-taneous velocity in the i th direction and ∆ x i is the corresponding grid spacing at thatlocation. Note that the viscous time step restriction does not become significant in anyof the simulations performed here. For minimal-span channels at Re τ = 590 with noouter-layer damping, the maximum CFL number occurs in the outer layer due to thelarge streamwise velocity. We see the time step improves by approximately 20–24% withthe use of the outer-layer damping (table 3, L + y ≈ Re τ = 2000,the time step for the smooth wall with outer-layer damping is unaltered compared tothe case without damping. This is due to the use of the cosine mapping to define thewall-normal grid, which produces an excessively fine grid near the wall especially for highReynolds number cases. As a result of such a small ∆ z , the wall-normal CFL numberis now maximum and this occurs in the near-wall region, which is independent of theouter-layer damping. This clustering of points near the wall is a known issue of the co-sine mapping for high Re τ simulations (Lenaers et al. M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral ID Re τ L x h L + x L + y L z h N x N y N z ∆ x + ∆ y + ∆ z + w ∆ z + h ∆ t + S ∆ t + R Temporal sweep (at initial conditions)MS6 590 2 π Table 4: Initial conditions for temporal sweep in a minimal-span channel. Refer to table1 for definitions.a systematic, controlled study. If a more appropriate mapping was used for these highfriction Reynolds numbers (e.g. that in Lee & Moser 2015), we would expect a similar orbetter improvement in time step as at Re τ = 590.The most appropriate position for z d is difficult to establish from the above data. Aminimum of z + d ≈
200 seems necessary to ensure that the damping does not interactwith the near-wall flow. It also seems reasonable that z d should scale in some way on z c , for larger domain widths. The three sets of simulations suggest that z d should beapproximately two times z c , so a tentative rule of thumb would be z + d & max(200 , z + c ).For an appropriate wall-normal mapping, this allows a 20–24% improvement in the com-putational time step. 3.4. Temporal sweep (table 4)
In this section, the temporal sweep ( § K i = 0) and a streamwise length of L x /h = 2 π (table 4).The initial friction Reynolds number is Re τ = 590 which reduces to Re τ = 180 via anadverse pressure gradient with different rates investigated. The roughness height is fixedat k = h/
40, so the roughness Reynolds number varies as 5 . k + .
15. Three differentsweeps were conducted in which the pressure gradient parameter at the end of the sweepis ∆ p,Re τ =180 = (0 . , . , . a , b ) shows the mean velocity for the slowest sweep, when thepressure gradient parameter (2.6) at Re τ ≈
180 was ∆ p,Re τ =180 = 0 .
03. Figure 12( c , d )shows the mean velocity for the fastest sweep, with ∆ p,Re τ =180 = 0 .
15. To obtain statisticsfor a desired Re τ , time-averaging is performed over a window of ±
10% of that value. Atthe initial friction Reynolds number of the simulation, Re τ = 590, both the sweeps (figure12 a , c ) show good agreement with the full-span channel in the near-wall region. As thespanwise width is L y /h = 0 .
6, then the critical wall-normal location is approximately z + c ≈ . L + y = 142 which agrees well with the figure. As ∆ p is not significant yet (figure2 b ) then there is little difference between the two sweeps.Differences start to emerge at the end of the sweep when the friction Reynolds numberis close to Re τ = 180 and the ∆ p is maximum. The slowest sweep (figure 12 b ) hasreasonable agreement very close to the wall below z + .
15, but the pressure gradienthas reduced the streamwise velocity compared to the full-span steady flow. However, thedifference between the smooth and rough-wall sweep flows (inset of figure 12 b ) agreesquite well with the full-span steady data. Taking the difference between two flows withthe same weak pressure gradient ostensibly still produces a correct estimate of ∆ U + .However, when the pressure gradient is increased to to ∆ p,Re τ =180 = 0 .
15 (figure 12 d )then larger differences are seen. The near-wall region has been noticeably changed, andthe difference in smooth and rough wall velocity (inset) is seen to overestimate the full- he minimal-span channel for rough-wall turbulent flows z + U + z + U + ( a ) ( b ) ∆ p,Re τ =180 =0 . z + U + S − U + R ∆ p,Re τ =180 =0 . z + U + S − U + R z + U + z + U + ( c ) ( d ) ∆ p,Re τ =180 =0 . z + U + S − U + R ∆ p,Re τ =180 =0 . z + U + S − U + R Figure 12: Mean streamwise velocity profile for ( a , c ) Re τ ≈
590 and ( b , d ) Re τ ≈ a , b ) slowest sweep( ∆ p,Re τ =180 = 0 .
03) and ( c , d ) fastest sweep ( ∆ p,Re τ =180 = 0 . § k = h/ U + is a direct result of the stronger(adverse) pressure gradient.The roughness function ∆ U + is now computed from the difference in smooth andrough-wall flows at matched U τ values. As seen in the insets of figure 12, the smooth-and rough-wall velocity difference is not constant with z + as in the steady data. Thisindicates that the flow has not been averaged for long enough at the desired Re τ value.As such, multiple sweeps from different initial conditions would need to be conducted toobtain a converged outer region. However, for the purposes of this investigatory study,the velocity difference is averaged from the crest of the roughness to the channel centre toobtain an estimate of ∆ U + . This is done for a range of friction Reynolds numbers, withthe resulting data plotted in figure 13 against the equivalent sandgrain roughness, k + s .Ideally, the temporal sweep data (lines) should agree with the steady flow data (blacksymbols). Indeed, the temporal sweep does show promise in this regard, especially forthe weaker pressure gradient cases. If we interpolate the steady roughness function datafrom Chung et al. (2015), then an overall relative error term can be approximated as ε = R | − ∆ U + steady / ∆ U + sweep | d k + s . The two slowest sweeps, with ∆ p = 0 .
03 and 0.07have similar relative errors of ε ≈
12% and 13%, respectively, while the fastest sweep isnoticeably larger with an error of ε ≈ M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral k + s ∆ U + Figure 13: (Colour online) Hama roughness function for the temporal sweep with mod-elled roughness at fixed k = h/
40, where k + s ≈ . k + . Symbols: black symbols, steadyflow data for the same modelled roughness (Chung et al. ∆ p,Re τ =180 = 0 .
03 (slowest sweep); , sweep data for ∆ p,Re τ =180 = 0 .
07; , sweepdata for ∆ p,Re τ =180 = 0 .
15 (fastest sweep).It appears that the sweep approach could work for determining the hydraulic behaviourof roughness. A pressure gradient parameter of ∆ p . ∆ p .
01. A value of ∆ p & .
15 results in pressure gradient effectsthat lead to inaccuracies in ∆ U + . Future studies could prescribe a constant value of ∆ p rather than linearly varying the bulk velocity, as this would ensure pressure gradienteffects remain negligible, while efficiently traversing the range of roughness Reynoldsnumbers.
4. Cost and Convergence
Quantifying the cost of the simulations is extremely important when it comes to eval-uating the benefit of minimal-span channels. There are two important costs to considerin these simulations; the cost of memory (or the number of processors required) and thecost of simulation time. The memory cost can be readily described by the size of thegrid: the number of spatial degrees of freedom for full-span simulations scale as Re / τ (Pope 2000, § Re τ & requires national-level high-performance com-puting facilities. The minimal channel, meanwhile, scales as k +9 / s (Chung et al. he minimal-span channel for rough-wall turbulent flows h and velocity U τ . This implies a simulation run time of T sim U τ /h ≈
10 (see, forexample, Hoyas & Jim´enez 2006; Lozano-Dur´an & Jim´enez 2014; Lee & Moser 2015) andthis long runtime is why DNS of full-span channels are so expensive. However, the largestcaptured eddies in the minimal-span channel will be of characteristic size z c , implying weneed T sim U τ /z c ≈
10 eddy turnover times of this z c -sized eddy. In this study, these aregenerally of the size h/z c ≈
10, implying that T sim U τ /z c ≈
10 would be only 1 turnovertime of the h -sized eddy (if it existed), representing an order of magnitude saving in time.This argument, however, does not make reference to the wall-parallel size of the domain.Lozano-Dur´an & Jim´enez (2014) suggested that channels with smaller domain sizes needto be run for proportionately longer, and showed that the statistical uncertainties in u ′ ( z = h ) of full-span channels decrease in a manner inversely proportional to the squareroot of the wall area. However, this result was for two channels that, from the point of viewof the current study, were full-span channels. The largest characteristic eddies in bothchannels were of size h , which explains the simple relationship between channel domainsize and simulation runtime. In contrast, the largest characteristic eddies in minimalchannels are of size z c which varies depending on the width of the channel. This meansthat such a simple relationship between channel domain size and run time is unlikely tohold. In the context of these z c -sized eddies, a full-span simulation could potentially havehundreds of eddies distributed throughout the domain, while a minimal-span simulationwill only have a few. For this minimal-span channel to obtain the same level of convergedstatistics as a full-span channel would presumably require the same number of these z c -sized eddies to pass through the domain. It therefore becomes useful to count thenumber of z c -sized eddies present in the domain. These eddies have a spanwise widthof λ y,z c = L y ≈ z c / . . z c and, as discussed in the previous section, a streamwiselength of approximately three times the span, λ x,z c ≈ λ y,z c = 7 . z c . Hence, the numberof z c -sized eddies at any instant is N instant = L x λ x,z c L y λ y,z c · L z h = L x . z c L y . z c · L z h . (4.1)The first factor is an approximation of the number of z c -sized eddies present on one wallof the channel. The other factor, L z /h , is present as a half-height channel ( L z /h = 1)will have half as many z c -sized eddies as a standard-height channel ( L z /h = 2).Over the course of the simulation, these z c -sized eddies will grow and decay in a processtermed ‘bursting’, so it becomes necessary to consider the time scale of these eddies. Thebursting process in the buffer layer was investigated by Jim´enez et al. (2005), who showedthat the time scale depended on the friction Reynolds number. Low Reynolds numberflows have a long bursting period of T + b ≈ Re τ ≈ T + b ≈
400 for Re τ & L + y . ⇒ z + c . a ) shows this buffer-layer bursting period, along with the current data forbuffer-layer minimal channels ( L + y . Re τ = 590. Good agreement is seen,with the bursting period at this friction Reynolds number being approximately T + b ≈ et al. (2005), definedas the weighted logarithmic average of the frequency spectra of the wall-shear stress.In the logarithmic layer, Flores & Jim´enez (2010) determined that the bursting periodscales with distance from the wall according to T b U τ /z c = 6, which was supported byHwang & Bengana (2016). The bursting period of the largest eddies in the minimalchannel, T + b , therefore depends on both z c and Re τ . Narrow minimal channels whichonly capture the buffer layer ( z + c . Re τ -dependence on the burstingperiod (figure 14 a ), while wider minimal channels which capture part of the logarithmic4 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral Re τ T + b z + c T + b Re τ ≈ Re τ ≈ Re τ > ( a ) ( b ) T +buffer T +log Figure 14: ( a ) Buffer-layer bursting period T + b against friction Reynolds number. Solidline denotes data of Jim´enez et al. (2005), symbols are current data for buffer-layer mini-mal channels. ( b ) Bursting period T + b against wall-normal critical height z + c . Symbols: ××× , Re τ = 590 standard-height channel; +++ , Re τ = 590 half-height channel; ⊕⊕⊕ , Re τ = 2000half-height channel. Lighter grey symbols denote shorter streamwise channel lengths.Line styles: , buffer-layer bursting period at different Re τ ; , log-layer burstingperiod T b U τ /z c = 6 (Flores & Jim´enez 2010; Hwang & Bengana 2016).layer ( z + c & z c . This is summarisedin figure 14( b ), which also shows the present data. The behaviour for intermediate z + c values between buffer and logarithmic flows (grey band) is unknown and likely dependson both z + c and Re τ . Here, the inertial motions of the logarithmic layer are scarce andare likely competing with the Re τ -dependent motions of the buffer layer, with neithercompletely dominating. This distinction becomes irrelevant for higher Re τ channels asthe buffer-layer bursting period saturates to T + b ≈ T b , and the instantaneous number of z c -sizededdies, N instant , gives the total number of z c -sized eddies that exist over the simulationruntime T sim , C ⋆ = T sim T b N instant = T sim T b · L z h · L x . z c L y . z c . (4.2) C ⋆ is simply a running count of how many z c -sized eddies have been sampled, and in-creases with simulation runtime. To investigate whether this is the correct measure forstatistical convergence, the statistical uncertainty of streamwise velocity in different mini-mal channels is analysed. There are various ways of quantifying the statistical uncertainty;here we will estimate it using the standard error of the mean. Instantaneous snapshotsof u ( x, y, z, t i ) at time t i T sim are spatially averaged in the wall-parallel plane to get h u ( z, t i ) i . The temporal average of these N snapshots gives the overall mean, U ( z, T sim ) = 1 N N − X i =0 h u ( z, t i ) i (4.3)The standard error of this mean for statistically independent samples can be computedas ǫ = s/ √ N , where s is the unbiased standard deviation of the signal h u ( z, t i ) i . Thequantity ǫ has the same units as u , so can be non-dimensionalised to ǫ + . A 95% confidenceinterval can then be given in which we are 95% certain the true mean falls within U + ± . ǫ + (Benedict & Gould 1996). The roughness function, being the difference in twoestimated velocities, therefore has a 95% probability of falling within the range ∆ U + ± he minimal-span channel for rough-wall turbulent flows T sim U τ /h ǫ + − − C ⋆ ǫ + − − ( a ) ( b ) (cid:0)(cid:0)(cid:0)(cid:0)✠ Inc. L + x Figure 15: Standard error of the velocity h u i at z + ≈
47 as a function of ( a ) large-eddyturnover time T sim U τ /h and ( b ) number of z c eddies C ⋆ (4.2). Line styles: black, full-span channel; grey, minimal-span channel with z + c ≈
47 where darker grey line indicatesincreasing length L + x (table 1). Dashed line in ( a ) shows ( T sim U τ /h ) − / , the expectedtrend in standard error, and dashed line in ( b ) shows 0 . C ⋆ ) − / . p ǫ +2 s + ǫ +2 r , where subscripts s and r refers to the uncertainties in the smooth-walland rough-wall velocities.In turbulence simulations the samples may not necessarily be statistically independent,so the assumption that ǫ = s/ √ N would not hold. Following Trenberth (1984) andOliver et al. (2014), an effective number of statistically independent samples N eff = N/T can be defined, where T = 1 + 2 N − X k =1 (cid:18) − kN (cid:19) ρ ( k ) (4.4)is the decorrelation separation distance, ρ being the autocorrelation function of h u ( z, t ) i normalised on its variance. The standard error can then be estimated as ǫ = s/ p N eff ,where s = 1 N − T N − X i =0 ( h u ( z, t i ) i − U ( z, T sim )) . (4.5)If the samples were truly independent, then the autocorrelation ρ should be near zero,resulting in T = 1 ⇒ N eff = N and the familiar unbiased standard deviation for s is obtained (with divisor N −
1) . The autocorrelation function can be noisy, espe-cially for small N , prompting Oliver et al. (2014) to develop an open-source code to fitan autoregressive time series model to a given signal. However we have found that forminimal channel simulations N is sufficiently large such that the autocorrelation defini-tion appears to be adequate, with T usually being O (1). The coarse-graining approachfor estimating the standard error for correlated samples, detailed in the appendix ofHoyas & Jim´enez (2008), was also attempted although no appreciable differences weredetected, again because the samples are nearly independent.Increasing the simulation runtime T sim reduces the statistical uncertainties as evi-denced in figure 15( a ), which shows the standard error of the velocity at a representativewall location z + = 47. The very short minimal-span channels (light grey lines) only haveone or two z c -sized eddies present, and the bursting nature of only a few eddies greatly in-creases the statistical uncertainties, requiring a very long runtime. The dashed line shows( T sim U τ /h ) − / , indicating that the uncertainties decrease inversely proportional to the6 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral C ⋆ ǫ + − − C ⋆ ǫ + − − ( a ) z + c ≈ Re τ = 590 z + c ≈ Re τ = 590( b ) C ⋆ ǫ + − − z + c K ( z + c ) − ( c ) z + c ≈ Re τ = 2000 ( d ) Figure 16: Standard error of the velocity as a function of the number of z c -sized eddieseddies C ⋆ (4.2) for ( a ) z + c ≈ Re τ = 590 (table 3), ( b ) z + c ≈ Re τ = 590 (table 1),( c ) z + c ≈ Re τ = 2000 (table 3). Line styles are same as figure 15. Dashed line shows ǫ + = K ( C ⋆ ) − / , where ( a ) K = 0 .
37, ( b ) K = 0 .
23, and ( c ) K = 0 .
23. Coefficient K isshown in ( d ) as a function of z + c . Symbols: +++ , Re τ = 590; ××× , Re τ = 2000. Dashed lineshows K = 33 /z + c .square root of time. This is essentially the same result as Hoyas & Jim´enez (2008) andLozano-Dur´an & Jim´enez (2014) who argued the decrease was inversely proportional tothe square-root of wall area, the difference here being that we have considered the snap-shots as a sampling of time rather than multiples of wall area.When the standard error is instead considered as a function of C ⋆ (figure 15 b ) we obtaina better collapse, indicating that counting z c -sized eddies as in the C ⋆ formulation isindeed an appropriate measure to use. The data show a collapse with ǫ + = K ( C ⋆ ) − / ,where here it was determined K = 0 .
75 for z + c ≈
47. This relationship shows thatobtaining a desired error tolerance requires capturing a certain number of z c -sized eddies( C ⋆ ). Increasing the streamwise domain length does not affect this measure, other thanto simulate more eddies at once.Figure 16 shows that there is not a universal number of z c -sized eddies that need tobe simulated to obtain a desired ǫ + , that is, the coefficient K in ǫ + = K ( C ⋆ ) − / isa function of z + c . Figure 16( a ) shows that for z + c ≈ K = 0 .
37, while figure 16( b )with z + c ≈
142 has K = 0 .
23. Figure 16( d ) shows that this coefficient decreases with z + c ,which implies that a wider minimal channel will have a smaller statistical uncertaintythan a narrower one, for the same C ⋆ value. This is possibly due to the way in whicheddies aggregate in wider channels. As the channel widens and z + c increases, the largest he minimal-span channel for rough-wall turbulent flows z c -sized eddy increases, however there remains a hierarchy of smaller eddies below thislargest eddy. The uncertainty ǫ + of the mean velocity would have contributions from allof these eddies and so would depend on z + c . The reason this effect is not already capturedin C ⋆ is that this quantity only considers the largest z c -sized eddy, whereas ǫ + ostensiblyconsiders all eddies up to and including those of z c size. The data in figure 16( d ) showwhat appears to be a − K ( z c ) = 33 /z + c .Recall that the 95% confidence interval of the roughness function is ∆ U + ± . p ǫ +2 s + ǫ +2 r .If both smooth- and rough-wall simulations were to have the same uncertainty, then thisinterval can be rewritten as∆ U + ± . ǫ + ≈ ∆ U + ± . K ( C ⋆ ) − / ≈ ∆ U + ± . C ⋆ ) − / /z + c . (4.6)The advantage of this formulation is that setting a desired error tolerance in ∆ U + en-ables the simulation run time (and hence CPU hours) to be determined a priori . To seethis, the spanwise channel width, L y & max(100 ν/U τ , k/ . , λ r,y ) (Chung et al. L x & max(3 L y , ν/U τ , λ r,x ) ( § z c = 0 . L y . The user must then determine a desired error tolerance, ζ , for ∆ U + ± ζ . This error tolerance will depend on the application; riblets have a rough-ness function − ∆ U + . U + ≈ C ⋆ through ∆ U + ± . C ⋆ ) − / /z + c . From here,(4.2) can be solved for the simulation time T sim using the minimal channel dimensionsand bursting period T b . The number of CPU hours is thenCPU hours = T + sim ∆ t + T N CP U , (4.7)where ∆ t + is the computational time step, so T + / ∆ t + is the number of time stepsperformed during the simulation. T is the average wall clock time taken to perform asingle step and N CP U is the number of processors employed in the computation. The timestep ∆ t + can be estimated if the streamwise CFL number at the channel centre is unity,i.e. U h ∆ t/ ∆ x = CF L = 1, which gives ∆ t + = ∆ tU τ /ν = ∆ x + U τ /U h . If ∆ x + ≈
10 and,for a full-span channel, the centreline velocity is U h /U τ ≈
25 then we have ∆ t + ≈ . T and N CP U will depend on the code being used. The full-span channel usedin this study, run for
T U τ /h = 10 large-eddy turnover times with ∆ t + ≈ . T ≈ . N cpu = 64 requires approximately 15800 CPU hours. For a minimal channelwith h/z c ≈
10 and 95% confidence interval of ∆ U + ± .
1, the above analysis leads toa requirement of only 1420 CPU hours (for N CP U = 64 and T = 2 .
5. Minimal channel applied to pyramids
Computational setup
Following the insights discussed in this paper, a finite-volume code is used to deter-mine the roughness function of square-based pyramids. This is the same code as inChan et al. (2015) and Chung et al. (2015), based on the work of Mahesh et al. (2004)and Ham & Iaccarino (2004). No outer-layer damping will be used ( K = ) and since abody-fitted grid is used, no roughness forcing is required ( F = ). Three pyramids aresimulated with different heights, but all with nominal Re τ = 840 and angles between8 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral
Case Re τ k + t λ + α L + x L + y N x N y N z ∆ x + ∆ y + ∆ z + w ∆ z + h C f,mod ∆ U + Sm 198 834 - - - 1180 197 288 48 320 4.1 4.1 0.14 6.2 0.0053 -Sm 313 850 - - - 1201 316 288 72 320 4.2 4.4 0.14 6.3 0.0051 -Sm 396 851 - - - 1204 401 288 96 320 4.2 4.2 0.14 6.3 0.0050 -40 198 848 40.4 200 22 1200 200 288 48 320 4.2 4.2 0.14 6.3 0.015 8.040 198 f 837 39.2 197 22 1184 197 504 84 560 2.3 2.3 0.24 3.6 0.015 7.940 198 2 879 41.7 207 22 1243 414 288 96 320 4.3 4.3 0.15 6.5 0.015 8.363 313 883 66.3 329 22 1249 329 288 72 320 4.3 4.6 0.24 6.5 0.021 10.280 396 868 82.4 409 22 1227 409 288 96 320 4.3 4.3 0.24 6.4 0.024 11.0SF09 4020 39.9 198 16 - - - - - - - - - - 6.4SF09 3900 62.7 310 16 - - - - - - - - - - 7.0SF09 4290 87.5 433 16 - - - - - - - - - - 8.7DCO16 670 39 147 28 - - - - - - - - - - 8.9HKS11 3520 63.3 442 16 - - - - - - - - - - 8.2
Table 5: Description of the simulations performed with the finite volume code for thesquare-based pyramids. k + t is the trough-to-peak pyramid height, λ + the pyramid wave-length, α the pyramid edge angle, and C f,mod = 2 U τ /U b,full the predicted full-span skinfriction coefficient (see text). Other entries are same as table 1. All simulations conductedusing a half-height channel with no outer-layer damping. SF09 refers to Schultz & Flack(2009), HKS11 to Hong et al. (2011), and DCO16 to Di Cicca & Onorato (2016).the pyramid edge and the horizontal of α = 22 ◦ . Each pyramid will be referred to byits nominal trough-to-peak height, k + t and wavelength λ + , as a unique identifier givenby the notation k + t λ + . Hence the case 40 198 has a height k + t ≈
40 and wavelength λ + ≈ ◦ are given in table 5. This makes the present pyramids steeper thanSF09, however their study included pyramids with α = 35 ◦ and observed no appreciabledifference in ∆ U + , suggesting this discrepancy is insignificant. The smallest pyramid alsohas a similar k + t to Di Cicca & Onorato (2016), herein referred to as DCO16, who had k + t = 39. Their pyramid slope angle was α = 28 ◦ and their pyramids were separatedfrom each other by a small gap of 32 viscous units. However, given that the roughnessheight is matched, these difference shouldn’t result in a substantially different rough-ness function. Additionally, the case 63 313 matches the roughness height of Hong et al. (2011) (herein HKS11). The cases that are simulated here, and those available in theliterature, are summarised in table 5. The expected full-span skin-friction coefficient, C f,mod = 2 U τ /U b,full , is also given in this table. Due to the altered outer-layer velocityprofile of the minimal-span channel, the composite profile of Nagib & Chauhan (2008)is fitted to estimate the full-span velocity profile for z > z c . This therefore allows foran estimate of the full-span bulk velocity U b,full to be obtained. The estimated smooth-wall coefficient value of C f,mod ≈ .
005 is in good agreement with the empirical fit ofDean (1978), in which C f = 0 . Re b ) − / = 0 . Re b = 2 hU b /ν is the bulkReynolds number. This shows that the minimal-span channel framework can still be usedto estimate the full-span skin-friction coefficient, which is commonly used by engineers.Now that the dimensions of the roughness are known, the channel domain size can he minimal-span channel for rough-wall turbulent flows x + y + z + . . − . x + y + z + . . − . − . . ( a ) ( b ) Figure 17: (Colour online) Pyramid geometry and body-fitted mesh for pyramid with k + t = 40, λ + = 198 at Re τ = 840 for ( a ) normal mesh 40 198 and ( b ) finer mesh40 198 f. Only every ( a ) fifth or ( b ) tenth wall-normal node is shown, and only everysecond spanwise and streamwise node in ( b ). Contour shows time-averaged total drag(viscous + pressure) force per unit area.be determined. Following the guidelines of Chung et al. (2015), the spanwise domainwidth must satisfy L y & max(100 ν/U τ , k t / . , λ r,y ). For the smallest pyramids (case40 198), we have L + y & max(100 , , L y = λ . An additional case with the same channel dimensions as 40 198 is simulated (case40 198 f), but with a finer computational grid to ensure that all the turbulent scales areresolved. This spanwise width of L y = λ means that the critical wall-normal location z c = 0 . L y for these pyramids is set at z + c = 80 ≈ k + t . The roughness sublayer maybe larger than 2 k t , so a channel with a larger width of L y = 2 λ is also simulated (case40 198 2), so that the critical location is now z c = 4 k t . For the larger pyramids (cases63 313 and 80 396), the pyramid wavelength is again the limiting constraint, so L y = λ . The streamwise domain length was investigated in § L x & max(3 L y , ν/U τ , λ r,x ). For the smaller pyramids (40 198), we have L + x & (594 , , L x = 6 λ isselected which in viscous units is L + x = 1188. Case 63 313 is also limited by the secondconstraint. For the largest pyramids ( L + y = 396), the first constraint ( L x & L y ) requiresa streamwise length of L x = 3 L y ⇒ L + x = 1188.The length of time to achieve a converged result can now be estimated since thespanwise width L y and hence critical wall-normal location z c = 0 . L y have been set. Forthe smallest pyramid case (40 198), the spanwise width is L + y = 198 so z + c = 79. If we wishto be 95% confident that the roughness function is in the range ∆ U + ± .
1, then from (4.6)the number of z c -sized eddies that must be observed are C ⋆ = (91 . / (0 . z + c )) ≈ T sim through (4.2), which requiresestimating the bursting period T b . Since the channel is relatively wide and z + c = 79, it’slikely to have the beginnings of a logarithmic layer (figure 14 b ), which implies T b U τ /z c =6 ⇒ T + b = 6 z + c ≈ T + sim ≈ . × . As the nominal streamwisegrid spacing is known to be ∆ x + ≈ .
2, and if we assume a centreline velocity of U + h ≈ t + = ∆ x + /U + h ≈ .
14 for a CFL number of unity at the centreline.Finally, using (4.7) with an assumed wall clock time per step of T = 6 . N CP U = 128, we expect to use around 52,000 CPU hours. A similar computation forcases 63 313 ( z + c ≈ C ⋆ ≈
53) and 80 395 ( z + c ≈ C ⋆ ≈
33) predicts 53,000 and0
M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral L x = 2 π , L y = π , L z = 2 h and T sim U τ /h = 10) with pyramids would likely require at least 2 × CPU hours using the current code, more than a full order of magnitude greater than theminimal channel.The pyramids are aligned so that the trough between neighbouring pyramids is 45 ◦ tothe mean flow direction. However, because hexahedral cells are used which are alignedwith the mean flow direction, the trough line falls across opposite vertices of the cell.This would mean the cell needs to be split into two pentahedrals to fit the roughnessgeometry. To mitigate this issue, the trough of the pyramids are rounded off and faceted.The resulting trough has a radius of curvature of approximately half the pyramid crest-to-peak height, R = 0 . k t . The crest of the pyramid retains its sharp edge, as evidencedin figure 17. The total volume of fluid in the channel with pyramids matches that of asmooth-wall channel, which would ensure a collapse of flow statistics in the outer layer,if it existed (Chan et al. Pyramid results
Figure 18( a – c ) shows the mean velocity profile for the pyramid roughness data. The insetsshow the difference in velocity between smooth-wall and pyramid-roughness channels, andit can be seen to be relatively constant from the pyramid crest to the channel centre.This indicates that the velocity profile for the smooth-wall and rough-wall are matchedfrom the crest (apart from the offset ∆ U + ), or that the roughness sublayer for pyramidsis small. The roughness function is obtained by averaging the difference in smooth- andrough-wall velocities over k + t + 5 < z < z + c .For the smaller pyramids ( k + t = 40) in a narrow channel, we obtained a roughnessfunction of ∆ U + = 8 . U + = 7 . U + = 8 .
3. While slightly larger than theprevious two values, this difference is predominantly due to the increase in Re τ for thischannel. The roughness Reynolds number, k + t ≈ .
8, is almost 5% over the target valueof k + t = 40, which means the roughness function would also be larger. The otherwisegood agreement between this case and the narrower case 40 198 supports the conclusionthat the roughness sublayer is less than 2 k t in height above the pyramid crest. If theroughness sublayer was larger than 2 k t , then the narrower minimal channel with z c = 2 k t could not fully capture the roughness effects and hence would have a different roughnessfunction to the wider minimal channel. Moreover, this agreement with case 40 198, when L y = λ y , suggests that the interaction of flow structures generated between neighbouringroughness elements is small, as this case only has a single repeating element in thespanwise direction. This agrees with our previous study (MacDonald et al. d -type surfacesthat these spanwise effects would become significant. However, this would require furtherstudy and is outside the scope of this paper. The above three values of the roughnessfunction can be compared with values of 6 . ±
10% from SF09 and 8 . d ).For the case with k + t ≈
63, the roughness function was found to be ∆ U + = 10 .
2, which he minimal-span channel for rough-wall turbulent flows z + U + z + U + ( a ) ( b ) k + t ≈ k + t ≈ z + U + S − U + R z + U + S − U + R z + U + k + t ∆ U + ( c ) k + t ≈
80 ( d ) z + U + S − U + R Figure 18: ( a – c ) Mean velocity profile for pyramid roughness. Line styles: full-spanchannel, Re τ = 950 (Chung et al. a ) k + t ≈
40, ( b ) k + t ≈
63, ( c ) k + t ≈ +++ indicate wall-normal critical location z c = 0 . L y . Inset shows difference in smooth-walland pyramid-roughness velocity. ( d ) Roughness function against trough-to-peak pyramidheight, k + t . Symbols: N , SF09 data in table 5; △△△ , other SF09 data for α = 16 ◦ pyramids; (cid:4) , DCO16 data; (cid:7) , HKS11 data; , Nikuradse (1933) sandgrain data; N , current data(40 198, 63 313, 80 396); ◮ , 40 198 2; ◭ , 40 198 f.is larger than 7 . ±
10% from SF09 and 8.2 from HKS11. For the largest pyramids, where k + t ≈
80, the roughness function was found to be ∆ U + = 11 .
0, which is again largerthan the value of 8 . ±
10% from SF09. As a result, the equivalent sandgrain roughnessappears to be k s ≈ . k t , compared to k s ≈ . k t from SF09. Note that DCO11 wouldpredict an even larger equivalent sandgrain roughness than k s ≈ . k t , given their larger∆ U + at k + t ≈ S = ( { u ′ w ′ } , − { u ′ w ′ } , ) / h u ′ w ′ i , was developed inRaupach (1981) in which the Reynolds shear stress is conditionally averaged on sweeps( u ′ > w ′ <
0, quadrant 4) and ejections ( u ′ < w ′ >
0, quadrant 2), where {·} denotes the conditional average. The hyperbolic hole region H = 0. From this statistic,it can be seen in figure 19( a ) that sweeps are the dominant means of momentum transferwithin the roughness canopy, in agreement with Raupach (1981). These sweeps are muchmore dominant than in a smooth wall, which has a slight preference for sweeps very closeto the the wall ( z + . M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral z + ∆ S − . − .
25 0 0 .
25 0 . − z / k t ∆ S − . − .
25 0 0 .
25 0 . . . ( a ) ( b ) ✻ Inc. k + t Figure 19: (Colour online) Difference ∆ S between sweep and ejection stress contributions(Raupach 1981) as a function of ( a ) z + and ( b ) z/k . Line styles: , smooth wall(sm 396); , k + t ≈
40 (40 198); , k + t ≈
63; , k + t ≈
80. Symbols: +++ , z + c ; ××× ,pyramid crest. C ⋆ ǫ + − − C ⋆ ǫ + − − ( a ) ( b ) Figure 20: Standard error of the velocity as a function of the number of z c -sized ed-dies eddies C ⋆ (4.2) for ( a ) z + c ≈
79, and ( b ) z + c ≈ k + t ≈
40 pyramid; , k + t ≈
80 pyramid. Horizontal dashed line shows ǫ + =0 . / .
77 = 0 . C ⋆ =(91 . / (0 . z + c )) , the expected C ⋆ value required to obtain this tolerance level (4.6).both rough-wall and smooth-wall flows collapse with ejections being slightly larger thansweeps. The outer-layer is associated with ejections being dominant, with ∆ S tendingto − z = h . However, in the minimal-span channel this statistic stays approximatelyconstant above z > z c , with ∆ S being close to zero. If instead the wall-normal coordi-nate is normalised on k t (figure 19 b ) then we see that all three rough-wall flows nearlycollapse, although the case with the smallest height ( k + t ≈
40, red dashed line) there isa slight difference, with sweeps remaining dominant at the roughness crest.Figure 20 shows the statistical uncertainty of U + ( z = z c ) as a function of C ⋆ , thenumber of z c -sized eddies. This is the same format as figure 16( a – c ). Importantly, thepredicted number of z c -sized eddies that need to be observed to obtain an error toleranceof ∆ U + ± . C ⋆ value. While this means an he minimal-span channel for rough-wall turbulent flows § Re τ and roughness geometry then what thepredicted C ⋆ empirical relation was developed with, showing a certain robustness to themethod.The roughness function values of 8 . ∆ U + .
11 measured here place these pyramidcases in the fully rough regime, while the results and recommendations from the previoussections ( §
6. Conclusions
A series of numerical experiments were conducted to investigate the fundamental dy-namics of the minimal channel, with an emphasis on more efficiently simulating rough-wall flows in order to obtain the roughness function ∆ U + . These experiments were per-formed using a finite difference code.The streamwise length was investigated in § L x = 3 L y , where the spanwise width is obtained fromthe guidelines in Chung et al. (2015), L y & max(100 ν/U τ , k/ . , λ r,y ). However, for verynarrow channels this implies a very short streamwise domain length, which may not beable to sustain the turbulent motions. The small spatial domain also makes obtainingconverged statistics a very time-consuming process. As such, minimal channels wouldbenefit from a minimum streamwise domain length of around 1000 viscous units and soa guideline for setting L x is L x & max(3 L y , ν/U τ , λ r,x ), where λ r,x would be somecharacteristic streamwise length scale of the roughness.Two alterations to the outer-layer flow were investigated to see whether they affectedthe healthy near-wall flow. A half-height (open) channel was used in § dU/dz . When dU/dz was instead artificially set to zero in theunphysical region above z c , the pressure fluctuations from this altered velocity field werethen in much better agreement with the full-span channel. A benefit of using a half-heightchannel is that only half the number of gridpoints are required, at the cost of runningthe simulation for twice as long to reach the same level of statistical convergence. It isup to the user to balance this trade off between memory and time.A forcing model was applied to the outer layer in § z d . This is similar tothe filtering employed by Jim´enez & Pinelli (1999), although the present damping modelforces U ( z > z d ) = U ( z = z d ). This allowed for an improvement in the computationaltime step of around 20–24%, although this did not extend to higher friction Reynoldsnumbers in the present simulations due to the grid using a cosine mapping in the wall-normal direction. The present data suggest that the location where the forcing startsshould be z d & max(200 ν/U τ , z c ). The first constraint is present for very narrow chan-4 M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral nels where z + c is close to the wall, as the forcing will interfere with the near-wall flowotherwise.A temporal sweep was conducted in § ∆ p,Re τ =180 = 0 .
03 and 0.07 could reasonably predict some of the ∆ U + vs k + s curve. The largest value of ∆ p = 0 .
15 resulted in substantial flow changes comparedto the steady flow and so are not feasible. The early work of Perry & Joubert (1963)showed that ∆ U + was independent of an applied pressure gradient, and the presentwork goes on to provide a bound on ∆ p for when this holds. This also suggests thatpressure gradients, normally applied by a spatial acceleration of the flow in experimentalfacilities, can equally be simulated by a temporal acceleration. This is similar to theapproach of Kozul et al. (2016) who represented a spatially developing boundary layeras temporally developing instead.An eddy-counting argument was developed in § z c -sized eddies, C ⋆ , arecounted over the duration of the simulation, which involves estimating their characteristiclength and time scales. For a desired level of uncertainty, ζ , in the roughness function,∆ U + ± ζ , the expected number of z c -sized eddies that need to exist over the course of thesimulation was estimated as C ⋆ = (91 . / ( ζ · z + c )) as in (4.6). This means that for a known z + c , the user can determine how many z c -sized eddies need to be observed through thisempirical relation. The simulation run time, and hence number of CPU hours required,can then be estimated a priori . This would enable researchers to determine beforehandhow best to allocate a limited number of CPU hours when studying rough-wall flows.A case study of square-based pyramids ( §
5) was used to illustrate the minimal-spanchannel framework and the insights gained in this paper. The viscous dimensions of thepyramids were set to match those of Schultz & Flack (2009). The roughness functionsfrom these simulations were overestimated compared to those of Schultz & Flack (2009),although underestimated compared to that of Di Cicca & Onorato (2016) who had asimilar roughness geometry. The predicted C ⋆ value required to obtain a desired level ofstatistical uncertainty was shown to be in reasonable agreement, if not somewhat conser-vative, with the data from the pyramid simulations. As such, the estimated CPU hoursrequired for the simulation can be predicted reasonably accurately before performingthe simulation. For this pyramid roughness, the minimal-span channel used nearly 20times less CPU hours than the estimated cost of a full-span channel, as well as usingsubstantially fewer CPUs. Acknowledgements
This work was partly funded through the Multiflow program by the European Re-search Council. Computational time was granted under the Victoria Life Sciences Com-putational Initiative, which is supported by the Victorian Government, Australia.
REFERENCESdel ´Alamo, J. C., Jimenez, J., Zandonade, P. & Moser, R. D
J. Fluid Mech. , 329–358.
Benedict, L. H. & Gould, R. D.
Exp. Fluids (2), 129–136. he minimal-span channel for rough-wall turbulent flows Bernardini, M., Pirozzoli, S., Quadrio, M. & Orlandi, P.
J. Comput. Phys. (1), 1–6.
Borrell, G.
Busse, A. & Sandham, Neil D.
J. Fluid Mech. , 169–202.
Chan, L., MacDonald, M., Chung, D., Hutchins, N. & Ooi, A.
J. Fluid Mech. , 743–777.
Chin, C., Ooi, A. S. H., Marusic, I. & Blackburn, H. M.
Phys. Fluids (11), 115107. Choi, H. & Moin, P.
J. Comput. Phys. (1), 1–4.
Chung, D., Chan, L., MacDonald, M., Hutchins, N. & Ooi, A.
J. Fluid Mech. , 418–431.
Chung, D., Monty, J. P. & Ooi, A.
J. Fluid Mech. , R3.
Dean, R. B.
J. Fluids Engng. (2), 215–223.
Di Cicca, G. M. & Onorato, M.
J. Turbul. (1), 51–74. Flack, K. A., Schultz, M. P. & Shapiro, T. A.
Phys. Fluids (3), 035102. Flores, O. & Jim´enez, J.
J. Fluid Mech. , 357–376.
Flores, O. & Jim´enez, J.
Phys.Fluids (7), 071704. Ham, F. & Iaccarino, G.
Annual Research Briefs 2004 , pp. 3–14. Center for TurbulenceResearch, Stanford University/NASA Ames.
Hama, F. R.
Trans. Soc.Naval Arch. Mar. Engrs , 333–358. Hamilton, J. M., Kim, J. & Waleffe, F.
J. Fluid Mech. , 317–348.
Handler, R. A., Saylor, J. R., Leighton, R. I. & Rovelstad, A. L.
Phys. Fluids (9), 2607–2625. Hong, J., Katz, J. & Schultz, M. P.
J. Fluid Mech. , 1–37.
Hoyas, S. & Jim´enez, J. Re τ =2003. Phys. Fluids (1), 011702. Hoyas, S. & Jim´enez, J.
Phys. Fluids (10), 101511. Hwang, Y.
J. FluidMech. , 264–288.
Hwang, Y.
J. Fluid Mech. , 254–289.
Hwang, Y. & Bengana, Y.
J. Fluid Mech. , 708–738.
Jeong, J., Hussain, F., Schoppa, W. & Kim, J.
J. Fluid Mech. , 185–214.
Jim´enez, J.
Annu. Rev. Fluid Mech. , 173–196. Jim´enez, J.
Phys. Fluids (1), 065102. Jim´enez, J., Kawahara, G., Simens, M. P., Nagata, M. & Shiba, M.
Phys. Fluids (1),015105. M. MacDonald, D. Chung, N. Hutchins, L. Chan, A. Ooi and R. Garc´ıa-Mayoral
Jim´enez, J. & Moin, P.
J. Fluid Mech. , 213–240.
Jim´enez, J. & Pinelli, A.
J. Fluid Mech. (1), 335–359.
Jones, M. B., Marusic, I. & Perry, A. E.
J. Fluid Mech. , 1–27.
Kim, J.
J.Fluid Mech. , 421–451.
Kim, J. & Moin, P.
J. Comput. Phys. (2), 308–323. Kim, K. C. & Adrian, R. J.
Phys. Fluids (2), 417–422. Kline, S. J., Reynolds, W. C., Schraub, F. A. & Runstadler, P. W.
J. Fluid Mech. (04), 741–773. Kozul, M., Chung, D. & Monty, J. P.
J. Fluid Mech. , 437–472.
Lee, M. & Moser, R. D. Re τ ≈ J. Fluid Mech. , 395–415.
Lenaers, P., Schlatter, P., Brethouwer, G. & Johansson, A. V.
J. Comput.Phys. , 108–126.
Leonardi, S. & Castro, I. P.
J. Fluid Mech. , 519–539.
Lozano-Dur´an, A. & Jim´enez, J. Re τ = 4200. Phys. Fluids (1), 011702. Lundbladh, A., Berlin, S., Skote, M., Hildings, C., Choi, J., Kim, J. & Henningson,D. S.
TRITA-MEK
MacDonald, M., Chan, L., Chung, D., Hutchins, N. & Ooi, A.
J. Fluid Mech. , 130–161.
Mahesh, K., Constantinescu, G. & Moin, P.
J. Comput. Phys. , 215–240.
Marusic, I., Monty, J. P., Hultmark, M. & Smits, A. J.
J. Fluid Mech. , R3.
Mizuno, Y. & Jim´enez, J.
J. Fluid Mech. , 429–455.
Munters, W., Meneveau, C. & Meyers, J.
Phys. Fluids (2), 025112. Nagib, H. M. & Chauhan, K. A.
Phys. Fluids , 101518. Nikuradse, J.
Oliver, T. A., Malaya, N., Ulerich, R. & Moser, R. D.
Phys. Fluids (3), 035101. Patel, V. C.
J. Fluid Mech. (01), 185–208. Perry, A. E. & Joubert, P. N.
J. Fluid Mech. (02), 193–211. Pope, S. B.
Turbulent Flows . Cambridge University Press.
Raupach, M. R.
J. Fluid Mech. , 363–382.
Raupach, M. R., Antonia, R. A. & Rajagopalan, S.
Appl. Mech. Rev. (1), 1–25. Schultz, M. P. & Flack, K. A.
Phys. Fluids , 015104. Seddighi, M., He, S., Vardy, A. E. & Orlandi, P.
Flow Turbul. Combust. (1-2), 473–502. Sekimoto, A., Dong, S. & Jim´enez, J. he minimal-span channel for rough-wall turbulent flows tionary and homogeneous shear turbulence and its relation to other shear flows. Phys.Fluids , 035101. Spalart, P. R. & McLean, J. D.
Phil. Trans. R. Soc. A (1940), 1556–1569.
Spalart, P. R., Moser, R. D. & Rogers, M. M.
J. Comput. Phys. (2), 297–324. Toh, S. & Itano, T.
J. Fluid Mech. , 249–262.
Townsend, A. A.
The Structure of Turbulent Shear Flow , 2nd edn. Cambridge UniversityPress.
Trenberth, K. E.
Mon. Weather Rev.112