Statistical Structure and Deviations from Equilibrium in Wavy Channel Turbulence

The structure of turbulent flow over non-flat surfaces is a topic of major interest in practical applications in both engineering and geophysical settings. A lot of work has been done in the fully rough regime at high Reynolds numbers where the effect on the outer layer turbulence structure and the resulting friction drag is well documented. It turns out that surface topology plays a significant role on the flow drag especially in the transitional roughness regime and therefore, hard to characterize. Survey of literature shows that roughness function depends on the interaction of roughness height, flow Reynolds number and topology shape. In addition, if the surface topology contains large enough scales then it can impact the outer layer dynamics and in turn modulate the total frictional force. Therefore, it is important to understand the mechanisms underlying drag increase from systematically varied surface undulations in order to better interpret quantifications based on mean statistics such as roughness function. In this study, we explore the mechanisms that modulate the turbulence structure over a two-dimensional (2D) sinusoidal wavy surface with a fixed amplitude, but varying slope. To accomplish this, we model the turbulent flow between two infinitely wide 2D wavy plates at a friction Reynolds number, Re τ = 180. We pursue two different but related flavors of analysis. The first one adopts a roughness characterization flavor o f such w avy s urfaces. T he second o ne focuses o n u nderstanding the n on-equilibrium near surface turbulence structure and their impact on roughness characterization. Analysis of the different statistical quantifications show strong dependence on wave slope for the roughness function indicating drag increase due to enhanced turbulent stresses resulting from increased production of vertical velocity variance from the surface undulations.


Introduction
Surface undulations can have significant impact on turbulent boundary layers both in the atmosphere as well as in engineering applications.In particular, engineering applications such as internal flows in pipes and turbomachinery, external flows over fouled ship hulls [1], wind turbine blades and other aerodynamic surfaces are common examples.In the atmospheric side, while most roughness such as grass and shrubs are very small, there exist medium to large scale roughness in the form of tree canopies, man made structures and hills.The ubiquitous nature of such flows has made understanding their dynamics a necessity.A significant amount of research has been devoted to understanding turbulent flows over pipe roughness, for example, the work of Darcy [2] nearly two hundred years ago, in the early half of twentieth century by Nikuradse [3], Colebrook [4] and Moody [5] and more recently by various research groups [6][7][8].In the last two decades, fundamental investigation of turbulent flows over uniform roughness embedded in flat surface has been undertaken through a series of experimental studies [9][10][11][12][13][14][15] as reviewed in [13,16].In addition there has been extensive simulation-based research of turbulent boundary layers over systematically designed roughness using direct numerical simulation (DNS) [8,17,18] and large eddy simulation (LES) [19].There has also been interesting recent work on reproducing Nikuradse-type sand grain roughness using DNS at moderately high Reynolds numbers [20,21].
Through this extensive and growing body of literature, the underlying goals and fundamental questions remain consistent, namely, how to estimate flow drag over a given roughness topology at a specified Reynolds number or flow rate.From a geophysical perspective, the goal is to model the outer layer dynamics and understand the turbulent coherent structures within the roughness sublayer that impact man-made applications in the lower atmosphere [22][23][24].From a computational standpoint, the question is one of modeling the effective dynamics within the roughness sublayer to bypass the complexity of resolving the roughness elements.
Significant early attempts to answer some of the above questions were the work of Nikuradse [3] and the subsequent extension by Colebrook [4] to relate flow drag with roughness.Both these efforts classify roughness as hydraulically smooth, transitional or fully rough regimes depending on the relationship between drag and roughness scales.In the fully rough regime, drag is independent of the Reynolds number and depends only on the roughness scale whereas in the transitional regime both of these are important as per [3,4].These ideas are summarized in the popular Moody diagram [5].A more generic quantification of roughness induced effects applicable across different classes of flows is the Hama roughness function [25], ∆ u + which is commonly aligned with the classical view of rough wall turbulent boundary layers.Specifically, the classical view is that roughness influences the turbulence structure only up to a few roughness lengths from the mean surface location while the outer layer flow is unaffected except for a modulation in the velocity and length scales -a rough wall extension of Townsend's Reynolds number similarity hypothesis [26].Therefore, this notion of 'wall similarity' [27] implies that shape of the mean velocity in the overlap and outer layers is unaffected (relative to a smooth wall) by the roughness.This phenomenology is mostly consistent with observations as per Jiménéz [16], but exception do exist.Quantitatively, the roughness function represents the downward displacement in the mean velocity profile plotted in a semi-log scale indicative of the increased drag from the surface inhomogeneities.Combined with Townsend's wall similarity hypothesis, ∆ u + represents the shift in the intercept used to describe the logarithmic region of the mean velocity profile as Log law for smooth wall where, κ is the von Kármán constant, u + is the averaged streamwise velocity over a rough surface and y + is the wall coordinate.Normalization is done using the inner layer variables such as friction velocity, u τ and kinematic viscosity, ν expressed as u + = u u τ and y + = yu τ ν .
conditions for the existence of outer layer similarity.In their work, outer layer similarity was consistently encountered even for moderately high Reynolds numbers over uniform three-dimensional rough surfaces with reasonably large roughness elements.In fact, little to no deviations in outer layer similarity was observed for δ/k 20 and δ/k s 6 where k, k s are the mean roughness height and equivalent sand grain roughness respectively.Importantly, it was reported that there exists no critical roughness height (i.e. a height corresponding to a sharp transition) as the influence of roughness size on outer layer statistics is more gradual.These trends appear to breakdown for two-dimensional roughness (especially periodic) which are known to generate stronger vertical disturbances due to the absence of significant spanwise motions in the roughness sublayer.This two-dimensional surface effect is clearly observed in higher order statistics and less so for the mean velocity profiles.Volino et al. [28,29] clearly illustrate this using experiments with transverse two-dimensional bars and three-dimensional cubes as roughness elements (δ/k s ≈ 2 − 3) and flow friction Reynolds numbers, δ + ≈ 2000.Subsequently, Krogstad and Efros [30] show that flow over two-dimensional roughness elements with higher δ/k and at higher Re τ (larger scale separation) generate outer layer similarity just like three-dimensional roughness.Therefore, the details of the roughness topology along with the roughness scale and the flow Reynolds number modulate turbulence structure.It is only their relative importance that changes across the different regimes.
The frictional drag from the surface is also strongly influenced by the surface topology and not just the roughness scale.The correlation of friction coefficient with Reynolds number in the Moody diagram [5] is parameterized for mean roughness height.By leveraging the existence of outer layer similarity, the shift in the mean profile or roughness function is used as a surrogate for prediction of increase in friction drag over a given rough surface.The roughness correlations of Nikuradse [3] or Colebrook [4] relate respectively.The Nikuradse expression is calibrated for the uniform sand grain roughness while the Colebrook relationship is designed for commonly occurring surfaces.However, there exists many examples where ∆ u + does not depend solely on k + such as a turbulent flow over two-dimensional transverse (or wavy surfaces as reported in this work) bars [31] separated by distances comparable to the height.
Perry et al. [32] classify such cases as "d-type" roughness in contrast to "k-type" roughness where ∆ u + scales with k + .In the case of closely spaced transverse bars, there exist vortical flow cells between each of these elements thus causing the turbulent flow to skim over which makes ∆ u + depend little if any on k + .
Schultz and Flack [12] systematically study turbulent flow over three-dimensional pyramid elements of different heights and inclination angles to understand the role of roughness slope in addition to k on the drag.Their results clearly indicate that smaller slopes produce significant deviation from the uniform sand roughness behavior with ∆ u + changing slower with k + in the transitionally rough regime.Further, these deviations get stronger with increase in inner scaled roughness height.Nakato et al. [33] report similar observations for sinusoidal wavy surfaces with slopes greater than ∼ 6 • mimicking the uniform sand roughness behavior.Physical insight for these observations is available from the numerical work of Napoli et al. [17] who superposed sinusoidal waves to generate a corrugated two-dimensional rough surface.
From their work, the anomalous relationship between ∆ u + and k + at small slopes or 'waviness' regime is attributed to the dominance of viscous drag over form drag. On the contrary, in the 'roughness' regime involving higher slopes (as in [3,14,20]), the form drag dominates.To characterize these deviations, they design a slope dependent roughness parameter termed as effective slope, ES, that represents the average surface slope magnitude over a given sampling region.In their findings, ∆ u + varies linearly with ES for ES < 0.15 (beyond which ∆ u + is constant for a given k + ) whereas Schultz and Flack [12] report that the transition happens at ES < 0.35.Therefore, ES is an additional 'waviness' parameter along with k + that modulates ∆ u + , i.e. ∆ u + = f (k + , ES).Of course, one can build a rich enough parameter space in addition to ES and a to learn f using advanced data science methods.
In this study we explore the structure of near-wall turbulence and deviations from equilibrium flat channel turbulence in the waviness regime using direct numerical simulation of wavy channel flow at a friction Reynolds number, Re τ = δ + ≈ 180.The simulation infrastructure uses higher-order spectral like compact schemes [34] for both advection and diffusion terms while a third-order multi-step method is used for time integration.The wavy surface is represented using an immersed boundary method [35,36] similar to many other efforts [18,21].The focus of our current analysis is to better understand the mechanisms underlying the drag increase at small slope angles dominated by viscous drag.For the sinusoidal two-dimensional surfaces considered in this study, the effective slope, ES is directly related to the non-dimensional ratio of the amplitude (a) and wavelength (λ) of the sinusoid, i. on the three-dimensional surface roughness studies of Flack et al. [13,14], but may not be adequate for the two-dimensional wavy surfaces used in this study.Therefore, analysis in this work will focus on assessing the extent of outer layer similarity and the relationships between roughness function, effective slope and roughness/wave height.In addition, we delve into the nature of roughness induced deviations on higher order turbulence statistics and their production mechanisms in order to generate a process level understanding.Note that the roughness Reynolds numbers used in the work (a + ≈ 13 − 14) fall within the transitional regime.
The rest of the article is organized as follows.In section 2, we describe the numerical methods, simulation design, quantification of statistical convergence and validation efforts.In section 3, we present the results from the analysis of outer layer similarity, roughness induced drag quantification.We further characterize how the turbulence structure, namely, the flow stresses, components of the Reynolds stress tensor and the different production mechanisms are modulated by the wavy surface undulations.In section 4 we summarize the major findings from this study.

