# Exact coherent states for grooved Couette flows

###### Abstract

Recent progress indicates that highly symmetric recurring solutions of the Navier-Stokes equations such as equilibria and periodic orbits provide a skeleton for turbulence dynamics in state-space. Many of these solutions have been found for flat-walled plane Couette, channel, and pipe flows. Rough-walled flows are of great practical significance, yet no recurring solutions are known for these flows. We present a numerical homotopy method to continue solutions from flat-walled plane Couette flow (PCF) to grooved PCF, demonstrated here at a Reynolds number of 400, to act as a starting point for similar continuation to rough-walled flows. Loss of spanwise homogeneity in grooved PCF reduces continuous families of solutions (identical up to translational shifts) in flat-walled Couette flow to multiple, discrete families in grooved Couette flow; this can manifest in turbulence as spatially anchored exact coherent structures near the wall, so that turbulent statistics reflect symmetry-restricted structure of exact recurring solutions. Furthermore, the vortex-streak structures characteristic of these equilibria are squeezed out of the grooves when the groove-wavelength is smaller than the characteristic spanwise size of the structures. This produces reduced shear stress at the wall even at the low Reynolds numbers considered, and the mechanism is consistent with the drag reduction observed in some riblet-mounted turbulent flows.

## 1 Introduction

For most of the twentieth century, turbulence was studied through the lens of statistics despite knowledge of the existence of recurring structures. Townsend’s attached eddy hypothesis (Townsend 1980) was an early attempt to describe turbulence in terms of coherent structures and relate their scaling to observed statistics. These eddies were suspected to be the building blocks of turbulence (Perry & Marušic 1995), but they were ad-hoc structures without basis in the Navier-Stokes equations (NSE). Advances in computational power allowed computation of so-called exact coherent structures (ECS), starting with plane Couette flow (PCF) equilibria of Nagata (1990). Minimal channel DNS of Jiménez & Moin (1991) confirmed the significance of the interaction between streamwise streaks and vortices seen in the equilibria. Waleffe (1997) provided a detailed picture of the recurring motions that occur in wall-bounded turbulence; the linear instability and non-linear breakdown of streaks produces vortices, which in turn feed into the streaks.

In the past decade, many exact solutions of the Navier Stokes equations have been found and their connections in state-space have been mapped (for instance Gibson & Cvitanović 2008; Eckhardt et al. 2008; Kawahara et al. 2012; Park & Graham 2015; Willis et al. 2017). A plausible description of the underlying dynamics at asymptotically high Reynolds numbers has also been proposed (Hall & Smith 1991; Hall & Sherwin 2010) highlighting the inviscid mechanism of vortex-wave interactions. Today, these solutions — equilibria, travelling waves, periodic orbits, and relative periodic orbits — are understood to shepherd turbulent trajectories in state-space.

While rapid progress has been made towards understanding turbulence dynamics in terms of exact recurrent solutions of the NSE, such efforts are yet to involve wall-roughness, despite the practical significance of rough-walled turbulence. Studying the variation in the structure and statistics of exact solutions with changing wall geometry (from smooth to rough) can offer insight into the interactions of wall-roughness with turbulent flow. Such an approach can complement the work to date on rough-walled turbulence (Jiménez 2004; Flack & Schultz 2014; Squire et al. 2016), where turbulence structure and statistics are examined using experimental and numerical methods, to provide a quantitative basis for characterizing surface topologies and relating them to measures such as equivalent sand roughness.

We consider the special case of riblets (longitudinal grooves) to serve as a starting point for the extension of ECS from smooth-walled to rough-walled geometries; we continue the equilibria of Nagata (1990), which were computed independently and hosted on the website www.channelflow.org by Gibson (2014). In addition to the geometrical simplicity, riblet-mounted flows are particularly interesting because of the drag reduction observed for certain arrangements (Walsh 1990; Dean & Bhushan 2010; Bannier et al. 2016). This drag reduction is accompanied by a localization of the vortex-streak structure near the wall (Lee & Lee 2001). We expose the drag reduction and localization to be related to the influence of riblets on symmetric ECS.

The original solutions of Nagata (1990) were obtained using homotopy from Taylor-Vortex flow to PCF. We use homotopy to extend the equilibria of Nagata (1990) from flat-walled PCF (smooth-walled PCF) to riblet-mounted PCF using a domain transformation method. Henceforth, we refer to flat-walled PCF as flat PCF, and riblet-mounted PCF as grooved PCF. We also distinguish the base flow (and solutions obtained from continuing this to grooved PCF) as the laminar solution, while the 3-dimensional solutions of Nagata (1990) and their continuations are referred to as equilibria.

### 1.1 Flow symmetries

#### 1.1.1 Continuous symmetries

Smooth-walled PCF involves infinite, parallel plates moving in opposite directions. This allows for solutions that are invariant under certain symmetry transformations; we follow the notation of Gibson et al. (2009) to represent these symmetries. We use non-dimensional coordinates along the streamwise, wall-normal, and spanwise directions respectively, non-dimensionalized by half of the distance between the walls.

