Wavy Walls, a Passive Way to Control the Transition to Turbulence. Detailed Simulation and Physical Explanation

: Inducing spanwise motions in the vicinity of solid boundaries alters the energy, mass and/or momentum transfer. Under some conditions, these motions are such that drag is reduced and/or transition to turbulence is delayed. There are several possibilities to induce those spanwise motions, be it through active imposition a predeﬁned velocity distribution at the walls or by careful design of the wall shape, which corresponds to passive control.In this contribution, we investigate the effect that wavy walls might have on delaying transition to turbulence. Direct Numerical Simulation of both planar and wavy-walled channel ﬂows at laminar and turbulent regimes are conducted. A pseudo laminar regime that remains stable until a Reynolds number 20% higher that the critical is found for the wavy-walled simulations. Dynamic Mode Decomposition applied to the simulation data reveals that in these conﬁgurations, modes with wavelength and frequency compatible with the surface undulation pattern appear. We explain and visualize the appearance of these modes. At higher Reynolds numbers we show that these modes remain present but are not dominant anymore. This work is an initial demonstration that ﬂow control strategies that trigger underlying stable modes can keep or conduct the ﬂow to new conﬁgurations more stable than the original one.


Introduction
Flows over wavy boundaries are of great relevance in both geophysics and engineering [1,2], e.g., the dynamics of sand dunes in deserts and sediment dunes in rivers or the interaction of the oceanic wind and surface waves and their impact on energy-generating devices. Flows over an undulated boundary are more involved than those past flat surfaces, but at the same time offer richer phenomena: for a specified Reynolds number, an adequate choice of the amplitude and wavelength of the undulation may lead to either an enhancement or an inhibition of heat, mass and/or momentum transfer [3].
Several authors have already addressed flows over wavy configurations; in the incompressible [1,4] and in the compressible [5][6][7] regimes. In either case, the induced pressure gradient created by the wall waviness modifies the flow behavior. Larger amplitude undulations lead to scenarios with multiple separations and reattachments. The intermittent nature of the locally separated area beyond a surface peak affects the turbulent production; under certain conditions the coherent structures can be inhibited, and streaks might even disappear then.
The relative angle between the wavy pattern and the bulk velocity also plays a role on the flow behavior. Chernyshenko [8] hinted that a fluid moving obliquely over a sinusoidal undulated surface might experience a reduction in drag, which would derive from the spanwise motions in the wall vicinity induced by the traverse pressure gradient. Indeed, the drag reducing effect of spanwise/streamwise periodic velocity enforced at the wall had been already investigated using Direct Numerical Simulation (DNS): [9][10][11] showed how the spatial Stokes layer disrupts the turbulent production cycle, resulting into a positive shift of the logarithm layer (much as riblets do, [12]). They went further as to identify combinations of forcing amplitudes, wavelengths and frequencies for the wall velocity leading to net-energy savings of up to 23%. However, this impressive figure is obtained for velocity perturbations of the order of the bulk velocity, hence difficulty in the practical implementation of this active control approach. This motivated the interest in wavy walls, a fully passive means to achieve drag reduction. Departing from DNS data [10], Chernyshenko derived a first estimate of the angle and the wavelengths involved. The study was completed later using DNS [3] considering different Reynolds numbers, wave lengths and amplitudes (a sketch of the configuration is shown in Figure 1). The spanwise shear-strain profiles obtained were shown to resemble closely those associated with in-plane wall motions. Unfortunately, since the spanwise velocities were moderate, only a limited drag reduction efficiency was attained: a maximum net drag reduction of about 0.7% was reported at Re τ = 360 for a wavy wall with amplitude 20 + , wavelength 920 + and flow angle θ = 70 • , where the + represents wall units. However, the authors mention that: "although net drag reduction levels of 1-2% were observed for the three Reynolds number tested (180, 360 and 1000), the value is not quantifiable with any greater degree of precision".
Despite the mild drag reduction attained, the facts that no power supply is required, and the easier manufacturing and posterior maintenance make of wavy walls an interesting technique. However, the majority of the studies here referenced are in turbulent conditions and to the author knowledge, the effect of the waviness in transition to turbulence has not been reported yet in detail. As noted in [13], well-designed perturbations introduced in a boundary layer can potentially inhibit the growth of the most unstable disturbances. This has been recently confirmed by [14], who experimentally described how spanwise perturbation introduced by a 3D wavy-walled roughness is more effective in delaying transition than its 2D counterpart (transition Re = 120,000 vs. 50,000). This experiment is performed in boundary layer flow and a not detailed explanation of the reason of this delay is concluded.
To shed some light about the physical reasons, in this contribution, we investigate the effect of wavy walls in delaying the transition to turbulence. To do that, we conduct DNS of planar and wavy-walled channels, and compare the results. Dynamic Mode Decomposition (DMD, [15]) is applied then to the DNS channel data [16,17], in an attempt to discern the effect of the wall undulation on the flow field. The document is structured as follows: next section describes the methodology employed. Section 3 describes the different computations conducted and offers a first analysis. Section 4 investigates the impact that wavy wall has on transition delay. The conclusions are stated then in the last section.  [3]. Most computations conducted considering a slanted inlet velocity profile (b) with θ = 70 • ; one verification computation has been conducted on a slanted domain (c). Additionally, in (b), one of the spectral domains used for discretization is shown.