Numerical Methods
In this study, we adopt a customized in-house version of the Incompact3D [34] code framework to perform our DNS study.The dynamical system being solved is the incompressible Navier-Stokes equation for Newtonian flow described in a Cartesian co-ordinate system with x,y,z pointing to streamwise, vertical and spanwise directions respectively.The skew-symmetric vector form of the equations are given by ∇.u = 0. ( Here f represents the body force, p represents the pressure field.The fluid density (ρ) is considered unity for this incompressible fluid as we solve these equations in non-dimensional form.We denote the advection-diffusion term by F for simplicity.Naturally, the above systems of equation can be rewritten to generate a separate equation for pressure.The system of equations are advanced in time using a 3 rd order Adam-Bashforth (AB3) time integration with pressure-velocity coupling using a fractional step method [37].For the channel flow, the body force term, f is dropped.The velocity is staggered by half a cell to the pressure variable for exact conservation of mass.A 6 th Order Central Compact Scheme (6OCCS) with quasi-spectral accuracy is used to calculate the first and second derivative terms (contained in F) in the transport equation.The pressure Poisson equation (PPE) is solved using a spectral method by applying Fast Fourier Transform on the elliptic equation to generate an algebraic equation.the right hand side of the PPE is computed using a quasi-spectral accuracy using 6OCCS and then transformed to Fourier space.To account for the discrepancy between the spectrally accurate derivative for the pressure gradient and a quasi-spectral accuracy for the divergence term, the algorithm uses a modified wavenumber in the pressure solver.
A major downside to the use of higher order schemes as above is the representation of the complex geometry.In particular, the boundary conditions for higher order methods are complex and hard to implement without loss of accuracy near the surface.In this work, we adopt an immersed boundary method (IBM) framework to represent the complex surface shapes.In the IBM the surface representation is accomplished through an added body force term to the governing equations while the background grid can be a simple Cartesian grid.Therefore, this approach saves significant grid generation effort, but is prone to inaccuracies.In this study, we leverage the higher order IBM implementation in Incompact3D using the direct forcing method requiring reconstruction of the velocity field inside the solid region.This can be illustrated using the schematic in figure 1. Figure 1(a) denotes the solid nodes in red, fluid nodes in blue and the interfacial nodes in green.The solid curve represents the continuous shape of the fluid-solid interface.The IBM framework aims to enforce zero velocity at the interface through a velocity field reconstruction in the red solid nodes so that the 6OCCS gradient computations are unaffected.Therefore, the key to the accuracy of this approach is the velocity reconstruction step inside the solid region (red nodes in the schematic) using information at the blue and green nodes.The numerous different IBM implementations [36] differ in the details of this velocity reconstruction.In the current study, we adopt the one-dimensional higher order polynomial reconstruction as reported in [38].This reconstructed velocity field is directly used to estimate the derivatives in the advection and diffusion terms of the transport equation.An illustration of this approach is shown in figure 1(b).Using this 1D polynomial reconstruction, one estimates different solid region velocity fields when computing the derivatives along the different directions (x,y and z).This is an advantage as well as a disadvantage.However, for the purposes of this study, this approach has shown to be reasonably accurate as described in Section 2.3.

Simulation Design
We carry out four different simulations of turbulent flows in flat and wavy channels with different steepness levels (ζ), but with the same peak wave height (a) as shown in figure 2(a).We define an average wave steepness, ζ = 2a/λ, where λ is the wavelength.In our study ζ varies from 0 − 0.022 corresponding to zero, one , one and one half and finally, two waves over the streamwise length of the domain.For all these cases, care was taken to ensure that the friction Reynolds number is maintained to a nearly constant value of ∼ 180.The corresponding bulk Reynolds number, Re b = u b δ ν , for the flat channel case is ∼ 2800.
For the wavy channel turbulent flows with the same effective flow volume and mean channel heights, maintaining the same flow rate (or bulk velocity) increases the corresponding friction Reynolds number, V elocity Figure 1.Illustration of 1D polynomial reconstruction based on Lagrangian polynomial.The solid black curve in (a) represents the fluid-solid interface, red triangular markers represents solid gridpoints where the reconstruction is performed using the fluid gridpoints shown as filled blue circular markers along with the target quantity on the interface marked as filled green circular marker.To retain stability the gridpoints represented as the empty blue circular markers just above the interface is ignored from the reconstruction computation.Dotted black rectangle shows the direction along which the 1D reconstruction is performed as the gridpoints under consideration is enclosed by this rectangle.In (b) a velocity curve is shown which has zero enforced value on the surface (at filled green circular marker).Using the values on the three gridpoints marked as filled blue circle, we extend the curve by computing values on the solid region (at red triangular markers).
boundary conditions are enforced while a uniform grid distribution is adopted.In wall normal direction, no slip condition representing the presence of the solid wall causes inhomogeneity.To capture the viscous layers accurately, a stretched grid is used.The grid stretching in the imhomogeneous direction is carefully chosen using a mapping function that operates well with the spectral solver for the pressure Poisson equation.The different inner scaled grid spacings are also included in Table 1.  1. Tabulation of different design parameters for the simulations such as: wavelength (λ), amplitude (a) and steepness (ζ = 2a λ ) of the wavy surface, friction velocity (u τ ), Reynolds numbers (Re) based on boundary layer height (δ) and different velocities expressed as the subscripts ('cl'=centerline velocity, 'b'=bulk velocity, 'τ'=friction velocity) and the grid spacing in different directions ('∆x'=streamwise, '∆z'=spanwise, '∆y w '=wall normal near the wall, '∆y cl '=wall normal near the flow centerline).Superscript '+' refers to inner scaled quantity (scaled with respect to dynamic viscosity (ν) and friction velocity (u τ )).

Convergence of Turbulence Statistics
In order to quantify the convergence of the simulation and ensure statistical stationarity of the turbulence, we consider the streamwise component of the inner scaled mean spatial and temporally averaged horizontal stress that includes both the mean viscous and Reynolds stress components as x,z,t .Here x,z,t represents the averaging operation with subscripts denoting averaging directions.In the limit of statistically stationary and horizontally homogeneous turbulence, τ H,x (y) can be approximated to a linear profile, 1 − y δ as derived from the mean momentum conservation equations.We estimate a residual convergence error Res as whose variation with y/δ is shown in figure 3.  We note that this error is sufficiently small for the flat channel (ζ = 0) with magnitudes approaching 0.01 near the surface and much smaller in the outer layers.The plot also shows similar quantifications for wavy channel turbulence data with large residual errors near the surface.This is not surprising given that closer to the wall, the turbulence structure is known to deviate from equilibrium due to deviations from horizontal homogeneity.In fact, such deviations from equilibrium phenomenology will be expounded further in the later sections of the article.Nevertheless, we show here that farther away from the surface, the mean horizontal stress approaches equilibrium values as an indicator of stationarity.

Assessment of Simulation Accuracy
We perform a baseline assessment of the computational accuracy for the turbulent channel flow using an immersed flat channel surface before adopting it for more complex surface shapes.We compare the mean and variance profiles from the current DNS of immersed flat channel flow with the well known work of Kim, Moin and Moser [39] (KMM87 here onwards).This turbulent channel flow corresponds to a bulk Reynolds number, Re b ≈ 2800, mean centerline velocity Reynolds number, Re cl ≈ 3300 and a friction Reynolds number, Re τ ≈ 180.KMM87 used nearly 4 × 10 6 (128 × 129 × 128) grid points and solved the flow equations by advancing modified variables, namely, wall-normal vorticity and Laplacian of the wall-normal velocity without explicitly considering pressure.They adopt a Chebychev-tau scheme in the wall-normal direction, Fourier representation in the horizontal and Crank-Nicholson scheme for the time integration.In our work, we adopt a spectrally accurate 6 th order compact scheme in space and a third order Adam-Bashforth time integration as reported in [34].Figure 4 clearly shows that the inner-scaled mean (figure 4(a)) and root mean square of the fluctuations (figure 4(b)) from the current simulations match that of KMM87.We observe slight differences for the streamwise velocity fluctuation RMS in the outer layer which can be attributed to the improved resolution (and accurate time integration) in our simulations.The method employs a staggered grid arrangement for improved mass conservation.x,z,t (KMM87) x,z,t (KMM87) x,z,t (KMM87) x,z,t (Current)

Results
The primary goal of this study is two-fold: (i) to explore the non-equilibrium, near-surface turbulence structure over systematically varied sinusoidal undulations and (ii) characterize the roughness characteristics of such wavy surfaces.In addition, wherever possible, we quantify deviations from equilibrium phenomenology as evidenced in flat channel turbulence, assess the extent of outer layer similarity and relate to characteristic roughness induced effects as appropriate.The analysis can be realized using both instantaneous as well as averaged turbulence structure.In this article we focus on the streamwise-averaged or more commonly known as the 'double-averaged' turbulence structure which is a function of solely the wall normal distance.The term 'double-averaging' refers to the combination of averaging along homogeneous (z, t) and inhomogeneous (x) directions.For the spatial averaging we include both streamwise (x) and spanwise (z) spatial directions and for the temporal (t) averaging we include 2500 three-dimensional snapshots over 20 flow through times for the chosen friction Reynolds number.We use the notation u x,z,t to specify a quantity u being averaged over x, z and t.For the flat surface, horizontal homogeneity and stationarity implies that x, z and t are equivalentin the averaging operation, streamwise (x) and vertical (y) variability.This allows one to characterize the near-surface inhomogeneity along both directions.However, in order to quantify deviations from equilibrium and assess the impact of near-surface inhomogeneity on the turbulence we consider streamwise averaged statistics.

Streamwise Averaging of Turbulence Statistics
In this section, we focus on the deviations from equilibrium in turbulence structure using streamwise averaged statistics that depend only on   x,z,t  1).