Homogeneity in the streamwise and spanwise directions leads to admission of solutions that are invariant under continuous translation in the plane; the symmetry transformations for these translations become the continuous two-parameter group , denoted by

(1.0) |

This homogeneity also means that spatially periodic solutions can be sought, so that computations can be carried out over small domain sizes using the Fourier spectral method for discretization. Different box sizes produce different solutions; however, provided that the box-sizes are within certain ranges — of the order of the channel half-height — these solutions all seem to contain the vortex-streak structure.

While the homogeneity allows for relatively inexpensive, spatially periodic solutions, it also produces vast families of solutions that are dynamically similar. That is, every equilibrium is part of a continuous family of solutions spanned by for arbitrary , and can also have dynamically similar travelling waves associated with the family. Presence of such families could obscure our view of turbulent flows. Willis et al. (2013) developed a projection method called the “method of slices” to identify families of dynamically similar solutions with one representative solution for each family, thus reducing the number of solutions that would have to be considered, say, when developing low-order models for control. Specifically, flat-walled plane shear flows admit travelling wave solutions and relative periodic orbits, and the method of slices relates such solutions to dynamically similar equilibria and periodic orbits by projecting the former onto a hyperplane using a pre-defined template. Grooved flows, however, are inhomogeneous in the direction transverse to the grooves. The consequences of this loss of homogeneity are explored in section 3.1.

#### 1.1.2 Discretization

Owing to the existence of spatially periodic solutions for channel flows and PCF, the flowfield can be discretized using the Fourier spectral method in the wall-parallel directions. A periodic field for a periodic box of streamwise and spanwise lengths and is written as a sum of Fourier modes,

(1.0) |

where . The Fourier modes are truncated to some and to obtain a finite-dimensional representation of the field variable in and . Henceforth, hats on the coefficients of Fourier modes are dropped for convenience, with the understanding that integer-subscripted symbols represent Fourier coefficients.

The Fourier coefficients are functions of the wall-normal coordinate and time. Chebyshev collocation method is used for discretization along the wall-normal direction. The time dependence is usually accounted for by marching forward from an initial condition; but in the steady-state solver used in this work, this dependence is dropped.

#### 1.1.3 Discrete symmetries

Smooth PCF also admits discrete symmetries, which are restricted by the counter-moving walls to reflection across planes (), rotation about the -axis (), and point-wise inversion about the origin (), which is the product of the reflection and rotation. These symmetry transformations generate a discrete dihedral group , where

(1.0) | ||||

and is the identity transformation. The equations of flat PCF are invariant under the group , where stands for a semi-direct product, subscripts indicate translations along and sign changes in and , and subscripts indicate translations along and sign changes in .

The laminar solution (base flow, ) for the flat PCF system is invariant under every symmetry transformation in . Equilibria, such as those of Nagata (1990), are invariant under fewer symmetries. Streamwise travelling waves cannot exist when is imposed, spanwise travelling waves cannot exist when is imposed, and neither can exist when is imposed. More complex solutions appear as these symmetries are relaxed, with the fully turbulent state unlikely to satisfy any of the symmetries in (although discrete translational symmetries are imposed in simulations in the form of periodic boxes). The reader is referred to Gibson et al. (2009); Halcrow (2008) for a comprehensive treatment of invariance under discrete symmetry groups. They identify one particular order-4 group, , to be important to finding equilibria of flat PCF, where

(1.0) | ||||

are products of half-cell shifts and the discrete symmetries of . The solutions used for continuation in the present work are invariant under the symmetries of .

Due to the high-dimensionality of the state-space of flat PCF, and fluid flows in general, a large number of exact invariant solutions exist; in fact, an infinite number of solutions exist due to the continuous translational invariance of the governing equations. A small set of dynamically important solutions has to be identified in order to realistically describe and predict turbulence. The symmetry groups , , and other specific subgroups of offer one way of reducing the infinite stationary points of the NSE to a small set, while the method of slices described earlier offers a complementary way to constraining continuous symmetries.

## 2 Methodology

### 2.1 Domain transformation for grooved PCF

A general case of wall roughness periodic in can be expressed as a Fourier series,

(2.0) |

This grooved geometry can be mapped to a flat-walled geometry using a domain transformation, (for similar usage, see Kasliwal et al. (2012); Moradi & Floryan (2013)), so that the Fourier spectral method is easily applied in the transformed domain. Mapping partial derivatives from the transformed domain to the physical domain of section 2.1 is, in general, rather cumbersome. However, when and , i.e., when the width between the top and bottom walls remains constant, the required domain transformation is

(2.0) |

which produces simpler relations between the partial derivatives in the physical domain of section 2.1 and the transformed domain.

We investigate the special case of longitudinal grooves (), with the additional simplification that . The top and bottom walls of the PCF are given by

(2.0) |

and the grooved PCF is mapped to a flat PCF geometry as

(2.0) |

The partial derivatives in the two domains relate as

(2.0) |

### 2.2 Iterative solver

A primitive variable formulation is used in a steady-state solver to extend the equilibria of Nagata (1990) to grooved PCF; the solver is outlined in appendix A. The velocities along are non-dimensionalized by the speed of the walls; the reference frame is chosen such that the walls have velocities at . The pressure is non-dimensionalized by , where is the density of the fluid. The equations to be solved are

