Identification of flow structures in fully developed canonical and wavy channels by means of modal decomposition techniques

A Dynamic Mode Decomposition (DMD) of Direct Numerical Simulations (DNS) of fully developed channel flows is undertaken in order to study the main differences in flow features between a plane-channel flow and a passively “controlled” flow wherein the mean friction was reduced relative to the baseline by modifying the geometry in order to generate a streamwise-periodic spanwise pressure gradient, as is the case for an oblique wavy wall. The present analysis reports POD and DMD modes for the plane channel, jointly with the application of a sparsity-promoting method, as well as a reconstruction of the Reynolds shear stress with the dynamic modes. Additionally, a dynamic link between the streamwise velocity fluctuations and the friction on the wall is sought by means of a composite approach both in the plane and wavy cases. One of the DMD modes associated with the wavy-wall friction exhibits a meandering motion which was hardly identifiable on the instantaneous friction fluctuations.


Introduction
The future of aviation relies upon more efficient and greener technologies. The reduction of the turbulent skin-friction drag, via passive or active control, can lead to major improvements in efficiency and hence has become a major target of the aeronautical industry. Turbulentdrag-reduction technologies generally aim at reducing the turbulent shear stress in the near wall region. Several DNS-based studies (Bernard et al. [1]) provided direct evidence of a link between the production of Reynolds stresses and the quasi-streamwise vortical structures in the buffer layer. It was shown in Kravchenko et al. [2] that streamwise vortices tend to be located above and to the side of high-skin-friction regions, which suggests an intimate link with the wall shear stress and consequently skin friction.
In particular, imparting transverse in-plane wall motions was shown to yield considerable levels of drag reduction (up to about 45% in Quadrio et al. [3]). In an effort to emulate the effect of wall motions using a passive device, Ghebali et al. [4] studied the effect of large-scale oblique undulations on a solid (non-moving) wall. The rationale is that although not exactly identical to the in-plane motions, the wall undulations generate a spanwise shear layer that is quantitatively similar to that arising from the imposition of an in-plane standing wave. The results showed a maximum overall-drag reduction of 0.6% at Re τ ≈ 360. In this study, flow structures are extracted for a close configuration (with the parameters scaled in wall units) and compared to those found in a plane channel, at Re τ ≈ 180. In order to extract the flow structures, two different modal-analysis techniques are considered, namely Proper Orthogonal Decomposition (POD) and Dynamic Mode Decomposition (DMD) (cf. Taira et al. [5] for a general overview).
The POD technique benefits from a long history and numerous studies of laminar and turbulent flows can be found, e.g. in [6,7]. However, the developments in DMD techniques [8] are more recent; applications include laminar [9,10], transitional [11] and turbulent flows [12,13] with clear tonal frequencies, but the application of DMD to fully developed turbulent flows (with a broadband character) is still limited and is the focus of this paper. In this contribution, both POD and DMD are applied to assist in proposing links between the flow structures and friction drag.

POD and DMD
Dynamic Mode Decomposition (DMD) techniques [8,14,15] have increased in popularity over the past years and have been applied to a variety of flows, using data from numerical simulations and experiments (e.g. [10,15,16,12]).
First, a brief summary of the snapshot-based DMD technique, as proposed by Schmid [8], is presented. Given a sequence of instantaneous flow fields numbered from 1 to N (e.g. taking one or all recorded variables), the following data matrix can be constructed: where the subindex and superindex identify, respectively, the first and last time instants of the sequence. The data is 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 − 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 linearised Navier-Stokes operator) such that v(t j+1 ) = Av(t j ). Equation (1) can then be rewritten as a Krylov sequence (Saad [17]): For an ordered sequence, equation (2) can be equated to equation (1): which can alternatively be written in matrix form as: Next, the Singular Value Decomposition (SVD) of the matrix V 1 N −1 = UΣW H is obtained, where the superscript H denotes the conjugate transpose. The left singular vectors -the columns of U -correspond to the POD modes of the input data sequence. Thus, the snapshot-based DMD algorithm requires the calculation of the POD modes. The SVD of the snapshot matrix is then inserted into equation (4), which yields AUΣW H = V 2 N . The reduced matrix S = U H AU associated with the initial system described by A, can be rewritten using the former equality as: The reduced matrix S is the projection of the matrix A onto the POD space contained in U, and previously obtained through a SVD [8].
Once the reduced matrix S has been calculated, the reduced DMD modes y i can be obtained, as well as the associated eigenvalues µ i (i.e. growth rates (µ i ) and frequencies (µ i ) mapped to the unit circle) of the reduced system by solving the eigenvalue problem Sy i = µy i . The approximated eigenmodes of the matrix A can then be recovered via a projection onto the original space, using φ i = Uy i . Eventually, the growth rates and frequencies in the complex half-plane can be recovered from the eigenvalues as: λ i = log(µ i )/∆t s .
Various extensions of the DMD algorithm have been proposed after the original method described by Schmid [8]. In particular, the sparsity-promoting version of Jovanović et al. [19] is also leveraged in this study, in an attempt to reduce the number of representative modes.