Description of the Flow Solver
All the flow computations conducted in this work used our in-house High Order Discontinuous Galerkin Spectral Element solver, extensively described and validated in [18][19][20][21]. Direct Numerical Simulation of transitional and turbulent flows requires tracking accurately a wide range of spatio-temporal scales: high order spectral element methods, given its low diffusion and dispersion errors, are excellent candidates for this task.
Our solver employs a multidomain approach to discretize the compressible Navier-Stokes equations: with U denoting the vector of conserved variables (density, momentum and energy). The flow domain Ω (Figure 1) is partitioned into non-overlapping hexahedral elements Ω e (one of them shown in Figure 1). Inside every element, the unknown U is expanded on a spectral basis constructed as the tensor product of Gauss-Legendre polynomials of degree p x j . The different subdomains are connected in a weak sense by a Roe Riemann solver.
For the channel flow configurations considered in this work, periodic boundary conditions are considered along the streamwise and spanwise homogeneous directions x-y. A constant unitary mass flow rate is enforced through a homogeneous, time-dependent volumetric force/source term S that is imposed at the Navier-Stokes equations: where k u = cos(θ) − ρu k and k v = sin(θ) − ρv k penalize the instantaneous imbalance between the computed space average flow rate (ρu k , ρv k ) and the desired (cos(θ), sin(θ)). The constant α is taken as 0.9 in all the simulations considered in this work. This source term is a mathematical artefact (although equivalent to a physical force or external pressure gradient) that forces the "average" streamwise and spanwise momentum components to be, respectively, (cos(θ), sin(θ)). It is important to highlight that although the integration of the Navier-Stokes in time and space is accurate with our methods, at the initial stages of the computation, the average momentum is far from the prescribed, and that difference can be positive or negative. This explains why the source terms (Equation (1)) could oscillate until the flow achieves its statistically steady state and the momentum converges (see e.g., Figure 2, more details in Section 2.1).
The source term Equation (1) driving the flow is intimately related to the friction. In that sense, it will be used has a direct sensor to detect transition to turbulence, since drag undergoes a large increase when the flow turns turbulent. For the planar channel, this relation is direct, whereas for the wavy wall the total drag force splits into friction and pressure contributions. Friction and pressure may be further segregated through integration over the wavy surface and projection along the flow-direction, yielding the drag coefficients D f and D p , respectively. In this work, however, we are interested only in global drag, regardless the origin. Accordingly, we resort to the temporal average of the source term to compare the different configurations here studied: with || U b || the average bulk velocity.  (2)). Observe that during the initial steps of the simulation the average momentum is far from the prescribed (cos(θ), sin(θ)), and taking into account that the source term is divided by ∆t, (Equation (1)), this produces unphysical strong oscillations of the drag.
To conclude, for the wavy-walled channel simulations, adiabatic wall boundary conditions are enforced at surfaces given by the expressions z w = A w sin(2πx/λ) ( Figure 1). An amplitude of A w = 0.1 has been considered in this work, equivalent to A + w = 18 in wall units at Re τ = 180. The computational setup has the mesh aligned with the λ x /λ y direction (wall crests/troughs along x, see Figure 1b). which makes the computation statistically more homogeneous. The spectral accuracy of the scheme is maintained through mapping from elements Ω e to a reference hexahedronΩ using transfinite interpolation.