(2.0) | ||||

where . The boundary conditions are

(2.0) |

The Reynolds number is based on the dimensional half-width and the dimensional speed of the walls.

The velocities and pressure are mapped from the grooved PCF domain with walls given in section 2.1 to a flat PCF domain according to section 2.1; the boundary conditions in the transformed domain are

(2.0) | |||||

The fields are discretized using the Fourier spectral method along the streamwise and spanwise directions and the Chebyshev collocation method along the wall-normal direction in the transformed domain . Section 2.2 is discretized by writing the partial derivatives in in terms of partial derivatives in using section 2.1. A major consequence of the domain transformation is that the Fourier coefficients of the spanwise derivatives involve multiple Fourier modes:

(2.0) |

where is the -dependent Fourier coefficient for mode . This additional inter-modal interaction contains all the information due to the grooved walls, and is the trade-off for the simplified boundary conditions of section 2.2.

The discretized equations are solved using the Newton-Raphson method (see appendix A). A full-rank matrix inversion (SVD-based) is used in the iterative correction stpdf. The iterations are initialized with exact solutions for flat PCF: EQ1 (lower branch) and EQ2 (upper branch) of Nagata (1990), taken from the solutions database of Gibson (2014). Our steady-state solver defines instantaneous time-derivative and divergence of the velocity field as the residual (), and the residual norm is defined as the integral

(2.0) |

Solutions for grooved PCF are computed with 35 wall-normal Chebyshev nodes and 825 Fourier modes: and along the streamwise and spanwise directions respectively. Spatial accuracy of the solutions is defined as the residual norm of the solutions when extrapolated on a grid with 70 Chebyshev nodes and 3185 Fourier modes: and .

## 3 Results and Discussion

The wavelength of the grooves, defined by , is set to coincide with the spanwise periodic length of the original flat PCF equilibria: . The streamwise length of the periodic box is also set according to the equilibria to be continued: . All solutions are for . The flat PCF solutions EQ1 (lower branch) and EQ2 (upper branch) are shown in figs. 1(b) and 1(a); contours of streamwise velocity and quiver arrows for cross-stream velocities are plotted. The continued solutions for grooved PCF are shown in figs. 1(d) and 1(c) at . The vortex-streak interaction in grooved PCF can be seen from the isosurfaces of zero streamwise velocity and streamwise vorticity shown in fig. 3. The residual norm for these solutions is lower than and their spatial accuracy is .

Solutions for grooved PCF look similar to their flat PCF counterparts, except for some distortion in the streamwise vortices. Streamwise streaks are accompanied by streamwise vortices centered on the low speed streaks. This is not surprising; the interaction of streamwise streaks and vortices is known to be a hallmark of wall-turbulence and has been shown to exist irrespective of the kind of wall-bounded flow (Waleffe 1997; Blackburn et al. 2013).

### 3.1 Phase restriction

The flowfield shown in fig. 1(a) is just one member of a continuous family of solutions that arise due to spanwise homogeneity of flat PCF. This solution satisfies the shift-reflect symmetry about the plane . We refer to this solution as the zero-phase solution and denote it . Let us denote a flowfield that is obtained by translating along by some shift as . This flowfield still satisfies a shift-reflect symmetry, but about the plane instead of ; we say that these solutions have non-zero phases, with the phase being . The Fourier coefficients, denoted for mode , of can be related to the Fourier coefficients of as

(3.0) |

so that spanwise translation of flowfields can be seen as phase-shifts in the Fourier coefficients.

Let us also define the continuous family of flowfields obtained by translating along as , and the discrete family of flowfields obtained by translating along through integral multiples of as . The argument for is dropped for convenience. Every is a subset of . The members of are identical to each other due to the periodicity of .

For flat PCF, because of its spanwise homogeneity, every member of is an equilibrium of the NSE. When is used as the initial iterate for finding a grooved PCF equilibrium with groove amplitude , it converges to (shown in fig. 3 for ). However, when a different member of is used as the initial iterate, convergence is not likely. We find that the iterations converge only for flowfields , or for flowfields sufficiently close to members of ; the converged solutions for the latter case are identical to the solutions obtained from continuing the corresponding member in . Failure of iterations to converge indicates that the initial iterate is not in the vicinity of an equilibrium. We take this as indicating that most members of do not have a corresponding equilibrium in grooved PCF.

The significance of can be attributed to a discrete symmetry of PCF. The original equilibria, EQ1 and EQ2, are invariant under the discrete symmetries , , and . The symmetry is a composition of reflection along the plane , and a half-cell shift along ; the reflection and -shift commute. Because of the spanwise homogeneity of flat PCF, it admits solutions with reflectional symmetry about any plane, such as a flowfield which is shift-reflect invariant about the plane . By contrast, grooved PCF has reflectional symmetry only about particular planes, as can be seen in fig. 1. In , reflectional symmetry holds only for planes and . So, any equilibrium of grooved PCF can have reflectional symmetry only about these two planes. Initial iterates from conform to this symmetry even if they do not exactly satisfy the NSE, while flowfields in do not even satisfy the symmetry. Our results show that the solutions obtained by continuing also satisfy the shift-reflect symmetry .