Numerical simulations
The instantaneous realisations analysed in this study were generated by DNS using an in-house code solving the incompressible Navier-Stokes equations. The finite-volume discretisation of the momentum fluxes is second-order accurate, and the solution fields are integrated in time by means of a fractional-step method with a third-order gear-like scheme. The projection step is effected by solving a pressure-Poisson equation, and stable velocity-pressure coupling is ensured by the use of Rhie-and-Chow interpolation. Non-planar geometries are simulated with a body-fitted mesh. The flow is driven by a spatially homogeneous pressure gradient in the streamwise and spanwise directions, adjusted at each time step in order to enforce a unit streamwise bulk velocity and zero spanwise bulk velocity. The fluid viscosity is set so that the bulk Reynolds number is Re b = U b h/ν = 2820, corresponding to a friction Reynolds number Re τ = u τ h/ν ≈ 180. Each simulation is initiated with a laminar flow, onto which random perturbations are added while the viscosity is temporarily reduced. The time advancement is continued until a fully developed turbulent state is reached, corresponding to the time instant t = 0.

Simulation set-up
The simulation parameters are shown in table 1 for the plane and wavy channels.
The computational domain used for the wavy channel is shown in figure 1 and comprises one single undulation at an angle of 20 degrees to the main flow direction (corresponding to θ = 70 • in [4]), of amplitude A w = 0.15h and wavelength λ + x = 15h + ≈ 2700. The distance between the upper and lower walls is constant, and the walls are located at The overall drag in this case is very close to the baseline as a consequence of the pressure drag balancing or even negating the reduction in friction, but benefits from a large amplitude which was chosen to increase the effect on the flow structure in order for the differences to be larger. As a consequence, this flow configuration may shed a more contrasted light than a small wave slope would, on two different phenomena: 1) the reduction in friction, and 2) the effects associated with the pressure drag responsible for a significant drag penalty. In view of the very fine resolution required to accurately measure the drag level, as shown in [4], this study does not aim at quantifying the precise level of drag, but rather seeks to understand the underlying physics by identifying stronger changes in the coherent structures.

Databases
Extended databases consisting of instantaneous snapshots of the velocity and pressure field were generated for the two simulations. As shown in figure 1 for the wavy channel, the raw database consists of all the cells contained in the lower half of the solution domain, saved every 100∆t. In order to cut the memory storage required during the post-processing, the 3D volumes were filtered by discarding data points and retaining information at regular space and time intervals ∆ s i , shown in table 2.
The consistency between the friction coefficient calculated from the driving pressure gradient and that calculated from the reduced data base is shown in figure 2 for a reduced volume of the snapshot database. The reduction in the number of statistical samples arising from the partial use of the data contained in the numerical domain results in an increased standard deviation for the friction coefficient in time with respect to that of the driving pressure gradient which accounts for the entire volume simulated. The total stress on the wall can be calculated as: where S w is the surface of the wall, e x is the streamwise direction unit vector, n and t are, respectively, the local unit normal and tangential vectors, and u t is the tangential velocity component. However, only the frictional part (without the pressure contribution) is considered in this study, and its calculation contains a certain degree of approximation, as the wall-normal velocity gradient is approximated by ∂u/∂y. Although a more careful calculation would be required for a precise measure of the friction, it is reasonable to expect that the approximations made entail an error level lower than or comparable to the spatial-discretisation error, and make little difference in the main flow dynamics.