Outer Layer Similarity and Mean Velocity Profiles
As the mean channel height (for wavy geometry) is kept constant across all the different steepness, ζ, the observed changes in the mean statistics are only due to surface effects and not the outer layer dynamics.
In figure 6 we show the inner-scaled, double averaged streamwise and vertical velocity for the different cases.The prominent observation for the streamwise velocity is an upward shift (downward shift in the  In particular, the vertical velocity is asymmetric with respect to the symmetric wavy shape as seen from isocontours of time-averaged mean vertical velocity shown in figure 7. We observe that the upward and downward slopes display varying tendencies due to presence of adverse and favorable pressure gradients which decelerate (push the flow upward) and accelerate (push downward) the flow as expected.
However, the extent of upward deceleration dominates the downward acceleration which breaks the symmetry of the flow patterns around the wave crest.This asymmetry increases with ζ resulting in stronger net vertical flow in the lower buffer layer (figure 6(b)).
In spite of these near surface deviations, the dynamics outside the roughness sublayer tend to be similar when normalized and shifted appropriately.anything, the deviation is slightly higher near the surface in the roughness sublayer.

Quantification of Mean Velocity Gradients and Inertial Sublayer
The normalized mean streamwise velocity gradients identify the different regions of the turbulent boundary layer and are especially useful to quantify the extent of the inertial sublayer (or the logarithmic region) and the von Kármán constant.In this study, we estimate the normalized premultiplied inner-scaled mean gradient, γ = y + d u + x,z,t dy + as shown in figure 8 It is easy to see that γ = Φ/κ.We observe that the Φ profiles for different ζ mimic the characteristic equilibrium structure starting from zero at the wall followed by a peak at the edge of viscous layer and subsequently, a gradual decrease in the buffer layer to a value of one in the inertial sublayer.
This clearly indicates outer layer similarity.In fact, there exists an overall shape similarity in Φ hinting at the potential for universality if only the appropriate scales at the different regimes can be identified.
The origin of the 'overshoot' or near-surface peak is well known and is related to the inconsistency from normalization of the mean gradient using inertial scale variables closer to the surface (viscous layer) where the physically relevant characteristic length scale is followed by a decrease as one approaches the inertial sublayer.In the inertial layer, Φ = κy/L il with L il varying linearly with y as per law of the wall (resulting in Φ and γ assuming constant values).
In this context, we see that as the friction velocity, u τ increases with ζ (see Table 1), the viscous length scale, L v decreases resulting in faster growth of Φ = κy/L v in the viscous layer, but over a smaller height that scales with L v .This is consistent with figures 8(a) and 8(b) which show that the magnitude of the peak at the viscous-buffer layer transition decreases with increase in ζ.In addition, we observe an upward (rightward) shift in the log region (i.e.region of nearly constant Φ and γ) with ζ.Taken together, the above observations, namely the upward shift in the log region (figure 8(a)) and the smaller peak in Φ with increase in ζ, indicate that the buffer layer becomes increasingly thicker for steeper waves.The 'buffer layer' is known as a region of high turbulence production [40] where both the viscous and Reynolds stresses are significant.Therefore, the expansion of the buffer layer with ζ is a consequence of the turbulence production zone expanding due to the wavy surface.This is evident from figure 9 where the decay in turbulence kinetic energy (TKE) production is slower for higher ζ in the buffer region (y + ≈ 10 − 50) in both inner-scaled (figure 9(a)) and dimensional (figure 9(b)) forms.We expect this trend to be even stronger in the presence of significant separation at larger values of ζ.This trend is consistent with prevalent understanding of classical rough wall boundary layers at high Reynolds numbers, especially in the lower atmosphere where the roughness elements of size a + 50 − 100 tend to completely destroy the viscous layer [16] if not most of the buffer layer.In our studies, a + ≈ 13 for the different ζ (see Table 1) and only modulates the buffer layer.A related observation is that the vertical location of the inner scaled peak turbulence production (y + ≈ 12) does not change with ζ, but the magnitude decreases.This is not surprising as for ζ > 0, there exists other sources of turbulence generation, i.e. from the surface roughness or undulations which contributes to the total friction.

Characterization of the Roughness Function and Roughness Scales
A common way to assess the influence of the wavy surface on turbulence structure is to quantify the effective drag and its influences on the flow structure.While the increase in friction velocity for a fixed Re b (apparent from Table 1) is a natural way to quantify the increased drag, estimating the downward shift in the mean streamwise velocity profile (figure 6(a)) is another approach and often used to characterize the effective roughness scales.It is well known that the logarithmic region in the equilibrium flat channel turbulent boundary layer (TBL) is given by where the additive constant B is typically estimated to fall within the range, ≈ 5.0 − 6.0 and depends on the details of the buffer and viscous layer for a given simulation or measurement.The flat channel data in the current work provides an estimate of ≈ 5.6, possibly due to a combination of the friction Reynolds number regime and simulation algorithm.In the presence of surface undulations of scale a, we observed from the earlier discussion that the log region underwent a upward shift due to an expanding buffer layer.
As per [13,16,25], the influence of these buffer layer modulations on the log layer shift is characterized in terms of a modified logarithmic profile for rough-wall turbulent boundary layers (TBLs), where ∆ u + x,z,t is defined as the roughness function.The roughness function, ∆ u + x,z,t can be related to the characteristic "equivalent" sand grain roughness, k s as and the characteristic roughness length, k 0 as It is easily seen that k 0 = k s e −8.5κ .While k s and k 0 are used to quantify the non-equilibrium 'roughness' effects near the surface, they mostly cater to complex roughness such as grasslands, urban canopies or sand grain type surfaces.Of course, k s corresponds to a case where the buffer layer dynamics is significantly modified by the roughness while k 0 corresponds to a situation where the buffer and.viscous layers are completely destroyed by the roughness.Therefore, such metrics do not represent the smooth, low steepness surfaces adopted in this work.For comparison sake, we also report the equivalent sand grain roughness, k s of Nikuradse [3]   Given that a + is nearly constant for all the cases while k + s and k + 0 show near linear growth confirm that such wavy surfaces do not fit the k-type roughness description [16,41].In addition, given that the boundary layer height in channel flow turbulence is fixed as a constant δ, the independence of k s , k 0 on δ indicate a departure from d-type classification [41].In fact such surfaces as considered in this work belong to the 'transitional' and 'waviness' regime as a + ∼ 13 does not represent a sufficiently large (i.e O(100)) roughness Reynolds number, Re a = au τ /ν.We have also reported the λ + values for the different cases in Tables 1 and 2.
The 'waviness' regime implies a surface that is very different from a Nikuradse roughness dominated by form drag caused by flow separation and vortical re-circulation zones within the roughness sublayer.
Therefore, strong waviness causes the drag (as estimated by the roughness function ∆ u + x,z,t ) to be smaller than the corresponding Nikuradse value for a given k + = a + .At low slopes (waviness) the overall drag is more dominated by viscous shear and less by the form drag.The opposite is true in the roughness regime.This is clearly illustrated in figure 11(a) where the correlation between ∆ u + x,z,t and k + = a + from Nikuradse [3] and Colebrook [4]   and Reynolds number (Re τ ).This transition has been correlated to the dominance of form drag over viscous drag [12,17].For the benefit of the reader, we have explicitly documented the roughness function correlations of Nikuradse and Colebrook used above in Appendix A.   Peer-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161

Characterization of Horizontal Flow Stress and Implications to Drag
The horizontal flow stress directly impacts the flow drag through the boundary layer and in turn the mean velocity profiles discussed above.The viscous flow stress τ V acting on a fluid particle is described including both spanwise and streamwise components as without overbars denoting their magnitudes.The decrease in magnitude of the inner-scaled viscous stress with ζ in the viscous and lower buffer layers is a consequence of the normalization using the averaged wall stress, u 2 τ which increases with wave steepness.We observe that the mean τ V is relatively unaffected, but its contribution to the total drag decreases with increase in ζ.In general, the mean streamwise flow near the wall slows down due to the presence of wave-like undulations (see figure 6(a)) which in turn reduces its gradient in the wall normal direction.This reduction in the average viscous stress is compensated by the non-zero vertical velocity and its variation along the streamwise and vertical direction.This explains why the net double-averaged (i.e. both temporally and spatially averaged) dimensional viscous stress sees very little increase in the viscous layer as seen from the dimensional stress profiles in figure 12(e).This observation clearly indicates that the increase in net wall stress (u 2 τ ) with ζ has its origins in the increase of Reynolds stress in the buffer-log layer transition as seen in figure 12(f) which is reflected in the total mean stress variation as well (figure 12(d)).Given that the non-dimensional roughness scale, a + ≈ 13 corresponds to the buffer layer, it is not surprising that the buffer layer shoulders much of the effect of increasing wave steepness.
However, the mechanism underlying increase in the peak double-averaged Reynolds stress with ζ will invariably depend on the structure of the attached (or detached) shear layers in the vicinity of the wavy surface resulting in a coupling between the viscous shear layers and buffer layer turbulence production.
The nature of this coupling will be further explored in the future.Of course, when the shears layers are detached as in a separated flow, the interactions could entail very different characteristics.

Characterization of Reynolds Stress Tensor and its Production
In the earlier discussions, we focused on the mean gradients and their impact on the horizontal flow stresses.In this section, we focus on the effect of changing ζ on elements of the Reynolds stress tensor and the turbulent kinetic energy that are borne out of the interaction between mean gradients and the Reynolds stress.We observed earlier (figure 8) that the peak in the mean gradients at the start of the buffer layer decreases with surface undulations which also impacts turbulence production (figure 9) in the lower buffer layer and in turn the individual components of the Reynolds stress tensor, u i u j x,z,t .In fact, the most noticeable deviations from equilibrium in wavy wall turbulence occur in the second order statistics.In particular, we observe in figure 13(a) the inner-scaled streamwise variance that peaks in the buffer layer and this peak value decreases with increase in ζ.In addition, the inner scaled profiles nearly collapse in the outer region for all ζ while the location of peak streamwise variance shifts upward as ζ increases.Given the lack of signifcant flow separation, this upward shift in the peak variance is modest, but noticeable.In figure 13 we identify the peak location for each curve with color matched horizontal lines so that the trends can be identified.Using this, we see a systematic upward shift of the peak value of u 2 + x,z,t with ζ in figure 13(a).Related research by Ganju et al. [42] has shown that this upward shift is tied to significant increases in roughness scale, a + which can cause very different dynamics around the wave including flow separation.Further, the effect of changing λ was reported to be minimal in their investigation.Their first observation is consistent with classical understanding of high Reynolds number rough wall turbulence [16,40] where the peak variances occur at nearly the roughness height, a.In fact, wall stress boundary conditions for large eddy simulation over rough surfaces are designed to model the same.In our study, the amplitude, a is fixed while the wavelength, λ is decreased in order to change ζ = 2a/λ.For a fixed bulk Reynolds number, the decrease in λ (or increase in ζ) increases the net drag and in turn the friction velocity, u τ .The resulting decrease in the viscous length scale, L v = ν u τ changes the inner-scaled wave height, a + = a/L v rather modestly from 13.07 − 13.92 when ζ doubles from 0.011 to 0.022 (see Table 2).Therefore, this systematic upward shift in the location of peak streamwise variance cannot be solely attributed to these very modest increases in a + .In fact, the wave steepness significantly impacts the buffer layer dynamics and in turn the variance distribution through the turbulence production mechanisms as delineated under.
We further dissect the above observations using the variance production (figure 14 ) with significant production in the roughness sublayer (y + a + ) followed by destruction above the roughness scale, a + that decays with height.This production and destruction process clearly represents deviation from equilibrium as its origins lie in the streamwise mean velocity gradient, d u z,t /dx being non-zero from horizontal inhomogeneity.
A key consequence of this inhomogeneity driven destruction process is that both the peak inner-scaled variance production and the peak variance decrease with ζ (figure 13(a)).However, it also turns out that the systematic upward shift observed for the peak in P u v  x,z,t (zoomed near the surface) and (c) vertical gradient of streamwise velocity, d u + x,z,t /dy + .The black horizontal line corresponds to the average of the maximum magnitude of u v + x,z,t for the different ζ.Note that the individual peak values were too close to each other to be shown separately.
In the absence of three-dimensional flow (both forced by a three-dimensional surface and induced by separation over a two-dimensional surface), the location of peak w 2 + x,z,t shows no clear monotonic trend although a consistent downward shift is observed for ζ > 0. This can be attributed to either the small amounts of separation observed in these flows (see figure 5) or due to conversion of the vertical variance produced from the surface undulation through the pressure-strain term.We expect the latter to be the likely mechanism although no quantification is provided work to support this hypothesis.As one would expect away from the surface, the v 2 + x,z,t and w 2 + x,z,t profiles across different values of ζ approach each other in the outer layer indicating that the effect of the surface undulations is concentrated closer to the surface.We nevertheless note that in this region of the TBL, v The mean inner scaled turbulent kinetic energy, TKE + , displays the cumulative effect of the individual variances as shown in figure 13(d).In particular, we observe an exaggerated upward shift (note the horizontal lines in figure 13(d)) in the location of peak k + in the buffer layer.This is caused by the combined effects of the upward shift in u 2 + x,z,t along with the downward shifts in both v 2 + x,z,t and w 2 + x,z,t .
Beyond the peak, the different curves nearly collapse in the outer layer although in the inertial logarithmic region, TKE + shows consistently higher values for the wavy turbulence cases due to systematically higher v 2 and w 2 for ζ > 0.

Vertical Turbulent Momentum Flux
In addition to the diagonal components of the Reynolds stress tensor discussed above, we also look at the dominant off-diagonal terms, namely, − u v + x,z,t and − w v + x,z,t as shown in figure 15.In the limit of high Reynolds number (i.e Re τ ≥ 4000), for channel flow turbulence over smooth flat surfaces, there exists a well defined log layer with nearly constant u v + [44].This nearly constant u v + layer is very narrow at low Reynolds numbers.Nevertheless, this study still lets us characterize the influence of the wavelike undulation on the turbulence structure.As discussed earlier, the location of the peak in u v + x,z,t (at y + ≈ 32) shows very little variation with no clear trend, but its magnitude increases with ζ.
The peak location in u v + x,z,t falls roughly between the peak values of u 2 + x,z,t and v 2 + x,z,t as shown in figures 13(a) and 13(b).This increase in peak value is not surprising as the wavy surfaces naturally increase in the positive peak of − u v + x,z,t is correlated with a small negative peak in the viscous layer at higher values of ζ, indicative of the surface shape induced mixing that is different from turbulence induced stress.This negative value of − u v + x,z,t close to the surface indicate that higher momentum fluid particles move away from the wall due to up slope part of the wavy surface.

Conclusion
In this work, we report outcomes from a DNS-based investigation of the turbulence structure and its deviation from equilibrium using high Reynolds number flow between two infinitely parallel plates with 2D wavy undulations.In particular, we set out to assess the influence of small wave slopes (with very little flow separation) on the turbulence structure and their correspondence to common roughness characterization that invariably deals with the high slope regime.To maximize the shape sensitivity on the flow structure, we operate in a transitional roughness Reynolds number (k + = a + ∼ 13) which is much smaller than the fully rough regime corresponding to (k + = a + 70).
The streamwise mean velocity structure indicates a characteristic downward shift (to higher y + ) of the logarithmic region of the TBL indicating increased flow drag with increase in wave slope, ζ.This is associated with a sustained upward vertical flow in the lower roughness sublayer and corresponding downward flow in the buffer layer.The strength of these vertical motions increases with ζ and have a dominant role to play in the near surface turbulence production processes.In fact, analysis of the mean non-dimensional streamwise velocity gradients and inner-scaled turbulence production show that the buffer layer expands with increasing wave slope, ζ indicating that the well-known equilibrium understanding of near-wall turbulence processes is modulated even for such highly shallow wavy surfaces.
In fact, characterization of the roughness effects from such shallow surface undulations with minimal flow separation is very different from the separation and form drag dominated Nikuradse type sand grain rough surfaces.Therefore, in our case, drag as represented by the roughness function turns out to be very weakly dependent on the wave amplitude (related to roughness height, k + = a + ) and more on the effective wave slope (2ζ).These conclusions are consistent with [12,17] where wavy surfaces in the high slope limit approach Nikuradse type roughness.Therefore, in such cases, one needs to model ∆ u + as f (a + , ζ).
In the absence of flow separation, the turbulence generation process within the roughness sublayer is primarily topology driven.For the two-dimensional surface undulations considered in this work, the differences in turbulence generation for different ζ originate in the production of small quantities of vertical velocity variance closer to the surface due to non-zero values for the vertical velocity gradients.In contrast, for the equilibrium TBL over a flat surface, horizontal homogeneity implies that the mean vertical velocity and its gradients are zero.However, the predominant generation of vertical variance still occurs through redistribution of the streamwise velocity variance (generated closer to the surface) through the return to isotropy term.Therefore, these two different production mechanisms interact to cause v 2 + x,z,t to peak at a larger value and closer to the surface (relative to a flat TBL) with increasing wave steepness, ζ.Similar trends are observed for the inner-scaled spanwise velocity variance.
The effect on the streamwise velocity variance, u 2 + x,z,t is more complicated.Specifically, the inner-scaled streamwise variance, u 2 + x,z,t (and consequently, TKE + ) displays a distinct upward shift in the location of this peak value for increasing wave steepness.In addition, the normalized peak variance magnitude decreases with ζ.Our analysis shows that even though surface-driven changes impact the various production terms ( P u u 11 + x,z,t and P u v 11 + x,z,t ) for u Overall, in the absence of significant flow separation, the surface undulations seem to generate vertical turbulence fluctuations closer to the surface which in turn modulates the entire near-wall turbulence structure.In the presence of steeper 2D waves with significant flow separation or three-dimensional wavy surfaces, the complexity will be enhanced due to the generation of spanwise fluctuations near the surface.

Figure 2 .
Figure 2. Schematic illustration of the cartesian grid with the immersed boundaries of different shapes in (a) and a close-up of the buffer region in (b).The solid thick curve represents the wave for λ = 4π and the dashed line for λ = 8π 3 .A similar setup is used for other surface shapes as well.

Figure 3 .
Figure 3. Quantification of statistical stationarity for the different DNS datasets using the residual of mean horizontal stress from 2500 samples collected over ∼ 12 δ u τ .
a) Normalized and averaged streamwise velocity distribution in wall coordinates