#### 3.1.1 Bifurcation

has two non-identical members: and . Let us denote their continued solutions at some amplitude by and respectively. Our results show that ; the members of bifurcate into two distinct families of solutions in grooved PCF: and , with distinct values for scalars such as energy density and bulk dissipation rate. The difference arises because, while grooved PCF has reflectional symmetry about both and , the walls at are displaced along and the walls at are displaced along . This bifurcation can also be interpreted in terms of continuing along positive and negative groove amplitudes.

#### 3.1.2 Spatial anchoring

Willis et al. (2013) argue that the presence of continuous families of solutions could obscure our description of turbulent flow and describe a method of slices to systematically reduce a family of similar solutions to one representative solution; with this approach, travelling wave solutions reduce to equilibria, and relative periodic orbits to periodic orbits. Grooved PCF achieves this reduction by physically disrupting some of these symmetries. Spatially localized equilibria have also been found in flat PCF (for instance Gibson & Brand 2014; Chantry & Kerswell 2015). Based on our results for global solutions, we can speculate that any symmetric localized solutions, when extended to grooved PCF, would also be phase-restricted by the grooves.

In flat PCF turbulence, any of the solutions in a family could be visited; no member of a family is special. As a result, turbulence statistics would only show spatial averages over an entire family of solutions. The phase-restriction of solutions in grooved PCF, however, implies that turbulence statistics should reflect the spatial structure of the solutions.

### 3.2 Wavelength dependence

We now introduce multiple grooves within the periodic box of size , with the walls given for these cases by

(3.0) |

The term is included to retain the fundamental periodicity along .

Figure 4 illustrates continuation of EQ1 to three different geometries — one, two, and three grooves per box — at a dominant groove amplitude of . EQ2 is plotted for the same cases in fig. 5. The two groove per box geometry does not satisfy the shift-rotate symmetry or the shift-invert symmetry ; these cases were run at a resolution of 20x35x32 grid points () to converge to residual norms and accuracies for EQ1 and for EQ2. The three groove per box satisfies , , and , and these cases were run with 28x35x32 grid points () like the one groove per box case, and have residual norms similar to the one groove per box case. The vortex-streak structure remains the same for the two and three grooves per box cases as well. However, the large scale distortion seen in the one groove per box case is not as pronounced in the other two cases.

#### 3.2.1 Mean velocity

The mean velocity profiles for all three wavenumber cases are shown in fig. 6, with the profiles for grooved PCF shown as a difference from the mean velocity for flat PCF. In the profiles of the one groove per box cases, the deficit in the mean velocity (compared to ) is increased by increasing groove-amplitude. In contrast, the profiles of the two and three groove per box cases bridge the velocity deficit in the mean velocity of both EQ1 and EQ2, and the mean velocity grows closer to the profile as the amplitude increases. The difference in mean velocity for EQ2 for the two groove per box is not zero at the centerline; this is because of the loss of the shift-rotate symmetry (see section 1.1.3) for the two groove per box case.

The bridging of the velocity deficit is stronger for EQ2 than EQ1, although the velocity deficit for EQ2 is much higher than that for EQ1 in the first place. This effect is also stronger for the three grooves per box case than it is for the two grooves per box case. This decrease in velocity deficit also leads to the shear stress at the wall for grooved PCF becoming closer to that for the laminar flow in flat PCF. The one groove per box case increases the shear stress at the wall compared to the flat wall case, while the two and three groves per box cases reduce the shear stress at the wall. While this hints at a drag reducing tendency for the two and three grooves per box cases, the profiles plotted in fig. 6 extend only until the tips of the grooves. They do not show the actual stresses at the wall. An exact quantification of the drag is provided next through plots of dissipation rate.

#### 3.2.2 Energy density and dissipation rate

The kinetic energy density , bulk dissipation rate , and power input of a velocity field of plane Couette flow are given as

(3.0) | ||||

where is the volume of the periodic domain, is an infinitesimal arc-length along the wall in the spanwise direction, and is the wetted surface area of each wall. These quantities are normalized so that for laminar flow in flat PCF, and . For equilibria, , , and are not time-dependent, and .

Following the equality of power input and bulk dissipation rate, the drag coefficient () can be related to the bulk dissipation rate () as

(3.0) |

For the constant- continuation discussed in this chapter, the drag coefficient differs from the bulk dissipation rate by a constant factor of ; so, it is sufficient to plot one of either the drag coefficient or bulk dissipation rate. In keeping with the literature, bulk dissipation rates are shown.

The kinetic energy density and bulk dissipation rate for flat PCF for the laminar solution and upper and lower branch equilibria are shown in table 1.

Turb. Mean | Laminar | EQ1 (lower) | EQ2 (upper) | |

0.087 | 0.1667 | 0.1363 | 0.0780 | |

2.926 | 1 | 1.429 | 3.044 |

The same quantities for the three wavenumber cases of grooved PCF are shown in fig. 7 for the laminar solution and the equilibria, EQ1 and EQ2. The plots for the one groove per box case are reproduced to facilitate comparison.