Dynamic Mode Decomposition
To uncover the physical feature causing the stabilization effect induced by the wavy wall, we make use of the Dynamic Mode Decomposition (DMD, [15]) technique, which we summarize here.
DMD works with a sequence of instantaneous flow fields numbered from 1 to n s (e.g., taking one or all recorded variables, denoted here generically by v at different time steps of the simulation). The following data matrix is then constructed: where the subindex and superindex identify, respectively, the first-and last-time instances of the sequence. The data are ordered in time, and separated by a constant sampling time interval ∆t s such that: t j+1 = t j + ∆t s for all j = 1, . . . , n s − 1. In the case of linear stability analysis and within the exponential growth region, it is possible to define a linear operator A (i.e., a numerical approximation of the linearized Navier-Stokes operator) such that v(t j+1 ) = A v(t j ). For non-linear systems, A represents the Koopman operator. Equation (3) can be rewritten as a Krylov sequence: For an ordered sequence, Equation (4) can be equated to Equation (3): which can alternatively written in matrix form: Next, the Singular Value Decomposition (SVD) of the snapshot matrix V n s −1 1 = LΣR T is used into Equation (6), leading to A LΣR T = V n s 2 . Further manipulation leads to the reduced matrix: The reduced matrix A, which is the projection of the matrix A onto the space spanned by the left singular vectors matrix L, conveys most of the information codified into operator A. Solving the eigenvalue problem Ay j = µ j y j provides the reduced DMD modes y j and the associated eigenvalues µ j . The approximated eigenmodes of the matrix A can then be recovered via a projection onto the original space, using relation φ j = Ly j . Relation λ j = log(µ j )/∆t s allows recovery of the growth rate (µ j ) ≡ λ j,r and the frequency (µ j ) ≡ λ j,i . Finally, note that the DMD decomposition allows reconstruction of the original data sequence as: In this contribution, the amplitudes α i are computed following the formulation in [22]. That is, the α i 's stem from the minimization problem in the Frobenius norm: where the columns in matrix Φ are the dynamic modes φ i , diagonal matrix D α contains the unknown amplitudes α i and T is a Vandermonde matrix whose columns are generated by the successive powers of the column vector µ k 1 , . . . , µ k n s −1 T , with k = 0, . . . , n s − 1. Since matrix U is unitary, it does not affect the norm in Equation (9), and the optimization problem actually solved is: with the columns in matrix Y the eigenvectors y i of matrix A.

DNS of Planar Walled (Canonical) Turbulent Channel Flow
To validate the current solver and our numerical methodology a Direct Numerical Simulation (DNS) of a turbulent Poiseuille flow between two parallel planes is performed. Similar to other authors, here we address a configuration at a Reynolds number based on the centerline mean velocity and channel half-height Re = U c h/ν = 2820. For this flows the associated Reynolds number based on the friction velocity u τ = τ w ρ corresponds approximately to Re τ = 180 [23,24]. The domain extent and the discretization details are summarized in Table 1, where the + indicate wall viscous scales. Please note that two different meshes are employed: a external (h) mesh consisting of n x × n y × n z elements and an internal spectral (p) mesh of order p x × p y × p z . For Re τ 180 values of ∆ x + , ∆ y + 10 − 15 and ∆ z + wall < 1 are proposed. To achieve this resolution, the multidomain mesh is stretched in the normal direction following a cosine law between −1 and 1. This leads to a wall normal distribution of points where for polynomial order 8, the first and second points are located at z + = 0.07 and 0.7 respectively. Finally, since our Discontinuous Galerkin solver discretizes the compressible Navier-Stokes equations, compressibility effects are minimized by setting a Mach number of 0.1 for all the computations described. In this regard, no distinction is made between velocities and momentums along the paper. For the planar wall case, the simulation is initialized from a laminar Poiseuille profile superposed with a random perturbation with amplitude of 50% of the central velocity. Statistics are obtained averaging over 300,000 time steps taken at ∆ t + = 100. The simulation takes approximately 12 h running in 80 Intel Xeon Gold 6230 cores at 2.1 GHz in the CESVIMA-UPM supercomputer center. Figure 3 shows the average and r.m.s. velocities, obtained by projecting along the edge aligned with the input streamline (θ = 70 • ).