POD analysis
Since the DMD analysis involves a projection onto the POD modes, the latter were studied in detail for the plane-channel case. The POD modes were calculated using the method of snapshots with all three components of the velocity field. The influence of the proportion of the database used in the POD analysis was mostly evaluated by trial and error. The spatial resolution was selected by considering a small subset of snapshots and varying the resolution in each direction separately. When every other snapshot is included in the data sequence, the POD modes are not sensitive to a shift in the data sequence by one snapshot, which suggests that the frequency of accumulation of the snapshots is sufficient. As a consequence, the snapshot sampling frequency in the wavy case was chosen to be half that of the plane case.
In all the results presented in this section, the data sequence comprises snapshots every 50∆t, and 4∆x × ∆y × 2∆z, in a volume of size L x × h × L z /4. Although increasing the number of snapshots seemed to preferentially yield particular mode shapes, not all of them were found to be robust against a spanwise translation of the volume of snapshots considered. This may be attributed to the limited streamwise extent of the solution domain, in which structures are recirculated via the periodic boundary conditions.
As observed by Nikolaidis et al. [7], the convergence of the POD modes for the turbulent plane channel is very slow. To circumvent this limitation, Nikolaidis et al. [7] took advantage of   5 show isosurfaces of the POD modes for each component of the velocity vector, for a sequence of 5817 snapshots spanning a duration of over 10, 000 u 2 τ /ν, which represents about 57 eddy-turnover times. Although no statistical convergence can be claimed even for this relatively long snapshot sequence, some broad similarities are found in the spacing of the rolls, relative to the findings of [7], where the spacing of some of the first streamwise POD structures was shown to be of about 70 + , 100 + , and 150 + , which is in line with, respectively, modes {5, 6}, {1, 3}, and 2. Moreover, most of the modes show a correlation between the streamwise and the wall-normal components of the velocity, with ejections correlating with low-speed streaks, and sweeps with high-speed regions, consistently with the findings of [7] who interpreted this correlation as arising from the action of the lift-up mechanism maintaining the streaks.

DMD of the velocity field (plane channel)
A DMD of the velocity field was undertaken in an attempt to identify a limited set of lowfrequency modes able to reproduce the mean Reynolds shear stress from the dynamical modes. Such a procedure was studied by Sayadi et al. [11] for a transitional boundary layer. For brevity, the results shown are limited to the plane-channel case since the procedure is not fundamentally different in the wavy case. First, a DMD of the velocity-fluctuation vector is performed, second a sparsity-promoting technique is applied, and third, the mean Reynolds shear stress extracted  Note that in all cases, although the raw snapshots are saved in single precision, performing the calculation of the DMD modes in double precision was found to be an absolute necessity to ensure that the eigenvalues were close to the unit circle.

DMD modes
The frequencies and corresponding amplitudes are shown in figure 6. Some peaks were identified, and their position in the spectrum was found to be robust against a change in the length of the data sequence considered (by comparing the results with a reduced number of snapshots of 101 and 301; this observation does not apply to the lowest frequency peak which is not resolved when 301 or fewer snapshots are used). One peak appears to be located at T + = O(100), which may be related near the wall to a spatial dimension of λ + = O(1000) using Taylor's hypothesis and the convection velocity at the wall u + ≈ 10 (see e.g. delÁlamo et al. [18]). This suggests that the modes surrounding this peak might be linked to the near-wall streaks, and it is remarkably close to the optimal actuation period for oscillating walls T + = 125 in Quadrio et al. [3]. Additionally, the large amplitudes of the low-frequency modes are brought to light, and may suggest the presence of long structures or reflect the occurrence of structure recirculation arising from the finite domain length.
The DMD modes corresponding to the peaks identified in figure 6 are shown in figure 7. Although the mode numbers do not bear any physical meaning, unlike with the POD where they are ranked by energy content, the mode number is used as a more convenient label to identify  Interestingly, the mode number 4583, corresponding to T + ≈ 100, displays elongated structures, with low-speed and high-speed regions aligned at an angle to the flow. The u-structures are inclined, starting from the wall and reaching the outer region. Again, the low-speed u-structures correlate with positive wall-normal velocities (ejections) whilst highspeed ones correlate with negative values (sweeps). The spanwise structures are observed to be much flatter and wider, and are tilted in both the x-and y-directions.

Sparsity promotion
The sparsity-promoting technique by Jovanović et al. [19] was applied to the DMD for the plane channel using a set of 301 snapshots with the spatial and temporal spacings shown in table 2.
The Matlab source code available from http://people.ece.umn.edu/~mihailo/software/ dmdsp/ was employed. As the penalty γ is increased, the number of modes with non-zero amplitude decreases, and so does the fidelity of the modal reconstruction, as quantified in figure 8(a) using the same norm as Jovanović et al. [19]. Figures 8(b) and 9 show the spectrum corresponding to the highest value of γ (smallest number of non-zero-amplitude modes) compared to the full spectrum. The broadband nature of the fully developed turbulent flow makes the sparsity-promoting algorithm fail to isolate few modes, and the performance quickly decays with the level of sparsity.