The curves for the laminar solution are discussed first to serve as a base case. All of the curves for the laminar case have zero derivative at because of the highly symmetric nature of the laminar solution. The kinetic energy density and the dissipation rate increase with increasing groove-amplitude, and this increase is stronger for the lower wavelength cases. The increase in kinetic energy is easily understood in terms of the increase in wetted surface area (per unit planform area). The sharper increase in kinetic energy for the three grooves case also follows from the increase in wetted surface area. Bulk dissipation rate, in the laminar case, is mostly due to the wall-normal derivative of the streamwise velocity. Increasing groove amplitude brings high speed regions (of opposite directions) closer, which in turn produces stronger shear and larger dissipation rates.

We now turn our attention to the kinetic energy density for EQ1 and EQ2 for the three geometries. For the two and three groove cases, the energy density increases, and we can expect this to be for the same reason as the laminar case. Indeed, in figs. 5 and 4, the sizes of the high speed regions (blue at the bottom wall and red at the top wall) do increase as they fill up the area within the grooves.

The anamolous case is the one groove case, where the energy density decreases with increasing amplitude. This can be explained in terms of the vortex-streak structure “eating into” the high speed regions near the wall. The vortex-streak structure is based around the critical layer (it is a critical layer because the equilibria travel with a phase speed of zero). The critical layer moving closer to the wall is accompanied by a shrinkage of the high speed regions, which leads to reduced energy density. It is worth noting that the decrease in energy density for flat PCF from laminar to EQ1 to EQ2 is also linked to the critical layer growing closer to the wall.

This mechanism for reduction in energy density is most clearly observed by comparing the flat PCF and one groove cases for EQ2 in fig. 5. At at the bottom wall and at the top wall, the critical layer for the one groove case is much closer to the wall than it is for flat PCF. Consequently, there is a significantly smaller high speed region at either locations in the one groove case. In the two and three grooves cases, although the critical layer comes very close to the tips of the grooves, multiple high speed regions exist within the grooves.

We have observed two phenomena that affect kinetic energy density: (i) increasing wetted surface area leading to increased energy density, and (ii) increased proximity of the critical layer to the walls leading to reduced energy density. Laminar flows in grooved PCF are dominated by the former effect since the critical layer is far from the walls. The one groove case seems to be dominated by the latter effect as the vortex-streak structure penetrates into the grooves. In the two and three grooves cases, the latter effect is not very pronounced as the structure can only reach up to the tips of the grooves, and the increased wetted surface area leads to an overall increase in energy density.

We also note the increase in kinetic energy density for laminar flows to be significantly smaller than that in EQ1 and EQ2 (for two and three grooves per box). This is because the increase in kinetic energy density is normalized by the kinetic energy density of flat PCF, and the kinetic energy density is largest for the laminar case.

The trends in the bulk dissipation rate can also be explained in terms of the two competing phenomena. Their actions are as follows: 1) Increasing wetted surface area leads to an increase in dissipation rate, and 2) Increasing the proximity of the critical layer to the wall increases the dissipation rate. In laminar flow, as with energy density, the trend in dissipation rate is also dominated by the wetted surface area.

For EQ2, the critical layer is very close the walls. We expect the trends in dissipation rate to be dominated by the location of the critical layer. For the one groove per box case, increasing amplitude brings the critical layer closer to the walls, leading to increased dissipation rate (as well as power input and drag coefficient). For the three grooves case, increasing amplitude pushes the critical layer outside the tip of the grooves, leading to reduced dissipation rate. For the two grooves case, with increasing amplitude, the dissipation rate first increases and then decreases. However, the dissipation rate for the two grooves case remains lower than that for the one groove case at all amplitudes. This suggests that for the two grooves case, at low amplitudes, the critical layer might penetrate into the grooves.

For EQ1, the critical layer is far from the walls. The one groove case sees increasing dissipation rate with increasing amplitude as expected. The two and three grooves cases also see an increase in dissipation rate with increasing amplitude, suggesting that the squeezing out of the critical layer from the grooves is either insignificant or absent.

This mechanism for “drag reduction” is consistent with the one discussed in the literature (for instance, see fig. 7 in Lee & Lee (2001)). A Reynolds number dependent optimal spacing of the riblets exists for drag reduction; however, drag reduction is seen for riblet spacing smaller than 25 wall units (Dean & Bhushan (2010)). The near-wall cycle of turbulence also has a spanwise scale of wall-units (for example, Jiménez & Moin (1991); Hamilton et al. (1995)). This suggests that the greatest drag reduction involves riblets whose wavelength is a quarter of the characteristic size of the vortex streak structure.