DNS of Wavy Wall Turbulent Channel Flow
The validation of the wavy wall configuration is performed by comparing our results with those obtained by [3], details also in Table 1. The same external mesh with 2000 elements has been reused, but the polynomial degree p x j has been progressively increased as the simulation proceeds. A degree of 8 has been finally set as a good compromise between accuracy and computational costs (see Table 2).
The simulation is initialized from a laminar Poiseuille profile adjusted to the passage section. No perturbation has been necessary to trigger transition: as observed in Figure 2, the flow tends towards a laminar solution with C d 1.120 · 10 −3 , which is slightly above the laminar drag for a planar configuration, C d 1.063 · 10 −3 . This increase is due to the pressure drag introduced by the undulated surface. The transition to turbulence is clearly visible around time 330, where the drag increases suddenly to finally stabilize around C d 4.222 · 10 −3 . The drag reduction from the planar wall channel is about 1.2% (Table 2), which is comparable to the 0.7% reported by [3] for the incompressible case.  Table 2 summarizes the sensitivity of drag reduction figures to polynomial order p x j : the specific values remain consistently around 1%, even when considering the computational setup B, Figure 1c. As a final comparison, an additional computation has been conducted on a domain longer along x, i.e., L x × h × L y = 10.2 × 2 × 10.2, with external mesh density n x × n y × n z = 20 × 20 × 10 and polynomial order p x j = 8 (thus doubling the number of degrees of freedom). The results obtained do not differ significantly from those in Table 2 and are not included here for conciseness. The consistency of the results obtained confirms then that at Ma = 0.1, the wavy wall considered brings a drag reduction which exceeds the pressure drag penalty and agrees with those obtained by [3].