Figure 4 .
Figure 4. Comparison of mean velocity and RMS velocity fluctuation between DNS of flat channel turbulent flow with IBM and the Kim et al. [39] DNS without IBM

Figure 5 .
Figure 5.Comparison of instantaneous flow separation for the different wave steepness, ζ.The wavy surface is denoted in cyan and the separation is denoted in red.

yδ
and ζ.To average along the wavy surface, we define a local vertical coordinate, y local,1 at each streamwise location with y local,1 = 0 at the wall.Its maximum possible value is the mid channel height and changes with streamwise location.We then perform streamwise averaging along constant values of y local,1

Figure 6 .
Figure 6.Inner scaled mean (a) streamwise velocity (b) vertical velocity and (c) defect velocity computed using local coordinate-based average.The thick lines represent averaging at constant y local,1 and the thin lines with markers represent averaging at scaled y local,2 .Three vertical straight lines correspond to the different a + for ζ > 0 (see Table1).

Figure 7 .
Figure 7. Spanwise and temporally averaged streamwise and vertical velocity over wavy surfaces in turbulent channel flow.

Figure 8 .
Figure 8. Variation of non-dimensional mean streamwise velocity gradients, (a) γ = y + d u + x,z,t dy + and (b) Φ = κy u τ d u x,z,t dy .The thin dashed black line in (a) corresponds to the mean γ valued 2.5604 computed based on y + = 60 − 110.A related quantification often employed to interpret near wall structure is the non-dimensional mean streamwise velocity gradient, Φ = κy u τ d u x,z,t dy whose variation with inner-scaled wall normal distance is shown in figure 8(b).It is easy to see that γ = Φ/κ.We observe that the Φ profiles for different ζ mimic the