The present results show a drag decreasing tendency (for EQ2, which is representative of the near-wall cycle of turbulence) for grooves whose wavelength is a half and a third of the characteristic size of the structure. And indeed, the drag reduction is greater when the riblet spacing is a third of the characteristic size than when it is a half. The present solver, outlined in appendix A, cannot solve for cases with smaller wavelengths owing inadequacy in spanwise resolution (limited to 32 grid points), but we can reason towards the observed optimal spacing by considering the two competing phenomena described earlier. Decreasing the wavelength (at a fixed amplitude) produces greater drag due to the increase in wetted surface area. So, the decrease in wavelength is only beneficial until the reduction in drag due to the squeezing out of the critical layer offsets the increase in drag due to the increase in wetted surface area. Figure 5 shows that the critical layer is already outside the grooves for the three grooves case (i.e., when the wavelength is a third of the characteristic size). Presumably, when the wavelength is reduced beyond a fourth or a fifth of the characteristic size, further decrease in wavelength would not have any significant effect on the critical layer. Beyond this point, reducing the wavelength would only see an increase in drag due to the increase wetted surface area.

We end this discussion by looking at the groove spacing (wavelength) in terms of inner units, , where is the wavelength of the grooves (non-dimensionalized by channel half-height) and is the friction Reynolds number. For the one groove per box case, , and for the two and three grooves cases, the groove spacing is and respectively. The friction Reynolds number may be expressed explicitly in terms of the bulk dissipation rate as follows.

(3.0) |

where is the dimensinoal friction velocity, is the speed of either walls, is the average (over plan area) of non-dimensional strain-rate at the bottom wall, is the average (over plan area) of non-dimensional strain-rate at the top wall, and and are the power input and bulk dissipation rate respectively.

For EQ2, the dissipation rate for the different cases lie roughly between and , with the flat wall case having . Using a nominal dissipation rate of , and , the friction Reynolds number is . Hence, in terms of wall-units, the spacing of the grooves is for the one groove case, for the two groove case, and for the three groove case. These spacings seem to be inconsistent with the observations in the literature, especially the drag increase at inner units. This apparent disparity is easily reconciled by taking into account the fact that the characteristic spanwise size of the vortex-streak structure itself is inner units, while the characteristic spanwise size of the near-wall cycle in turbulence is about inner units.

## 4 Conclusion

Equilibria (EQ1 and EQ2) have been successfully continued for PCF from the flat wall case to grooved walls. Discrete symmetries of these equilibria, , , and , have been exploited to allow computations on a 32x35x32 grid (streamwise, wall-normal, spanwise) for geometries with one and three grooves per spanwise period of channel half-heights. The geometry with two grooves per spanwise period (or any even number of grooves per period) does not admit and , and this case has been run at a lower resolution of 20x35x32 grid points. Geometries with more than three grooves per periodic box have not been investigated since the resolution with each groove becomes worse with increasing number of grooves.

The spanwise inhomogeneity of grooved PCF makes the relative spanwise positioning of solutions important. Computations initiated with iterates of non-zero phase fail to convergence, suggesting that such solutions may not even exist. This led us to speculate that symmetric solutions continued to symmetric geometries tend to get spatially anchored about locations that preserve the discrete symmetries. The fate of symmetric solutions when continued to asymmetric geometries remains an open question.

Grooved PCF with groove-wavelength of a half or a third of the vortex-streak structure being captured have been found to exhibit a drag reducing tendency for upper branch equilibria (EQ2 in the present work). This drag reduction is a result of the squeezing out of the critical layer of the vortex-streak structure by the grooves. If the wavelength is several times lower than this, we expect the drag to increase as the above effect saturates and the wetted surface area becomes the dominating factor. Both phenomena serve to increase the kinetic energy density.

Using the spacing of streamwise streaks in the near-wall cycle ( wall-units) as the relevant characteristic size, the present results are consistent with the observations of an optimal streak spacing wall units. This provides a strong motivation to extend exact invariant solutions to geometries with greater number of grooves per periodic box to find an optimal size for the grooves. However, the limitation imposed by the present numerical method on the resolution stops us from undertaking such an investigation. We hope that future studies involving more efficient numerical scheme for continuation of exact solutions can shed more light on the Reynolds number dependence of optimal riblet spacing.

Acknowledgement.— The authors acknowledge the use of the IRIDIS High Performance Computing Facility, and associated support services at the University of Southampton, in the completion of this work.

## A Steady state solver

The steady-state NSE in the physical domain of coordinates are given in section 2.2. The physical domain is mapped to a flat PCF domain using the domain transformation of section 2.1:

The -derivatives of are expressed as Fourier series in ,

(A 0) | ||||

The partial derivatives in relate to partial derivatives in as

(A 0) |

The variables, , , , and are discretized in the transformed domain (Fourier spectral method in and , and Chebyshev collocation method in ). The spatial derivatives in the physical domain in the steady-state NSE are expressed in terms of derivatives in the transformed domain; the Fourier coefficients for each term in the governing equation are now presented. Fourier coefficient of for the mode is denoted or , and the -dependence is not explicitly shown.

Fourier coefficients for the pressure gradient terms are

where is a linear operator representing differentiation with respect to , and is defined in eq. A 0.

The Fourier coefficients of the diffusion term, involving the Laplacian of the velocity, can similarly be expressed in terms of Fourier coefficients of the velocity field. The diffusion term is ; the Fourier coefficients of the Laplacian of streamwise velocity are

(A 0) | ||||

using Einstein summation notation for , , with . The Laplacians of the wall-normal and spanwise velocity involve similar expressions for their Fourier coefficients.