Reconstruction of the Reynolds shear stress
Using the large set of 5818 snapshots, the DMD modes of the streamwise and wall-normal velocities are used in order to estimate the corresponding mean Reynolds shear stress. Figure 10 shows a comparison of the Reynolds shear stress obtained from DNS and that obtained from the sum over a range of modes of the contributions from the product of the DMD modes of the streamwise and wall-normal velocities (cf. Sayadi et al. [11]). The good statistical convergence of the Reynolds stress from DNS is demonstrated by the fact that the dash-dotted line -corresponding to the DNS snapshot sequence -is smooth and does not feature any highfrequency oscillation away from the wall. Unlike in the transitional flow of [11] where only a few low-frequency modes sufficed to estimate the shear stress, it was found that a significant proportion of the modes is necessary in the fully-turbulent plane channel. This is a known limitation of the DMD algorithm associated with the absence of a clearly defined frequency in the dynamics, so that important structures are spread across a large range of frequencies, both temporally and spatially. This also echoes the difficulty faced to isolate a small number of modes with the sparsity-promoting technique in section 4.2.

Composite DMD
A composite DMD was performed in order to seek a dynamical link between the friction at the wall and other quantities such as velocity, pressure, or the Reynolds shear stress by appending the wall friction in the snapshot matrix. Thus, the snapshot matrix in the results shown in this section contained the streamwise fluctuations scaled with the bulk velocity, together with the fluctuations of the non-dimensional friction coefficient. Owing to resource limitations, only a small sequence of 301 snapshots is considered in the present analysis. However, as a consequence of the convergence issues evidenced in section 3, the modes obtained may be only relevant to the dynamics of the snapshot sequence selected, and not necessarily to that of the turbulent flow associated with the flow configuration. Figures 11 and 13 show the DMD spectra obtained for the plane and the wavy case using a sequence with similar spacing (cf. table 2) and composed of 301 snapshots. The amplitude of the lowest-frequency mode is found to be higher in the wavy case, and the frequencies are more evenly distributed than in the plane case. Some low-frequency modes are shown in figures 12 and 14 for the plane channel and the wavy channel, respectively. In general, low-speed regions     Figure 11. DMD spectrum for the plane channel, using the first 301 snapshots. Wall units are based on the actual friction velocity. correlate with low-friction regions, and conversely high-speed regions correlate with high-friction levels. Again, there seems to be a preference in the alignment of the low-and high-speed regions, making an angle of about 10 degrees with the main flow direction. Moreover, the threedimensional arrangement of the streamwise-fluctuation modes presents tilted structures (around the z-axis) with low-and high-speed regions overlapping each other. The DMD mode number 3 is shown for the wall friction in figure 15(b), and compared to an instantaneous view of the wall friction from DNS in figure 15(a). Although the instantaneous contours mostly show elongated structures aligned with the main flow direction, the DMD mode brings to light a meandering motion resulting from the presence of the oblique wavy boundary.  Figure 13. DMD spectrum for the wavy wall, using the first 301 snapshots. Wall units are based on the baseline friction.

Summary and outlook
Simulations of fully developed turbulent channel flows were performed, at Re τ ≈ 180, with plane and wavy walls. Instantaneous realisations of these flows were analysed by means of POD and DMD. The DMD revealed that most of the modes are necessary to represent the flow accurately, as a consequence of the broadband spectrum of the flow dynamics. Eventually, a composite DMD was undertaken in order to find dynamic correlations between the wall friction and the three-dimensional streamwise-fluctuation field. It was observed for the wavy wall that one of the wall-friction modes could enable an alternative representation of the effect of the wavy wall, highlighting a spanwise-wavy pattern that was not identifiable from the instantaneous wallfriction. Future work will include a deeper assessment of the selection of the snapshot sequence, gathering more data by leveraging parallel map-reduce algorithms. This would enable higher Reynolds numbers to be processed, and potentially shed light on the inner-outer interactions and their consequences on the turbulent skin friction. The composite DMD may be studied in more detail for the pressure fluctuations, or kinetic-energy production or dissipation, especially in the case of the wavy wall. Eventually, it may be interesting to consider very small regions in a Lagrangian framework [20] in order to reduce the dilution of the travelling structures across a broad range of modes.