Figure 9 .
Figure 9. Schematic illustrating the wall-normal variation of streamwise averaged production of turbulent kinetic energy in (a) inner variable non-dimensionalized and (b) dimensional (m 2 /s 3 ) forms.
and the characteristic roughness length, k 0 scaled by the inner-layer variables for different values of ζ.As expected, these different roughness metrics increase linearly with wave steepness as seen in figures 10(a)-10(c).The effective sand grain roughness assumes a non-zero value of ≈ 3.3 for ζ = 0 due to the upward shift caused by the viscous and buffer layers.Therefore, k + s ≈ 4 is indicative of a nearly smooth wall which in our study corresponds to ζ ∼ 0 − 0.01.The higher values of ζ considered in this work generate k + s ∼ 6 although no substantial flow separation is observed.As expected, the k + 0 is extremely small indicating that the flow is smooth enough to retain the viscous and buffer layers.Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1Peer-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161Roughness function variation in the log region.

Figure 10 .
Figure 10.Variation of the different roughness quantifications with ζ in (a), (b), (c) and wall normal variation of mean roughness function in (d).

Figure 11 (
b), also shows the data from Schultz and Flack[12] who performed experiments with systematically varied pyramid roughness elements of different slope.These data trends indicate that the slope transition from waviness to Nikuradse type roughness regime (denoted by a vertical line in figure11(b)) occurs between ES ∼ 0.25 − 0.4 (ζ ∼ 0.12 − 0.18) with possible dependence on the extent of separation between the surface and outer layer scales (δ/k) Mean roughness function variation with roughnessReynolds number in comparison with the correlation obtained from the data of Nikuradse[3], the Colebrook[4] relation and fully rough asymptote along with the current DNS (See appendix for the correlation equations).Napoli et al. (2008) Schultz et al. (2009) Current DNS(b) Comparison of variation of mean roughness function with effective slope (ES) reported in Napoli et al.[17] and Schultz et al.[12] with the current DNS.