The convection term is fully non-linear; the velocity is not decomposed into a base-flow and fluctuations. The Fourier coefficients of the convection term in the streamwise momentum equation are expressed using sums of products of Fourier coefficients of the velocity field as

(A 0) |

Note the extra terms in the triadic interactions which are restricted to modes of the form in flat PCF. The terms in the wall-normal and spanwise momentum equations have similar expressions.

The continuity equation is also projected onto the Fourier modes as

(A 0) |

In the transformed system, the walls are at . Dirichlet conditions are imposed on each Fourier coefficient for the velocities and the divergence as

(A 0) | |||||

However, the boundary conditions on are set as to allow for wall motion.

### a.1 Matrix equation and Iterative scheme

The discretized steady-state NSE are cast as a non-linear matrix equation

(A 0) |

where is the state-vector containing the Fourier coefficients of , , , and , is a matrix representing the action of the linear terms, represents the non-linear term (convection), and represents a forcing. The forcing term is used to impose boundary conditions (or mean pressure gradient if a channel flow is to be solved for).

The state-vector is built as

(A 0) |

where and are positive integers representing the streamwise and spanwise resolutions, and is used to represent the coefficients for a particular Fourier mode. The entries of and other components extend from to , including the boundaries. Each Fourier coefficient () has entries, and equations are solved. of these equations are the NSE (including continuity equation) at the internal nodes; each row of corresponds to one out of the four NSE imposed at one internal wall-normal location for one Fourier mode. The remaining 8 rows are used for the boundary conditions of eq. A 0. To impose boundary conditions, all entries of corresponding to the 8 boundary rows are set to zero, and the entries of are set to reflect the terms on the left in each of the boundary conditions in eq. A 0. The terms on the right in the boundary conditions in eq. A 0, either 0 or 1, are set in the forcing term, .

#### a.1.1 Newton-Raphson method

The matrix equation is solved using the exact Newton-Raphson method, which iteratively refines a given estimate of an exact solution as

(A 0) | ||||

until a convergence criterion is satisfied. Here, denotes the Jacobian of , and the state-dependent matrix is the Jacobian of only the non-linear term ; can be shown to relate to the non-linear term as

(A 0) |

The superscript on is the iteration number, while subscripts in eq. A 0 represent Fourier coefficients of the fields at a particular iteration.

Preliminary numerical experiments with computing laminar solutions for grooved PCF and channel flows showed no improvement in convergence rate using parametric continuation. Hence, the flat PCF solution, either EQ1 or EQ2, is used as the initial estimate for all grooved PCF cases; the transformed domain of coordinates where the variables are discretized is identical to the flat PCF domain.

The Jacobian is formed (i.e. stored in memory) and inverted using an SVD-based inversion routine (lstsq of SciPy, which wraps around LAPACK’s gelsd). This iterative scheme is much simpler, and specific, compared to that of Gibson (2014) and Viswanath (2007), whose state equation is based on a time- map instead of the steady-state NSE; in their case, the Jacobian of is only expressed in terms of its action as a matrix-vector product, and a matrix-free variant of the Newton-Raphson method involving a constrained hookstep is used. The resolutions used in the present work are limited due to the explicit storage and inversion of the Jacobian; nonetheless, it allows simple computations that can be verified by computing the residual norm of the solutions on finer grids, which is computationally cheap.

To improve convergence of the Newton iterations, a line search method is used at the end of each iteration. If is the solution to the matrix equation, then an optimal correction, along , is found using a binary search over as

(A 0) |

The convergence criterion is that the residual norm,

(A 0) |

goes below a prescribed tolerance. The tolerance was usually set to , but this could not be reached in some cases (two groove per box). To ensure that the steady solutions are grid-independent, an spatial accuracy () is defined as the residual obtained when the state-vector is interpolated to a grid of twice the number of nodes(/modes) in each coordinate direction.

The divergence-free condition is used directly as an equation in , and the iterations converge without issue. Since the algorithm does not involve time-marching, we found no particular reason to use a pressure Poisson equation instead of the divergence equation for velocity.