Delay of the Transition to Turbulence
After the validation of our methodology, we address now the main part of this work; the impact that the wavy wall has on transition to turbulence. All the computations reported in this section reuse the 10 × 10 × 20 element meshes with p x j = 8 considered in Section 3: since the Reynolds numbers are lower, a higher accuracy is expected. In order to make a fair comparison between both configurations, planar and wavy, the simulation is initialized from a laminar Poiseuille profile superposed with a random perturbation with amplitude of 50% of the central velocity.
Transition to turbulence in channel flows has been reported to occur for Re ∈ (900, 1000), with some turbulent spots appearing even for Re < 900 [25]. In the simulations conducted in this work, a value between Re 975 − 1000 has been found for the planar case, where the friction drag experiences a large increase from C d 3.08 × 10 −3 − Re = 900 for the laminar regime to C d 4.95 × 10 −3 − Re = 1000 for the turbulent regime. However, for the wavy-walled channel and θ = 70 • , this increase is observed at a higher Reynolds number Re 1175 − 1200, with C d 2.80 × 10 −3 and C d 5.00 × 10 −3 in the laminar and turbulent cases, respectively. Specifically, the transition to turbulent is delayed around a 20% in terms of the Reynolds number.
The effect of the wall undulation is shown in Figure 4, where it is clearly shown the differences between the average vertical velocities. The magnitude of the vertical velocity is an order of magnitude higher for the turbulent case, which is consistent with a higher forcing (drag) term and penetrates much deeper into the channel core compared with the laminar case. Indeed, inspection of the instantaneous vertical velocity, Figure 5, discovers a persistent, organized structure connected to the wavelength of the wavy wall.  The connection of the increased level of organization with the delay of the transition is investigated next using DMD. Two data sequences, laminar at Re = 1150 and turbulent at Re = 1200 respectively, are considered. Each sequence consists of 3D velocity vector fields acquired at a uniform time step of ∆ t s = 0.03 over the non-dimensional time intervals t ∈ (470, 510)-for Re = 1150-and (400, 510) for Re = 1200. For illustration, Figure 6 shows the vertical velocity at a tracer point located at z = −0.935. As expected, the laminar case shows a periodic and well-organized behavior with Strouhal St 0 = 2π/δt 15.9, as opposed to the irregular, chaotic, of the turbulent regime. Please note that time step related to this Strouhal (δt) is twice the time taken by the flow to traverse the wavelength λ x = 5.1, so that the wavy wall introduces passive excitation of a frequency given by the average velocity and the wavelength of the wavy wall. The critical issue here is that this excitation promotes a stable mode able that maintain the laminarity further. The application of the DMD algorithm returns n s − 1 = 1199 modes φ j , eigenvalues λ j and amplitudes α j per data sequence. As usual [16,17], the major difficulty is to discern among all the modes obtained those most relevant for the process investigated. According to Equation (8), the temporal evolution of φ j is controlled by the product α j e λ j t . Here, we introduce a new parameter γ j , defined as: and exploit it to categorize the modes according both to its amplitude α j and growth rate. Figure 7 summarizes the DMD analysis for the Re = 1150 data sequence. The γ j -λ j,i spectrum is shown in Figure 7a: apart from the average value, a dominant mode (highlighted in red) at λ i = 15.93 ≈ St 0 and its harmonics are clearly visible. The streamwise and spanwise components of this dominant mode are visualized in Figure 7b,c. The streamwise perturbation velocity presents a z-uniform, streaky structure, with maximum (minimum) acceleration at the inflection points before (after) the crest. The spanwise component, which is an order of magnitude smaller is clearly modulated and attain more extreme values in the regions close to the wall and the crest of the wave. This spatial structure is better visible in Figure 7d: the positive (crest are, close to the wall) and negative (in the troughs but far from the boundary layer) forces the fluid to move from the troughs to the crest in the spanwise direction. The Re = 1200 flow is more involved, as different structures begin to interact creating complex non-linear patterns. Inspecting the drag evolution (Figure 8a) in t ∈ (400, 510) -just before transition-three different phases can be distinguished. In region ∆T1 (in blue) the solution is still "quasi-laminar" and the mode with St 0 = 15.93 dominates the evolution (and the drag) of the fluid. However, inspection of the DMD spectra (Figure 8b) reveal other incipient modes with a non-negligible relevance (note especially the mode with λ i 3) which may start interacting with other modes. The region ∆T2 (in red), with a sharp rise in drag, shows that the transition to turbulence is taking place. The original mode at λ i = 15.9 remains present, but DMD picks up also new relevant features, λ i = 3, 5.2 and harmonics. As the time progresses, non-linear interactions gain more and more relevance. This effect becomes clearer in region ∆T3, where many new dominant nodes appear as a result of non-linear combinations, showing the path to turbulence.  (1)) before reaching statistical convergence, considering three successive temporal ranges. (b): associated DMD spectra.
As a final analysis, all the computations so far used θ = 70 • [3]. We have been able to show that a fluid moving obliquely to an undulated wall is subject -provided the parameters are well adjusted-to a Stokes layer that delays transition and brings drag reduction. The following question arises then naturally: how important is the sweep angle? To answer this question, three more computations at Re = 1150 and θ = 50, 60, 65 have been conducted. As Figure 9 shows, at θ = 50 the solution is turbulent whereas it is laminar at θ = 60, confirming the fact that a threshold in the skew angle is necessary to create enough spanwise velocity to trigger the stable mode.

Conclusions
In this contribution, we have investigated the effect that wavy walls have on delaying transition to turbulence. Direct Numerical Simulation of both planar and wavy-walled channel flows at laminar and turbulent regimes has been conducted, using our in-house solver based on Discontinuous Galerkin Spectral Element Method.
The validation of the numerical methodology has been performed on the simulation of a turbulent channel at Re τ 180, where we have observed a good agreement between our results and those reported on the bibliography, confirming the moderate drag reduction of about 1%.
In the wavy wall configuration, we have observed a pseudo laminar regime that remains stable until a Reynolds number 20% higher that the critical one of the planar walled channel. Dynamic Mode Decomposition applied to the simulation data revealed that in the wavywalled configuration, modes with wavelength and frequency compatible with the surface undulation pattern are triggered. The relative inclination of the flow over the wavy surface enhances a stable mode that is not shown for the plane wall case. We explained and visualized the appearance of these modes that are also present at higher Reynolds numbers, though they are not dominant anymore.
Certainly, this investigation is still preliminary (e.g.no optimal analysis of the wavelength has been conducted yet). However, the fact that the combination of a wavy surface and a slanted flow, allows the triggering of a stabilizing mode in a completely passive manner encourages to pursue further research into the application of wavy surfaces as passive flow control device.

Data Availability Statement:
The data that support the findings of this study are available from the corresponding author, upon reasonable request.