Figure 11 .
Figure 11.Variation of mean roughness function (a) with roughness reynolds number and (b) with effective slope in comparison with reported data from known literature.In summary, the modulations in the mean averaged first order statistics from wavy surface undulations manifest as: (i) increase in drag (through friction velocity u τ ); (ii) modified buffer region including (iii) a systematic upward flow in the buffer layer and a smaller downward flow at the lower logarithmic layer.To interpret the above effects better, we analyze the horizontal flow stress and components of the Reynolds stress tensor in the following sections.

Figure 12 (
Figure 12(a) shows the inner-scaled double-averaged horizontal stress magnitude, τ H felt by a fluid particle.We further split this into the inner-scaled viscous and turbulent parts, τ V and τ V respectively as shown in figures 12(b) and 12(c).In the viscous layer, the total stress is dominated by the viscous stress for the different cases A-D with different ζ varying between 0 − 0.022.The inner scaled mean viscous shear stress magnitude (figure 12(b)) decreases with steepness (figure 12(b)) in the viscous layer where it is nearly constant before decreasing across the buffer layer.Away from the mean surface level, in the buffer layer, the inner-scaled Reynolds stress magnitude grows (from near-zero values in the viscous layer) into a peak value at y + ≈ 35 (figure 12(c)) whose magnitude increases with ζ before collapsing over each other in the log layer.Overall, the viscous stress dominates in the viscous layer while the Reynolds stress grows through the buffer layer (a region where the viscous stresses continually decrease in importance) to peak at the buffer-log transition region.

Figure 12 .
Figure 12.The schematic shows the inner scaled mean (a) horizontal stress, (b) viscous stress and (c) Reynolds stress in the top row and the dimensional mean (d) horizontal stress, (e) viscous stress and (f) Reynolds stress in the bottom row.The vertical lines correspond to the different a + values.

Figure 13 .Figure 14 .
Figure 13.Inner scaled mean (a) streamwise variance, (b) vertical variance, (c) spanwise variance and (d) turbulent kinetic energy (TKE).The horizontal lines correspond to height with maximum value of the statistics along the profile.

Figure 15 .
Figure 15.Inner scaled mean (a) covariance u v +x,z,t , (b) covariance u v + x,z,t (zoomed near the surface) and (c) vertical gradient of streamwise velocity, d u + x,z,t /dy + .The black horizontal line corresponds to the average of the maximum magnitude of u v +x,z,t for the different ζ.Note that the individual peak values were too close to each other to be shown separately.

generate (figures 14
(a) and 14(d)) stronger u and v fluctuations.As shown in figures 15(a) and 15(b),

2 x
,z,t , the overall variance production, P 11 + x,z,t , does not show an upward shift in the peak.Therefore, plausible mechanisms for generation of this upward shift still point to the vertical variance being modulated through the pressure-strain term.
e. ES = 4a/λ = 2ζ where ζ is the steepness factor.In this research, ζ is deliberately varied from 0 to 0.022 (ES ∼ 0 − 0.044) which is nearly an order of magnitude smaller than the transition location (in ES) beyond which ∆ u + becomes constant for a given mean roughness height a + .For this range of slope parameter (ζ), there exists very little flow separation over the 2D surface and consequently, very little spanwise flow.The [3]ghness height or wave amplitude, a, is chosen to generate moderate scale separation, i.e. δ/a ≈ 15 and the ratio based on the equivalent sand roughness height (assuming the Nikuradse[3]form) turns out to be δ/k s ≈ 30 − 50.These values are normally sufficient to generate outer layer similarity based

periodic Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1
due to increase in u τ with wave steepness, ζ.However, this increment is as at most ∼ 10% in the current work and is not expected to influence our analysis significantly.The simulation parameters for the different cases are summarized in Table1.The simulation domain is chosen as 4πδ × 2.2δ × 4πδ/3 (including the buffer zone for the IBM) where δ is the boundary layer height.This volume is discretized using a resolution of 256 × 257 × 168 grid points.In the streamwise and spanwise directions, , to generate mean statistical profiles.A slight variant of the above is to use a rescaled coordinate y local,2 = y local,1 × δ δ local which stretches y local,2 everywhere between [0, δ] where δ is the mean half channel height.One averages over constant values of y local,2 to generate another set of double-averaged mean statistics.Both these approaches implicitly approximate the terrain as nearly flat with a large radius of curvature in a local sense and therefore, nearly homogeneous.This approximation works well when a δ << 1.In our study a δ = 0.07 which is an order of magnitude larger than the typical viscous length scale, L v = ν/u τ = 1/Re τ ≈ 0.0055, but smaller than the log layer (y + ≈ 50) with strong inertial dynamics.For the mean velocity results presented in this section, we compare both the averaging approaches to illustrate their closeness to each other.Specifically, we use thick solid lines to denote the mean profiles averaged over constant local coordinate, y local,1 and thin lines with markers to denote averaged quantities using scaled local coordinate, y local,2 .However, for the rest of our analysis, we average over y local,1 in a manner consistent with the literature.The different colors, namely, blue, green, red and lime are associated with different wave steepness, ζ = 0, ζ = 0.011, ζ = 0.017 and ζ = 0.022 respectively.

16 July 2019 doi:10.20944/preprints201907.0192.v1
Peer-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161 is indicative of slowing down of the flow near the wall from increased drag due to the wavy surface for a fixed mass flow rate (bulk Reynolds number).This would naturally result in higher centerline velocities and Re cl as seen in Table1in order to maintain the prescribed flow rate.The vertical mean velocity structure (figure6(b)) is consistent with this interpretation as the wavy undulations generate increasingly with increasing wave steepness, ζ (figure 6(a)) which increases with y + before showing near linear growth in the log-layer.This trend is well known for rough-wall turbulent boundary layers [16] and Preprints (www.preprints.org)| NOT PEER-REVIEWED | Posted:

preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1 Peer
With some analysis, one can easily show that Φ undergoes a linear growth as Φ = κy/L v near the surface (L v being a constant).In the buffer layer, one can similarly formulate Φ = κy/L bl with L bl increasing super linearly with y to cause the peak Preprints (www.-reviewedversion available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161 Table 2 compiles estimates of the roughness function, ∆ u + x,z,t y , averaged over the entire logarithmic region given by y + ≈ 60 − 120 as illustrated in figure 10(d) over which the values are nearly constant.For all the metrics reported in this work, we use data from the averaged profiles across constant values of the non-scaled local coordinate, y local,1 .

preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1
[3]]-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161curve.The wave slope dependence on the flow drag is evident from figure11(b) where ∆ u +x,z,t is shown against the effective slope, ES= 2ζ.We clearly see that our data follows the trend of Napoli et al.[17], i.e., ∆ u + x,z,t increases with ES until it asymptotes to a value dependent on k + = a + and the flow Reynolds number.This capping value can be estimated somewhat accurately from the Nikuradse curve[3]for sand grain roughness (this value is denoted by the horizontal line in figure11(b) for our DNS data) as the sand grains typically represent a high effective slope surface.
are compared with our current DNS data.We clearly see that for the current study with nearly constant a + , ∆ u +x,z,t increases with wave steepness ζ to approach the Nikuradse Preprints (www.

Table 2 .
Tabulation of estimated turbulence parameters, namely, von Kármán constants for the different cases and commonly used roughness parameters.
Similarly, the Reynolds stress is given by τ R = τ R xy î + τ R zy k, where τ R xy = − u v x,z,t and τ R zy = − w v x,z,t .The total horizontal stress is then τ H ∂z ).

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1
Peer-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161 /dy + respectively.It is to be noted that u v + x,z,t peaks at the edge of the buffer layer at y + ≈ 32 whereas the normalized mean gradient, achieves its maximum value near the surface.In addition, the location of the peak in u v +x,z,t (at y + ≈ 32) shows very little variation with no clear trend, but its magnitude increases with ζ all through the buffer and log layers.Contrary to this, the magnitude of d u + z,t /dy + decreases with ζ due to surface undulations whose influence decreases away from the surface (through the viscous and buffer layers).In summary, we understand that the combined influence of the

Preprints (www.preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1
[42]-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161(thepeakoccursatoraroundy+≈12)shown in figure14(a).In summary, we observe that more severe the surface inhomogeneities, smaller the rate of streamwise turbulence production through d u + z,t /dy + , d u + z,t /dx + and u v + x,z,t , but the location of peak production in wall coordinates remain unaffected.The effect of surface undulations on vertical variance profiles is opposite to that observed for the streamwise variance, i.e. the peak variance in the buffer-log transition region (y + ≈ 55) increase and shift downward with ζ in comparison to flat channel turbulence data.While the shift in the peak location is systematic, we note that the peak value in itself changes very little from ζ = 0.011 to 0.022 after making a sharp jump from ζ = 0.0 to 0.011.It is well known[40]that streamwise turbulent fluctuations are generated closer to the surface in the (buffer) layer and is then distributed to the other components through the pressure-strain term[43].Consequently, the vertical and spanwise variances peak further away from the surface, closer to the log layer.Among the two, the vertical variance is nominally expected to achieve maximum value further away from the surface due to the wall effect, i.e. vertical fluctuations are damped closer to the wall.In our investigations, the vertical variance peaks closer to the log layer (y + ≈ 55) while the spanwise variance peaks lower in the buffer-log transition (y = 0 − 0.022.In essence, what we are observing is that changes in λ with fixed a impacts the growth rate in the buffer layer with perceptible effect on the peak location and magnitude.This trend is maintained until strong separation related dynamics including detachment of the shear layer sets in at higher ζ.As shown in Ganju et al.[42], this causes a secondary peak in the buffer layer for both vertical and spanwise variance, possibly due to turbulence production within the separation bubble as well as above it.Such effects are absent in the current study as evidenced by just a single peak for the + ≈ 35) as observed in figures 13(b) and 13(c) respectively.Due to the presence of surface undulations (non-zero ζ), vertical velocity variance is produced closer to wall (or effective wall for wavy surfaces) in the roughness sublayer (i.e y + a + ) as seen from the inner-scaled vertical variance production, P 22 + x,z,t in figure 14(d).The extent of near surface production increases with wave steepness, ζ.To dissect further, P 22 + x,z,t is further split into + to grow faster (figure 13(b)) through the viscous and buffer layers to ultimately peak closer to the surface (in the buffer-log transition).The location of peak vertical variance is illustrated in figure 13(b) through coloured horizontal lines that show a clear downward shift from ζ figure 13(c)) through the viscous and lower regions of the buffer layer to ultimately peak in the buffer layer (y + ≈ 35).Intriguingly, we note a systematic shift in the peak location and magnitude only for ζ = 0.0 to 0.011, but little variation from ζ = 0.011 to 0.022.This may be attributed to the peak occurring farther away from the surface where the inhomegeneity effects are small.As we are dealing Preprints (www.

preprints.org) | NOT PEER-REVIEWED | Posted: 16 July 2019 doi:10.20944/preprints201907.0192.v1
Peer-reviewed version available at Fluids 2019, 4, 161; doi:10.3390/fluids4030161with mildly steep two-dimensional wavy surfaces in this study, we observe that w 2 x,z,t is not produced near the surface as evidenced by the production terms in the variance transport equation (not shown here) being nearly zero throughout the boundary layer due to d w z,t /dx = d w z,t /dy = 0.In spite of the quasi-two-dimensional wavy surfaces employed here, the above trends will breakdown in the presence of strong separation that can introduce three-dimensional flow patterns.As part of an ongoing research study we are exploring higher values of ζ to verify the above statement.