- Bannier et al. (2016) Bannier, A., Garnier, E. & Sagaut, P. 2016 Riblets induced drag reduction on a spatially developing turbulent boundary layer. In Progress in Wall Turbulence 2, pp. 213–224. Springer.
- Blackburn et al. (2013) Blackburn, H. M., Hall, P. & Sherwin, S. J. 2013 Lower branch equilibria in couette flow: the emergence of canonical states for arbitrary shear flows. Journal of Fluid Mechanics 726, R2.
- Chantry & Kerswell (2015) Chantry, M. & Kerswell, R. R. 2015 Localization in a spanwise-extended model of plane couette flow. Physical Review E 91 (4), 043005.
- Dean & Bhushan (2010) Dean, B. & Bhushan, B. 2010 Shark-skin surfaces for fluid-drag reduction in turbulent flow: a review. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 368 (1929), 4775–4806.
- Eckhardt et al. (2008) Eckhardt, B., Faisst, H., Schmiegel, A. & Schneider, T. M. 2008 Dynamical systems and the transition to turbulence in linearly stable shear flows. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences 366 (1868), 1297–1315.
- Flack & Schultz (2014) Flack, K. A. & Schultz, M. P. 2014 Roughness effects on wall-bounded turbulent flows a). Physics of Fluids (1994-present) 26 (10), 101305.
- Gibson & Cvitanović (2008) Gibson, J. F.and Halcrow, J. & Cvitanović, P. 2008 Visualizing the geometry of state space in plane couette flow. Journal of Fluid Mechanics 611, 107–130.
- Gibson (2014) Gibson, J. F. 2014 Channelflow: A spectral Navier-Stokes simulator in C++. Tech. Rep.. U. New Hampshire, Channelflow.org.
- Gibson & Brand (2014) Gibson, J. F. & Brand, E. 2014 Spanwise-localized solutions of planar shear flows. Journal of Fluid Mechanics 745, 25–61.
- Gibson et al. (2009) Gibson, John F, Halcrow, Jonathan & Cvitanović, Predrag 2009 Equilibrium and travelling-wave solutions of plane couette flow. Journal of Fluid Mechanics 638, 243–266.
- Halcrow (2008) Halcrow, Jonathan J 2008 Charting the state space of plane couette flow: equilibria, relative equilibria, and heteroclinic connections. PhD thesis, Georgia Institute of Technology.
- Hall & Sherwin (2010) Hall, P. & Sherwin, S. 2010 Streamwise vortices in shear flows: harbingers of transition and the skeleton of coherent structures. Journal of Fluid Mechanics 661, 178–205.
- Hall & Smith (1991) Hall, P. & Smith, F. T. 1991 On strongly nonlinear vortex/wave interactions in boundary-layer transition. Journal of Fluid Mechanics 227, 641–666.
- Hamilton et al. (1995) Hamilton, J. M., Kim, J. & Waleffe, F. 1995 Regeneration mechanisms of near-wall turbulence structures. Journal of Fluid Mechanics 287, 317–348.
- Jiménez (2004) Jiménez, J. 2004 Turbulent flows over rough walls. Annual Review of Fluid Mechanics 36, 173–196.
- Jiménez & Moin (1991) Jiménez, J. & Moin, P. 1991 The minimal flow unit in near-wall turbulence. Journal of Fluid Mechanics 225, 213–240.
- Kasliwal et al. (2012) Kasliwal, A., Duncan, S. & Papachristodoulou, A. 2012 Modelling channel flow over riblets: Calculating the energy amplification. In Control (CONTROL), 2012 UKACC International Conference on, pp. 625–630. IEEE.
- Kawahara et al. (2012) Kawahara, G., Uhlmann, M. & van Veen, L. 2012 The significance of simple invariant solutions in turbulent flows. Annual Review of Fluid Mechanics 44, 203–225.
- Lee & Lee (2001) Lee, S-J & Lee, S-H 2001 Flow field analysis of a turbulent boundary layer over a riblet surface. Experiments in Fluids 30 (2), 153–166.
- Moradi & Floryan (2013) Moradi, H. V. & Floryan, J. M. 2013 Maximization of heat transfer across micro-channels. International Journal of Heat and Mass Transfer 66, 517–530.
- Nagata (1990) Nagata, M. 1990 Three-dimensional finite-amplitude solutions in plane couette flow: bifurcation from infinity. Journal of Fluid Mechanics 217, 519–527.
- Park & Graham (2015) Park, J. S. & Graham, M. D. 2015 Exact coherent states and connections to turbulent dynamics in minimal channel flow. Journal of Fluid Mechanics 782, 430–454.
- Perry & Marušic (1995) Perry, A. E. & Marušic, I. 1995 A wall-wake model for the turbulence structure of boundary layers. part 1. extension of the attached eddy hypothesis. Journal of Fluid Mechanics 298, 361–388.
- Squire et al. (2016) Squire, D. T., Morrill-Winter, C., Hutchins, N., Schultz, M. P., Klewicki, J. C. & Marusic, I. 2016 Comparison of turbulent boundary layers over smooth and rough surfaces up to high reynolds numbers. Journal of Fluid Mechanics 795, 210–240.
- Townsend (1980) Townsend, A. A. 1980 The structure of turbulent shear flow. Cambridge university press.
- Viswanath (2007) Viswanath, D. 2007 Recurrent motions within plane couette turbulence. Journal of Fluid Mechanics 580, 339–358.
- Waleffe (1997) Waleffe, F. 1997 On a self-sustaining process in shear flows. Physics of Fluids (1994-present) 9 (4), 883–900.
- Walsh (1990) Walsh, M. J. 1990 Effect of detailed surface geometry on riblet drag reduction performance. Journal of Aircraft 27 (6), 572–573.
- Willis et al. (2013) Willis, A. P., Cvitanović, P. & Avila, M. 2013 Revealing the state space of turbulent pipe flow by symmetry reduction. Journal of Fluid Mechanics 721, 514–540.
- Willis et al. (2017) Willis, A. P., Short, K. Y., Burak Budanur, N., Farazmand, M. & Cvitanovic, P. 2017 Relative periodic orbits form the backbone of turbulent pipe flow. arXiv preprint arXiv:1705.03720 .